Revisiting transverse momentum broadening in dense QCD media

We reconsider the problem of transverse momentum broadening of a highly-energetic parton suffering multiple scatterings in dense colored media, such as the thermal Quark-Gluon plasma or large nuclei. In the framework of Moli\`ere's theory of multiple scattering we re-derive a simple analytic formula, to be used in jet quenching phenomenology, that accounts for both the multiple soft and hard Rutherford scattering regimes. Further, we discuss the sensitivity of momentum broadening to modeling of the non-perturbative infrared sector by presenting a detailed analytic and numerical comparison between the two widely used models in phenomenology: the Hard Thermal Loop and the Gyulassy-Wang potentials. We show that for the relevant values of the parameters the non-universal, model dependent contributions are negligible, at LHC, RHIC and EIC energies thus consolidating the predictive power of jet quenching theory.


I. INTRODUCTION
The propagation of highly energetic partons in a dense QCD medium is affected by multiple collisions with medium constituents that cause their transverse momentum distribution to broaden. In addition, medium-induced gluon emissions can be triggered by these random kicks leading to the formation of a parton cascade.
The latter is the main mechanism of parton energy loss in large QCD media [1,2] and its theoretical description has drawn a lot of attention in the last two decades. However, most analytic results used in phenomenology were obtained in certain kinematic limits [3][4][5][6][7][8] and a solid and quantitative understanding of in-medium jet modifications, as measured in experiment, requires further theoretical control over the full kinematic range of the transverse momentum dynamics. This is currently of particular importance as the field moves towards highly precise comparisons between data and theoretical predictions.
In the present manuscript we study the single particle momentum broadening distribution for the case of an energetic parton propagating through a homogeneous plasma brick of length L. A complementary study of the medium induced radiation mechanism, following a similar scheme to the one present in this paper, is discussed in [9,10] for the transverse momentum dependent gluon spectrum and in [11][12][13] for its integrated version.
Although substantial effort has been put into studying momentum broadening in increasingly more real-istic and complex scenarios [14][15][16][17][18][19][20][21][22][23][24], theoretical uncertainties in the simplistic setup considered here remain to be understood. The origin of these ambiguities is mainly associated to the infrared (IR) modeling of the parton-medium interaction. In addition, it is common in phenomenological studies to take the broadening distribution in some limiting form such as the Gaussian approximation that applies to the multiple soft scattering regime. Even though the latter does not capture the correct physics in the whole range of transverse momenta, its simple analytical form makes it suitable for theoretical calculations, e.g. [25,26], together with a straightforward implementation in jet quenching Monte-Carlo event generators, where a given momentum broadening probability has to be sampled [27][28][29].
The goal of this paper is twofold. First, we revisit the theory of multiple scattering by Molière [30,31] which provides a more accurate description of the broadening distribution over the full range in transverse momentum. In particular, the Gaussian behavior and the power-law Coulomb tail, pertaining to Rutherford scattering, are reproduced at small and large transverse momentum, respectively. On the other hand, we systematically study the role of non-perturbative modeling on the momentum broadening probability distribution. We show that, for realistic values of the parameters, resulting non-universal power corrections are negligible. This paper is structured as follows. Section II introduces the kinetic theory formulation of the momentum broadening probability distribution, along with the two medium models to be explored. Next, Section III reviews Molière's systematic approach to multiple scattering theory applied to momentum broadening. In Section IV we explore the dependence of the broadening distribution on IR modeling. In addition, we discuss our findings in the parameter space explored by LHC, RHIC and the upcoming EIC in Section V. We summa-rize our results and discuss possible future directions in Section VI. Complementary material is presented in Appendices A to D.

II. KINETIC DESCRIPTION OF MOMENTUM BROADENING
The probability for a highly energetic parton traversing a dense QCD medium to acquire a transverse momentum k from multiple scattering during a time L is denoted by P(k, L). Its time evolution can be formulated in kinetic theory and is given by [32,33] where γ(q) is the collision rate and is related to the inmedium elastic scattering cross-section to be discussed shortly. C R is the parton color factor in representation R, i.e., C R = C F = (N 2 c − 1)/(2N c ) and C R = C A = N c for a quark and a gluon, respectively. Here and throughout we used the shorthand notation q ≡ d 2 q/(2π) 2 for transverse momentum integration, and x ≡ d 2 x for transverse coordinate space integrals.
Physically, Eq. (1) has a simple interpretation: in an infinitesimal time step δt, the probability for a parton to end up with momentum k at L + δt equals the probability of starting with momentum k − q and acquiring momentum q during δt. In addition, one must subtract the probability of already starting with momentum k and diffusing to some other momentum mode. This structure ensures the unitarity of P and, consequently, its probabilistic interpretation. We remind the reader that the previous result is only valid in the strict eikonal limit, where the acquired transverse momentum is much smaller than the parton energy.
Turning to the collision rate γ(q), a closed expression exists in two distinct scenarios. In the case of a medium formed by static scattering centres with Yukawa-type interactions, γ(q) is given by the so-called Gyulassy-Wang (GW) model [34] where n represents the density of scattering centres and µ is the (GW) screening mass. When the QCD medium is in thermal equilibrium, γ(q) can be computed with Hard Thermal Loop theory (HTL) resulting in [35] γ HTL (q) = where g is the strong coupling and m D corresponds to the Debye mass whose temperature (T ) dependence is given, at leading order in g, by m 2 where n f is the number of active light flavours. While the UV behavior is common to both models, i.e., q −4 , the main difference is in the way the IR Coulomb divergence is regulated. In the limit q → 0 we have, γ GW (q) → const. while γ HTL (q) → q −2 . The possibility of mapping µ and m D together with other differences and similarities between these two models for γ(q), commonly used in phenomenology, will be discussed in Section IV. Defining S(x, L) as the Fourier transform of P(k, L), Eq.
(1) becomes local in position space and can be easily solved to give where in the last step we assumed that the medium is a homogeneous brick of length L and defined the socalled dipole cross-section [36] as Explicit formulas for v(x) in the GW and HTL models can be obtained using the integrals computed in Appendix B. They read and v HTL (x) = 2q 0 m 2 is the bare jet quenching transport parameter and α s = g 2 /(4π). Finally, by Fourier transforming Eq. (6) back into momentum space, the probability for a parton to acquire transverse momentum k reads where we subtract the no-broadening contribution that would result in a δ (2) (k) term that does not contribute at finite k (see [37]). This term only vanishes when v(∞) = +∞. From Eq. (8) and Eq. (9), we observe that neither GW nor HTL satisfy this condition, but rather saturate at large |x|. Regarding the normalization of the probability distribution given by Eq. (11), it is easy to verify, using v(0) = 0, that k P(k, L) = 1 − exp(−v(∞)L).
To the best of our knowledge, a closed analytic form for Eq. (11) when using either GW or HTL expressions for v(x) does not exist. It is the goal of the next Section to achieve an analytic representation of the transverse momentum probability distribution, P(k, L), by combining an expansion of Eq. (11) in universal and model dependent terms together with Molière's theory [30,31] of multiple scattering.

III. P(k)-DISTRIBUTION FROM MOLIÈRE'S THEORY OF MULTIPLE SCATTERING
Adopting the small dipole size approximation, Eq. (8) (Eq. (9)) can be expanded to linear order in µ 2 x 2 (m 2 D x 2 ), which we refer to as the leading power (LP) contribution to S(x). The choice of this terminology is motivated by the desire of distinguishing between universal power corrections of the form (Q 2 s0 x 2 ) n (with Q 2 s0 ≡q 0 L [13], see Appendix B) and non-perturbative power corrections of the form (µ 2 x 2 ) n ((m 2 D x 2 ) n ), which we refer to as next-to-leading power terms. To leading power (LP), we obtain v GW (x) =q and v HTL (x) =q We see that, at this level of accuracy, it is possible to relate the GW and HTL parameters by defining the following physical scale [38], More concretely, this approximation leads to Given the physical scale defined in Eq. (14), this contribution is model independent. The possibility of expanding Eq. (8) and Eq. (9) up to next-to-leading power order is explored in Section IV. To obtain an analytic form of P(k, L) from Eq. (15) additional assumptions concerning the parton's energy have to be adopted.
At very high energy, i.e. k 2 Q 2 s0 , the dipole's transverse size is given by |x| ∼ 1/|k| 1/Q s0 . Then, S LP (x) can be expanded to linear order The zeroth order term can be neglected as it does not contribute to the P(k, L) distribution. Thus, S LP (x) is proportional to the dipole cross-section, i.e. it is dominated by single hard (SH) scattering contributions. In this case, P SH (k, L) reads where used the fact that for k 2 The momentum broadening probability distribution given by Eq. (17) captures the expected Coulomb-like 1/k 4 behavior at high momentum transfers. The scale at which multiple scattering (MS) becomes important is encoded in the medium opacity parameter χ ∼ Q 2 0s /µ 2 * ∼ L/ mfp , where mfp is the in-medium mean free path. When χ 1, the medium is dilute and therefore single (rare) hard scattering events, as discussed above, dominate the contribution to P(k, L). Conversely, when χ 1, the medium is densely populated and multiple soft scatterings become the relevant mechanism for momentum transfer.
In the k 2 Q 2 s0 regime, the logarithm in Eq. (15) is slowly varying with x and can be regulated by a large momentum scale Q 2 ∼ Q 2 s0 , so that taking into account all orders in Q 2 s ∼ Q 2 s0 log(Q 2 /µ 2 ) one obtains a Gaussian representation of momentum broadening where the super-script refers to multiple scattering (MS). As mentioned in the introduction, Eq. (19) is the preferred option in widely used and successful jet quenching Monte Carlo event generators such as the Hybrid model [28] or the newly developed code by the Saclay group [29]. Despite the widespread phenomenological application of Eq. (19), it fails to accurately describe the hard 1/k 4 tail of P(k, L), thus missing the physics associated with single hard scattering events (see Eq. (17)) [39]. The importance of such contributions has been recently studied numerically in [40,41].
In this work, we propose to use an efficient expansion scheme developed by Molière in 1948 [30] in order to provide a simple, analytic formula that encodes the correct behavior of P(k, L) from small to large k. Molière's approach is based on expanding v(x) around the multiple soft scattering approximation v MS (x). This is achieved by splitting the Coulomb logarithm into two pieces: a large, but constant, logarithm and a small, x 2dependent term which is treated perturbatively. Using this scheme, the leading power (LP) dipole cross-section can be written as where Q 2 is known as matching scale, v MS (x) corresponds to the cross-section entering Eq. (19) and δv(x) can be considered a perturbation as long as Q 2 µ 2 * , such that log 1 . This decomposition leads to the following definitions of the relevant scales in the problem where Q 2 can be taken to be proportional to Q 2 s , i.e., Here a is a free parameter to be determined for each set of medium parameters [42]. Given a value of a, by inserting Eq. (22) in Eq. (21), one obtains the following transcendental equation where the choice a = 1 corresponds to Moilère's prescription [30]. We also define the effective transport coefficientq asq Using Eq. (20), one can expand Eq. (11) (at LP accuracy) around the MS solution in powers of δv(x) as where we integrate over x and sum over n, from n = 0 to infinity. Notice that we still have k P LP (k, L) = 1, since all terms in Eq. (25), apart from the (0) term, vanish after integrating over k and x. Formally, the series representation introduced above is asymptotically divergent. Nevertheless, a very good approximation of the exact solution can be obtained when the series is truncated before it diverges at n < n max ∼ Q 2 s0 /µ 2 * . For a broader discussion on the origin of this truncation, we refer the reader to [43], where Molière's expansion was explored in the Color Glass Condensate framework.
In order to recast Eq. (25) in a more compact form we define the dimensionless expansion parameter This allows us to re-write Eq. (25) as while the next-to-leading order (1) term yields [30] λf (1) where I 1 (x) and I 2 (x, a) are given in Appendix A and Ei is the exponential integral function [44]. Here we introduced the reduced Laplacian operator ∆ Combining the results for f (0) and f (1) , we obtain one of the main results of this paper (originally derived by Molière [30]): Let us now verify that the above result reproduces the expected asymptotic behavior at small and large transverse momentum. The (0) contribution matches exactly the MS solution from Eq. (19), and thus its limiting behavior is easy to analyse. At large momentum transfers, k 2 Q 2 s , s it becomes independent of k. More importantly, the Gaussian profile implies that the typical momentum transverse acquired due to momentum broadening k 2 typ ∼ Q 2 s . Therefore, P (0) (k, L) correctly captures the physics associated to multiple soft scattering at scales s . For the (1) term, we use two limiting forms of the Ei function. That is, when x → ∞, Ei(x) ≈ e x (1/x + 1/x 2 + 2/x 3 ) so that the large k behavior of the (1) term reads This result matches Eq. (17) and therefore the (1) term successfully encodes the hard 1/k 4 tail of the full P(k, L) distribution. As a consequence, it is physically preferable for phenomenological applications to use Eq. (30) instead of Eq. (19). On the other end, x → 0, Ei(x) ≈ γ E + log x and then which, up to a small constant logarithm, corresponds to the MS result (Eq. (19)) in this kinematic limit, suppressed by a power of λ. This is analogous to the small energy limit behavior obtained in [11][12][13] for the gluon emission spectrum. In Fig. 1, we evaluate Eq. (27) at (0), (1) order and their sum (0)+(1). These curves are compared to the exact numerical solution of Eq. (11) when plugging the GW dipole cross-section given by Eq. (8). A small value of the expansion parameter λ is chosen on purpose such that this figure represents a proof-of-concept of the proposed scheme in its regime of validity. In the multiple scattering regime, i.e. at small-k ⊥ , the (0) contribution dominates over the (1) term as expected from our asymptotic analysis. Nevertheless, the (0)+(1) curve shows a small discrepancy, to be quantified in what follows, with respect to the full GW result. The situation is improved at large-k ⊥ , where the (1) contribution correctly captures the power-law tail completely absent in the (0) scenario. This figure demonstrates how a purely analytic, two terms expansion given by Eq. (30) exhibits an excellent agreement with the numerically obtained P(k, L) using GW/HTL models for v(x).
The natural question arises as to what is the value of the λ-parameter at which the expansion fails to reproduce the GW result. This problem, together with the role of higher orders, is tackled in Fig. 2  is degraded when increasing λ both at low and large k. This result is expected as the larger λ gets, the less precise is to consider δv(x) as perturbative contribution in Eq. (20). The relevant values of λ for current and future colliders will be discussed in Section V. As shown in the bottom panel of Fig. 2 this discrepancy can be alleviated by adding extra terms in the expansion. In particular, when adding the (2) (see Eq. (25)) contribution for λ = 0.1 we find a ratio to the exact GW result close to one in the whole interval in k. Unfortunately, we were not able to find yet a general analytic expression for the n-th term in the series, thus higher orders have to be computed numerically.

IV. SENSITIVITY TO IR MODELING UP TO NEXT-TO-LEADING POWER ORDER REMOVE[VIA A TWIST EXPANSION]: GW AND HTL COMPARISON
In the previous Section, as a first step towards an analytic expression for P(k, L) that encompasses the main physical mechanisms, we have expanded the dipole cross-section to leading power order (see Eq. (15)). In order to assess the sensitivity of transverse momentum broadening to the non-perturbative infrared structure of a given model for γ(q), going beyond the leading power (LP) term is mandatory. In what follows, we fix the hard scale of the problem Q s0 such that we are only sensitive to the dependence of this expansion with respect to the infrared regulator µ * . Consequently, these additional terms are expected to modify the low momentum regime of P(k, L). Hence, this will be the explored region in the figures of this section. At this point, we would like to emphasize that the expansion in universal and non-universal terms is intrinsically different from the Molière one introduced in the previous section. In short, the former explores the infrared sector through non-universal contributions, while the latter is a perturbative expansion with model independent terms. A systematic study of these non-universal contributions is so far missing in the literature. In turn, GW and HTL models are typically treated as two independent descriptions of the QCD medium, thus ignoring the fact that they can be mapped onto one another at LP accuracy. The common practice in jet quenching phenomenology of treating these models as if they were unrelated can be problematic given that: i) in the absence of a map between the different IR regulators to the physical Debye mass, any comparison between results assuming different models is meaningless, ii) a quantitative and controlled understanding of the role of non-perturbative physics at the infrared scale cannot be reached.
The importance of mapping the IR scales involved in GW and HTL becomes apparent in Fig. 3 where the value of m D is fixed and the corresponding value for µ, to be plugged in Eq. (8), is obtained using Eq. (14). The momentum broadening probability distribution computed with these two medium models only shows significant differences at small-k ⊥ for values of m 2 D larger than ∼ 1 GeV 2 . Hence, once the appropriate matching between m D and µ is considered, the use of these two in-medium elastic scattering cross-sections leads to discrepancies ≤ 10% in the infrared sector [45]. This result is in agreement with previous calculations where the mapping was used [13,46].
To characterize the sensitivity of transverse momen- tum broadening to infrared physics, we proceed to expand v(x) as given by Eq. (8) (Eq. (9)) up to next-toleading power (NLP) order, such that v(x) takes the form where c 1 and c 2 are model dependent constants given in Appendix C for the GW and HTL models. We would like to point out that Eq. (33) is an explicit manifestation of the fact that NLP corrections are non-universal, since there is no self consistent way of mapping different models to some functional form representing what would be measured in experiment. In fact, even using the LP map introduced in Eq. (14), we clearly observe that universality is not regained (i.e. c 1 and c 2 are truly model dependent). Nonetheless, using such a map allows us to directly gauge the order of magnitude of the model dependent contributions to P(k, L). The implications of considering the LP or NLP approximations to GW/HTL potentials instead of the full result are shown in Fig. 4. Regarding the leading power term, it fails to reproduce the full result accurately for relatively small values of m D , similar to those explored in current colliders as will be discussed in the next Section. Two comments are in order concerning the next-toleading power term. First, focusing on the HTL case, its magnitude is remarkably small. This is an important aspect as it indicates a close to minimal sensitivity to nonperturbative assumptions. Second, when turning to the GW case this difference with respect to the LP is significantly more pronounced. The underlying reason is that the map between m D and µ through the physical mass µ * is valid only at LP order. Therefore, when extending its applicability to NLP, more severe deviations between GW and HTL are deemed to occur.
Further, it is possible to combine the next-toleading power expansion of the dipole cross-section (see Eq. (33)) with Molière's scheme, by shifting the expan- and continue treating δv(x) (see Eq. (20)) as a perturbation around it. In this way, the λ-expansion enables the description of P for all k, while it is possible to explore the dependence on IR modeling via the expansion beyond the LP contribution.
The first non-trivial term in Eq. (25) is thus promoted (to linear order in µ 2 * [47]) to the following simple expression, Equation (34) is no longer a function of solely Q 2 /µ 2 * (i.e. λ), but also depends on the IR regulator alone and on the constants c 1 and c 2 ; this is an explicit example of the non-universality attribute we have associated to these terms.
The role of the NLP contribution in Molière's expansion at first (1) order is displayed in Fig. 5 for both GW (top) and HTL (bottom) potentials. In all cases, the inclusion of the NLP term has an effect smaller than 5% being this contribution larger when increasing m D , as expected. This fact confirms the mild infrared dependence in our description of P(k, L). Further, the NLP contribution enhances the value of P(k, L) at small-k ⊥ thus reducing the discrepancy with the full GW/HTL result observed in Figs. 1 and 2.

V. RELEVANCE FOR PHENOMENOLOGY AT CURRENT AND FUTURE COLLIDERS
Up to now, we have performed a theory guided selection of the parameters involved in the description of the medium in order to highlight the relevant region of interest. That is, in Section III we have deliberately chosen small values of λ as required by Molière's scheme or large values of m D in Section IV to better illustrate the contribution of next-to-leading power terms.
In what follows, we turn our attention to more realistic scenarios such as the ones being currently explored by LHC and RHIC together with the upcoming EIC. The relevant parameters for these three colliders are provided in Table I. For a given value of the medium length, L, and its temperature, all parameters in this table are uniquely determined. More concretely, m D is obtained through the leading order HTL formula mentioned in Section II and Q 2 s is given by Eq. (23) wherê q 0 = 18πα 2 s T 3 , where we take α s = 1/π. In all three cases, we take the length of the medium to be L = 5 fm roughly corresponding to the radius of both Pb and Au nuclei. For RHIC and LHC, we use the temperature estimates T RHIC = 220 MeV and T LHC = 300 MeV, corresponding to an indicative average temperature for each setup in accordance with the values found in e.g. [48][49][50]. We also assume the medium can be described by the HTL model. Turning to the EIC case, the characteristic scale of the medium explored by the partonic probe is no longer the temperature, but the nucleon size that can be approximated as 1/Λ QCD . This provides a natural value for the infra-red regulator in the GW model, i.e. µ = Λ QCD = 200 MeV. Regardingq 0 L, we make use of the Color Glass Condensate effective field theory to estimate its value following [51], leading tô q 0 L = 0.35 GeV 2 . Note that this is in the ballpark of the value used in [52]. We consider in what follows quark jets with C R = C F = 4 3 . The results for Pb+Pb collisions at LHC energies are shown in Fig. 6. Regarding the large-k ⊥ region, the inability of the leading order term to capture the behavior of both GW and HTL models becomes apparent in the top panel. In contrast, the (0)+(1) and (0)+(1)+(2) contributions deliver a distribution identical to the GW result up to 5% degree of accuracy for k ⊥ > 20 GeV. This statement also applies to HTL because the difference among these models is negligible in this region, as shown in the middle panel. In the infrared sector, k ⊥ < 1 GeV, we observe deviations of order 20% which can be improved with higher order terms, although this sector's contribution to physical observables is typically negligible. Around the peak of the distribution, the disagreement between (0)+(1)+(2) order with GW is of the order of 20%. This result could be improved by adding higher order terms in Molière's expansion (see Eq. (25)).
Turning to RHIC energies, similar features to those at LHC are observed at large transverse momentum as shown in Fig. 7. An important remark concerning the small-k ⊥ sector is a slight increment of the GW vs. HTL disparity with respect to the LHC case, despite probing a smaller m D , see Tab. I. The underlying reason is that the value of m D has to be sufficiently separated from the hard scale in the problem. That is, the larger the value of the ratio Q s0 /m D is, the smaller the sensitivity to the infrared, model-dependent contributions. Indeed, this ratio is 20% bigger at LHC than RHIC, thus explaining the mild differences between GW and HTL for these two experiments.
Finally, the future Electron-Ion Collider at Brookhaven National Lab [53] will open a new avenue to study modifications in jet observables when compared to hadronic colliders. The reason is that highly-energetic partons will not encounter a thermalized QGP at the EIC, but rather a dense gluonic system. Multiple interactions with this over-occupied gluon state naturally leads to transverse momentum broadening. Its probability distribution is displayed in Fig. 8, where we have assumed e+Au to be the collision system, and the prospected top energy for this machine and the GW model is employed. We would like to point out that, in this case, we obtain that the relevant hard scale is Q 2 s = 1.75 GeV 2 . At the same time, the average transverse momentum received from the cold medium is shifted towards smaller values when compared to RHIC and LHC. This experiment will probe the largest value of λ and, therefore, the application of Molière's scheme is less successful than for the LHC and RHIC setups. Nevertheless, the evolution of λ as one changes the experimental set-up is slow and therefore this method is suitable for semi-quantitative exploratory studies.
The role played by the non-perturbative contributions at RHIC, LHC and EIC energies its suppressed when compared to the deviations due to the (1) expansion, as was shown in the previous Section. Therefore, we refrain from providing the numerical results obtained when considering these contributions.  Fig. 6 but for the future EIC. Note that we do not consider the HTL potential as a thermal medium at the EIC is not expected.

VI. CONCLUSIONS AND OUTLOOK
The work presented in this manuscript follows the current (and future) global effort towards a more precise quantitative description of jet quenching effects [11-13, 46, 48]. In particular, we have: i) re-derived a description of the single particle momentum broadening distribution, first introduced by Molière, that is able to reproduce the Gaussian behavior at small-k ⊥ together with the power-law tail through a simple analytic expression, i.e. Eq. (30); ii) explored the impact of nonperturbative modeling of the dipole cross-section on the broadening distribution; iii) studied the applicability of Molière's scattering theory for the description of single particle momentum broadening at the current experimental conditions explored by LHC, RHIC and the foreseen EIC.
Although the major analytic result of this paper (i.e. Eq. (30)) has been known for over 70 years, its application to phenomenological studies in jet quenching has, to a large extent, been ignored. This is particularly surprising given the widespread usage of the multiple soft scattering parametrisation (i.e. Eq. (19)) that is known to fail to describe the hard sector, as we have shown both analytically and numerically for meaningful parameter selections. It should be noted that Molière's proposal, if restricted only to the first two leading terms whose expressions are known analytically, requires a pronounced hierarchy between the hard and soft scales of the problem in order to converge to the exact result when considering GW/HTL potentials. In particular, we found that such a separation is reasonably satisfied at current LHC and RHIC conditions, but is far from being achieved at the expected EIC setup.
Another important finding of this work is that nonperturbative contributions, intrinsically associated to the description of the infrared structure of in-medium scattering, are found to be quantitatively small. However, we (again) emphasize that to achieve a meaningful comparison between medium models, matching to a physical set of parameters at leading power (LP) accuracy is required. For the three physically interesting cases explored in Section V, non-perturbative contributions are highly suppressed compared to the ones obtained from Molière's scheme.
The present results are of relevance for jet quenching phenomenology and have, for example, been recently applied to computing the nuclear modification factor R AA [54]. Interesting future studies could focus on the impact of the single hard scattering tail on jet substructure observables such as the groomed opening angle [25] or the transverse momentum of the hardest splitting in the jet tree [55], and be further related to the presence of Molière scattering centers in the Quark Gluon Plasma [19,56].
As mentioned in Section III, Eq. (30) can be readily employed in any Monte Carlo event generator that includes jet quenching effects. In addition, the formalism described in this paper can also be applied to small-x physics [57]. A possibility would be to use Molière's scheme to construct the unintegrated gluon distribution needed for particle production in p+A [43,58].
Using the map given by Eq. (14), the model dependent constants appearing in the generalized NLP dipole cross-section (Eq. (33) In this Appendix we outline how some of the results presented in Section III can be obtained solely using the kinetic description of momentum broadening introduced in Section II. In addition, we provide an equivalent formulation of Molière's expansion in the form of a generalized diffusion equation in momentum space.
One can rewrite Eq. (1) using Eq. (7) as ∂ ∂L P(k, L) = − q v(q)P(k − q, L) . (D1) In the high energy limit, single hard scattering dominates, which implies that only the first iteration of Eq. (D1) contributes and thus one can write [32] P SH (k, which is satisfied by Eq. (17). On the other hand, when multiple soft scattering dominates, v(x) is quadratic in x, which allows one to write Eq. (D1) as a diffusion equation with diffusion parameterq ∂ ∂L P MS (k, L) =q 4 ∇ 2 k P MS (k, L) , (D3) which can be solved to yield Eq. (19).
Combining the kinetic description of momentum broadening with Molière's approach leads to an hierarchy of (trivially) coupled diffusion equations with a source term. To explicitly see this, we use Eq. (20) and Eq. (25)  Eq. (D4) can, in principle, be solved order by order using Green's method, but, to the best of our knowledge, and considering the aims of this paper, there is no obvious advantages over the procedure followed in the main text.