Jet broadening in dense inhomogeneous matter

In this work, we study the jet momentum broadening in an inhomogeneous dense QCD medium. The transverse profile of this nuclear matter is described within a gradient expansion, and we focus on the leading gradient contributions. The leading parton is allowed to interact multiple times with the background through the soft gluon exchanges. We derive the associated final particle distribution using both the GLV opacity series and the BDMPS-Z formalism. We further discuss the modified factorization of the initial and final state effects and its consequences for phenomenological applications in the context of heavy-ion collisions and deep inelastic scattering. Finally, we present the broadening probability (describing the final state effects) in several limiting regimes, and give its numerical estimates for phenomenologically motivated sets of parameters.


I. INTRODUCTION
Jets are collimated sprays of particles, produced by hadronization and branching of an energetic quark or gluon (parton), which are often found in the final state of experiments on ultrarelativistic particle collisions. If before the hadronization stage the parton cascade develops in the presence of an underlying medium, produced in the same collision, the jet substructure gets modified due to the interactions with the background. The simplest manifestation of this process is the suppression of jet energy by matter, commonly referred to as jet quenching, which has attracted significant attention in the literature [1], for a more recent review see [2][3][4][5]. Due to their high sensitivity to the spacetime structure of the medium, jets provide a promising tomographic tool to study the real-time evolution of nuclear matter both in heavy-ion collision (HIC) and deep inelastic scattering (DIS) experiments, see e.g. [6][7][8][9][10][11][12][13][14][15][16][17] and references therein.
The jet-medium interaction can be successfully described within perturbative QCD (pQCD) supplemented with a medium model, which is usually based on a collection of medium-induced stochastic color fields [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32]. In this picture, an energetic parton interacts with the matter through multiple t−channel gluon exchanges. Although such interactions lead to a negligible energy depletion of the hard parton, they induce gluon bremsstrahlung resulting in an energy loss. Description of these processes within pQCD is in general a complex problem, and some simplifying assumptions are usually needed. However, the commonly considered eikonal or static source approximations lead to a decoupling of the medium evolution from the jet energy loss and substructure modification observables, for a discussion see [15]. Thus, to extract the properties of the underlying medium evolution and information about its structure from jet observables one needs the associated theoretical framework to be extended beyond the simplest physical regimes.
Earlier attempts to include the effects of the medium flow into the jet energy loss calculations go back to [33][34][35], where the medium dilution was considered, and to [36,37], where the flow was described within a phenomenologically motivated model with an additional momentum transfer. Transverse flow effects could also be partially accounted for from purely kinematic arguments, see e.g. [38][39][40]. Only more recently, the medium evolution effects on the jet-medium interactions were formally included into the Gyulassy-Levai-Vitev (GLV) opacity expansion framework [15]. It was shown that the medium flow and variation of its properties in the transverse directions can be treated within the same medium model of stochastic color fields induced by in-medium sources, if the sources are allowed to move and the medium properties encoded in their potentials and density are changing from point to point. In the developed formalism, the changes in the local properties of the matter are described within a gradient expansion analogous to the one used in hydrodynamics, commonly applied to describe the evolution of the quark-gluon plasma (QGP) produced in HIC. In this way, [15] extends the idea to describe interactions of a probe with a hydrodynamically evolving matter within the same gradient expansion introduced in holographic models for strongly interacting plasmas, see e.g. [41][42][43][44][45][46].
In the current manuscript, we further develop the jet-tomography toolkit, and include the effects of the leading hydrodynamic gradients at all orders in the opacity expansion.
The main result of this paper is the modification of the single parton transverse momentum distribution due to the evolution in a medium of finite longitudinal extension L. In the limit of a longitudinally uniform matter 1 it is given by (16) or (41), and reads where the potential V(x) describes the interaction of the parton with the matter, and . This potential is fixed by the medium model and can be related to an effective dipole cross-section. Here, we model the matter with a background field produced by color sources with number density ρ and screened at distances of order 1 µ . The formula (1) describes the momentum distribution for a final state parton, after it has been produced from an initial hard process (with large energy E), described by dN 0 d 2 pdE , and propagated through a static longitudinally uniform medium with finite transverse gradients ∇ρ and ∇µ 2 of the medium parameters. The primary effect of the gradient terms is to generate a non-trivial angular dependence in the resulting parton distribution.
The present manuscript details the derivation of (1) in two commonly-employed jet quenching formalisms and discusses its properties. In Section II, following the previous work done in [15], we provide a derivation of (1) in the GLV framework, performing a resummation of the associated opacity series. In Section III we show how the same result can be obtained within the amplitude-level resummed framework introduced by Baier, Dokshitzer, Mueller, Peigné, Schiff, and Zakharov (BDMPS-Z) [20,21]. Finally, in Section IV we explore the properties of (1) in a manner suitable for direct future applications 2 in commonly used jet quenching models [47][48][49][50][51]. We summarize our findings and discuss future avenues of research in Section V. Some additional technical details are included into two appendices.

II. THE GLV OPACITY SERIES
In this section, we derive the gradient corrections to the jet broadening order-by-order in opacity expansion, and resum the obtained series. As in [15], we will focus on the spatial gradients of the source density and Debye mass at zero medium velocity, using scalar QCD 3 to describe the underlying theory. Then, the medium-induced color field is static and reads where we use bold font for vectors in the transverse 2D space, v j (q) is a model-dependent medium potential of an individual source numerated with j, while t a j and (x j , z j ) are its color generator and spatial position. Notice that the particular form of the field is derived under an assumption of large source mass. In what follows, we will consider the Gyulassy-Wang (GW) model [18] for the potential where q ⊥ ≡ |q| and µ j is the Debye mass in the HIC context defined by the local medium properties around the ith source.
With this model for the medium, we can turn to the details of the jet-medium interaction.
We start with an initial parton distribution produced by a hard-scattering event. To form a jet the leading parton has to be highly energetic E p z , and we will generally work in the eikonal limit, assuming that any transverse momentum or in-medium characteristic energy scale µ are much smaller than the energy E. Propagating through the matter, jets get modified interacting with the medium-induced color field, and the momentum distribution is also affected. In the perturbative regime, the corresponding change in the momentum distribution can be studied order by order in the coupling g. One should also distinguish the coupling entering the in-medium potential and the emission vertex involving the energetic parton, see e.g. [5].
After the amplitude is obtained, it should be squared and averaged over quantum numbers before one can construct the final momentum distribution. This procedure requires one to further specify the medium model by defining the multipoint correlations of external fields.
Following the prescription commonly used in the pQCD considerations of the jet-medium 3 At eikonal accuracy the spin flips can be ignored, justifying the use of scalar QCD, see e.g. [15]. interactions, we treat the medium color field to be classical and stochastic. We also assume a color neutrality condition requiring that and take into account only pairwise averages, see e.g. [15]. This approximation is motivated by the fact that interference terms are suppressed in a random classical system. We leave the color representation of the sources free, but assume that all of them are in the same representation. Here d tgt is the dimension of the color representation of the sources ("target"), and CR is the quadratic Casimir of the opposite representation. The amplitude squared reads where we have identified the first several orders in the opacity expansion numerated with N , expressing them through the terms in the perturbation series M (p) = r M r with r counting the number of the in-medium field insertions. One should notice that the contributions to |M | 2 involving an odd number of external fields average to zero under our assumptions.
The contribution M r to the full amplitude with r external field entries, see Fig. 1, can be written as where t a proj is the parton ("projectile") color generator, p n = p f − N m=n q m , and p in = p 1 under this numeration. The sums and integrals should be understood as acting on the whole expression, including J (p in ) which depends on the momentum transfers q n .
Assuming for a moment that all the source positions are different, we can perform the q z integrals by residues, noticing that only the small poles of the propagators contribute.
Indeed, the large q z −poles scale as E, and when one of such residues is substituted into the integrand that leads to a sub-eikonal result going beyond our accuracy. Similarly to the original GLV formalism, we also assume that the potentials v j (q) are screened at some scales µ j (varying from point to point), resulting in exponentially suppressed contributions for a sufficiently large and dilute medium µ jn (z jn − z j n−1 ) 1. Then, where J (E, p in ) is the eikonal limit of the source function with p in substituted from the cor- comes from the small pole of the propagator p 2 n = 0, and without loss of generality we set z i 0 = 0, which could be thought of as the center of the source function.
Notice that while the LPM phases are sub-eikonal, they are enhanced by the large medium length, and, thus, should be kept explicitly, c.f. [5,15]. Now one has to square the full amplitude and average over quantum numbers to obtain the momentum distribution modified by the medium. However, looking at (5) one should notice that two fields in the pair can come from the same side of the cut. Then, due to the color neutrality, the general expression (7) involves one or more θ(0) which should be defined. Thus, we have to re-consider the derivation in the case when the amplitude involves two consequent interactions on the same source, or, in other words, we have to separately study the so-called double-Born (DB) diagrams in addition to the ones involving only direct single-Born (SB) interactions. The former are needed to ensure unitarity. If such a contact interaction takes place, involving two consequent insertions n and n + 1 in M r , then the color averaging results in j n = j n+1 , and the corresponding q z integrals read n,z +q n+1,z (zj n −z j n−1 ) where q q z,m is the small pole of p 2 n = 0 solved for q n,z , the combination q (p) n,z + q n+1,z is q n+1,z -independent, and the final expression is written under an assumption that N m=n+2 q z,m is sub-eikonal, as is the case after the q n+2,z integration, for additional details see Appendix A. It should be also mentioned that the full q n+2,z -dependence of the integrand results only in additional screened poles, and the corresponding integration is unaffected by the presence of the contact interaction. Finally, one should notice that we have omitted the q n,z -dependent entry of the form v q n−1 , q (p) n−1,z in the integrand since it cannot modify the q n,z integration, and after the first integration over q n,z its q n+1,z -dependence disappears. Each SB contribution to the amplitude squared should be supplemented with all the corresponding contact terms.
When the amplitude is squared, each contribution to the N th order in opacity involves 2N sums over the in-medium sources. The color averaging reduces the number of the sums to N : only two gluon exchanges are allowed, happening either on different sides of the cut (SB interactions) or on the same source on one of the sides (DB interactions). Commonly, at this step, the discrete sums are replaced by continuous averages with a source number density, i.e.
where the spatial integration goes over the medium volume, and we assume that the medium is large but keep its finite longitudinal length L explicitly.
If the system is uniform in the transverse directions, each x n integral acts only on the corresponding Fourier factor, resulting in where q n and q n are the two momentum exchanges in the averaged potential pair and the sign is different for SB and DB interactions. Thus, the number of the transverse momentum integrals is halved, while all the LPM phases cancel out in the jet momentum distribution.
Collecting the SB and DB terms, one finds the well-known result [27] for the amplitude squared at N th order in opacity: where V(q, z) is a specific combination of the in-medium color potentials, which enters the distribution at all orders in opacity, and z N +1 = L. We will refer to it as the dipole potential, since it can be related to the forward scattering amplitude for a color dipole. Here, 2CR is the full color factor with C proj 1 = t a proj t a proj , and the superscript (N ) indicates the order in the opacity expansion. Notably, the opacity series can be now resummed, since the convolution (10) reduces to a local product in the x-space. Indeed, introducing the jet distribution in the x-space we readily write where for simplicity we set ρ(z) = const, and, consequently, V(q, z) = V(q).
On the other hand, if the medium is inhomogeneous in the transverse directions, then the x n integrals cannot be simplified without further assumptions. As in [15] we focus on the leading corrections in the case of a slow x−dependence, when the thermodynamic parameters can be expanded in their transverse gradients. Then, the leading gradient corrections to the general transverse integral appear due to linear terms, such as where α is a vector index in the transverse 2D space.
In the absence of medium flow, the two medium parameters of interest are ρ and µ, and to the leading order in gradients they can be written as where for compactness we use ρ(z) ≡ ρ(0, z). Integrating the δ-function derivatives by parts, one may find that they act only on the LPM phases, while all other contributions are either suppressed within the eikonal expansion or cancel between complex conjugated contributions. This generalizes the observation in [15] for the broadening at the first order in opacity. We will discuss the details of this derivation in a separate Appendix A using the N = 2 case as an example. Here we only present the squared amplitude at N th order in opacity: where we have introduced notations ρ k , µ k , and V k to distinguish different sources of gra- is an operator generating gradient contributions, and the ordering of the z-integrals is enforced by the θ n,n−1 in (7) and (8). After the variations with respect to ρ k and µ k are performed, we again set the thermodynamic parameters of the same type to be equal and constant in z.
Further simplifying (14) and transforming to the x-space, we find and the opacity series can be again resummed, resulting in where V (x) = ∂ ∂µ 2 V(x). Thus, we have derived (1), which is one of the main results of this work. It gives the Fourier transform of the momentum broadening distribution up to the first order in gradients and to all orders in opacity. One can further use it to study the jet momentum broadening, which we proceed to do in Section IV.

III. THE BDMPS-Z FORMALISM
In this section, we re-derive the leading gradient effects on the jet momentum broadening within the BDMPS-Z approach. From a practical point of view, in this formalism the resummation of multiple field insertions is first performed at the amplitude level by constructing the dressed in-medium propagator. We obtain the in-medium propagator for an inhomogeneous medium, including the leading gradient contributions.
Since the interactions with the medium are dominated by tree-level gluon exchanges, the dynamics of the background field is dominated by the classical Yang-Mills equations, and it can be treated as a stochastic variable. In the BDMPS-Z approach, it is typically assumed that the statistics of the field take a white-noise form due to the large number of uncorrelated degrees of freedom in the medium, analogous to the McLerran-Venugopalan model [32,52,53]. This is equivalent to the assumption used in the previous section that only pairwise averages are non-negligible.
Thus, we again start with a model for the in-medium color field. The model commonly used in the BDMPS-Z approach can be conveniently summarized with where the scattering potential should be set to v a (q) = j e −i(q·x j +qzz j ) t a j v j (q) to coincide with the GW model used in the previous section. It should be also noticed that in general the individual potentials of the scattering centers can be left unspecified in both approaches, although one should be careful treating the DB interactions.
Having the form of the scattering potential, we proceed to rewrite the amplitude given in (7) in terms of an effective dressed propagator. For that, we first write (6) using the in-medium color field (17), then where the temporal components satisfy the constraint p 0 n = E, indicating that there is no energy transfer between the medium and the probe via soft gluon exchanges. It is convenient to Fourier transform the potentials, simplifying the p z -integration. This is equivalent to working in a mixed representation, where "time" dependence is made explicit, commonly employed in the BDMPS-Z related literature and analogous to old-fashioned perturbation theory. Then, (18) can be rewritten as where (x n , z n ) denote the interaction points, and should not be mixed with the source coordinates in the GW model, now hidden in v a .
Finally, noticing that p 1 = p in and p r+1 = p f , we can write the perturbative amplitude as a convolution between the source and a contribution to an effective single particle propagator Notice that the in-medium potential v a is screened by the Debye mass µ, and thus has a finite spatial support of size ∼ 1 µ . As a consequence, the z-integrals can be safely taken to run from the production point z = 0 to the end of the medium at z = L, which enters the single particle propagator.
The full effective propagator can be obtained by summing over the number of interactions where it is simple to check that in the case of vacuum propagation it reduces to the usual Feynman result Inserting this back in (20), we find that the amplitude reduces to the initial source function, as expected.
Instead of dealing with the series in (21), one can construct an evolution equation for G from (19) and (20), which takes the usual Schrödinger-like form This equation should be supplemented with an "initial condition", which can be obtained from the fact that at L = 0 there is no modification to the amplitude sourced by J(p in ). The solution to (23) with the corresponding initial condition is well known [54], and its x-space form can be written as a path integral where P indicates path ordering, and x L ≡ x(L) and x 0 ≡ x(0) are the boundary conditions for the trajectory. This effective propagator can be thought of as describing a massive nonrelativistic particle, moving from the initial position x 0 at "time" τ = 0 to the final position With the further assumption that the QCD emission vertices are unaltered in the medium, one can derive a set of effective Feynman rules using the propagator above and compute any quantum amplitude. As a consequence, in such a path integral formulation of the BDMPS-Z formalism, in practice one can just draw all the relevant time ordered Feynman diagrams including the medium and directly obtain the amplitudes, similar to more standard vacuum pQCD calculations. The squared amplitude, already averaged over the quantum numbers and medium configuration, can be easily expressed through an in-medium correlation function of two propagators, c.f. (20), as well as the distribution corresponding to the jet momentum broadening itself.
In order to compute (25), we first have to revisit how the averaging procedure is performed where we have used the color neutrality condition and taken the continuous limit of the distribution of scattering centers in the medium (as in the previous section) in order to make the dependence on the source density explicit. Before we proceed, one should notice that only a particular limit of this correlation function enters the amplitude (19) when it is squared and averaged, see e.g. [55]. In fact, the average in (26) is Fourier transformed to the momentum space in the amplitude squared, with the z-component of the momentum being sub-eikonal. Thus, the q z dependence in (26) can be neglected without affecting the final result at the accuracy level being considered. Using this simplification, we find where x in ρ(x, z) and µ 2 (x, z) has been replaced with −i ∂ ∂(q−q) acting on everything but the delta function δ (2) (q − q), and we assume that ρ and µ 2 are constant in the longitudinal direction. One should notice that all the previous steps in this section are unaffected by the inhomogeneity of the medium, and the gradient effects enter solely through the potential averages.
The two-point correlator of the in-medium color potentials now depends not only on the transverse size of the effective color dipole formed by the propagating parton in amplitude and conjugate amplitude (see Fig. 2), which is proportional to the difference |r − r|, but also on its center of mass transverse position r+r 2 . This is well expected, since translation invariance is now violated by the transverse gradients. The higher-order gradients will enter the BDMPS-Z construction similarly, with higher powers of the transverse position of the center of mass.
Given the leading gradient form for the pairwise average of two in-medium potentials, we now turn to the average entering (25). Using the two-point correlator (27), we can write the average of two Wilson lines as 4 where the contact terms arise from the pairwise averages of potentials coming from the same exponential, and we use the fact that the potential is Hermitian in the GW model. Notice that the exponential with the gradients in the argument should be treated as a series valid up to the first order at the accuracy of our consideration.
The position space form of the relevant correlator of two propagators reads where u ≡ r − r and w ≡ r+r 2 . The path integral is simple enough to be evaluated analytically, since the w-dependence enters only through the small gradient terms. At the first order in gradients the final answer is proportional to the value of the integrand on the solution of the effective classical equation of motion 5 (EOM). Thus, using the standard methods [54][55][56], we can reduce this expression to where u c (τ ) is the classical solution satisfying with appropriate boundary conditions. Notice that the structure in the denominator of (30) is a part of a Jacobian appearing in the path integral, along with the overall multiple. In c is linear in gradients 6 . At the zeroth order, the right-hand side of the EOM is zero, and one readily finds that the separation vector can only change linearly with "time" leading to the uniform broadening result in (12), see e.g. [55,56]. In turn, the leading gradient correction to the trajectory is purely imaginary and satisfies the trivial boundary conditions u The momentum scale p f in (25) corresponds to a measured quantity and has to be matched in the two propagators. For this particular projection of the two-point correlator as a perturbation around the trivial trajectory, see Appendix B. 6 Notice that the imaginary part of the classical EOM is proportional to the gradient terms. Thus, substituting u c (τ ) into the integrand is equivalent to acting on it with a shift operator changing the zero acceleration of the leading order trajectory with a complex function.
(30) we find where again the delta functions of complex arguments should be understood as a shift operator (with small imaginary shift parameter) acting on the corresponding delta function of the (leading) real part of the argument. Contracting this with the initial source functions in (25), one can find the corresponding squared amplitude where, as in the derivation of (10), we have assumed that J has at most a constant imaginary phase, as it is for a tree-level 2-to-2 process, and introduced the symmetric momentum We have also used that the leading shift of zero at the first order in gradients, since the first derivative of this symmetric function with respect to p in − p in is zero at p in = p in .
Thus, we can write the corresponding Fourier transformed distribution as The velocity constraint at τ = L allows to simplify the form of u c . In particular, in the case of vanishing matter gradients, it implies directly from (32) that the size of the effective dipole is preserved during its evolution [4]. In the present case, this constraint implies and one can use it to relate u 0 and u L .
At the leading order in gradients we find and, consequently, u c takes a remarkably simple form The leading order term is the well known constant solution. The net effect of including gradients is to change the effective trajectory of the dipole separation with a constant (imaginary) shift in the acceleration which is fixed by the form of V(x). One should also notice that the delta function constraining the velocity at τ = L has a non-trivial dependence on u 0 which enters the argument also through V u (0) c (ξ) , leading to an additional multiple: Combining all the previous results and expanding in the smallness of gradients, we conclude that the x-space form of the jet distribution reads and one can easily see that it exactly coincides with (16), obtained from the opacity series resummation. It should be mentioned that the two measure factors coming from the path integral and velocity constraint cancel, resulting in the simpler expression above. Thus, we have shown that the usual BDMPS-Z formalism can be extended to inhomogeneous backgrounds to the leading order in gradient expansion. Conceptually, the approach can be extended in a straightforward manner to include higher order gradient effects. However, we note that in such a situation the path integrals to be solved are technically more involved (e.g. the w dependence in (29) would no longer be linear).
Before we turn to the properties of the distribution (41), it should be also mentioned that its extension to the case of z-dependent medium profile can be readily obtained. After a straightforward generalization of the derivation in the longitudinally homogeneous case, one finds whereĝ (τ ) = ∇ρ(τ ) δ δρ + ∇µ 2 (τ ) δ δµ 2 .

IV. THE FINAL STATE DISTRIBUTION AND ITS PROPERTIES
In this section, we proceed to discuss the properties of the final jet momentum distribution given by (16) (or equivalently (41)). Heuristically, the leading effect of the matter gradients is that the broadening becomes anisotropic, and the final jet momentum distribution is direction dependent, see Fig. 3. In other words, propagating through an inhomogeneous matter, probes pick up an additional transverse momentum, and even on average it is nonzero due to the non-trivial matter structure. The modifications to the jet structure are thus qualitatively different from the ones observed in the case of a homogeneous background. In this section, we will analyze the properties of the distribution (16) and illustrate our results with simple numerical estimates for the broadening probability, using phenomenologically relevant parameters. Let us first notice that the broadening of the jet distribution has to be unitary -at a fixed energy, the number of jets (partons) cannot be changed and the in-medium propagation results only in a reshuffling of the underlying momentum modes.
To check that, one can consider the value of the x-space distribution at x = 0, which corresponds to the volume integral of the momentum-space distribution. At this point, the This modification depends on the angle θ between ∇T and the transverse momentum accumulated during the evolution, as shown by the dashed shape.
x-space dipole potential is zero V(0) = 0 as well as its transverse gradient ∇V(0) = 0, and the full distribution (1) is equal to d 2 p dN (0) d 2 p dE . Thus, the unitarity condition is satisfied by the broadened jet momentum distribution even in an inhomogeneous medium One could also notice that in (41) the initial distribution is not fully factorized from probe-medium interactions. Indeed, in the uniform case, the x-space form of the final jet distribution can be represented as a product of the initial distribution with an interaction factor. The gradient effects in turn enter not only through the probe-medium interactions, but also through a shift of the argument in the initial hard distribution dN (0) d 2 p dE . Thus, the relation between the final and initial distributions involves an operator where is a transform of P(p), the (broadening) probability for a particle of energy E to acquire transverse momentum p due to propagation in the medium for a distance L, and is a shift operator being an identity operator in the absence of gradients. One can immediately see that P(p) is self-normalized, and thatŜ (x) does not change particle number at a given fixed energy.
It is instructive to compare the resummed final state distribution to the one obtained at the first order in opacity by evaluating the leading moments: where F (p) is an arbitrary function. Following [15], we focus on a Gaussian initial distribu- where f (E) is an unspecified energy dependence, and w is the characteristic width. We also assume that the matter is uniform in the longitudinal direction. One generally expects that the terms linear in gradients cannot modify the scalar moments, so p 2k ⊥ = p 2k ⊥ ∇ρ=∇µ 2 =0 , and we may focus on the odd moments p p 2k ⊥ sensitive to the anisotropy 7 , see e.g. [15,17]. Starting with the simplest case of the averaged momentum p at the given energy, we find it convenient to use the opacity expanded form of the distribution (14) rather than its resummed form. Doing so, one finds that the moment is controlled by the following integral where p n = p f − N m=n q m , p in = p 1 , and α and β are the indices in the transverse 2D space.
We also notice that in the absence of the gradients the averaged momentum is zero due to the lack of a preferred direction in the problem, and the corresponding part of (14) gives zero. Thus, for the contribution at the N th order in opacity, we find where we have additionally used that d 2 q n V n (q n ) = 0 removing the w 2 term. For any N > 1 this expression is zero due to the same property of V n (q n ), since the integrand involves 7 For an earlier attempt to study the matter anisotropy with directional jet quenching, see [57].
only two powers of q i , and at least one of the dipole potential averages to zero. In turn, the N = 1 contribution is zero since it involves only the w 2 term. Thus, we find that p = 0 to all orders in opacity.
The first non-zero moment of the momentum broadening distribution at the first order in opacity is p p 2 ⊥ [15]. We can again perform the final momentum averaging, and find Thus, one can see that the only non-zero contributions may come at N = 1 and at N = 2 while the higher orders in opacity decouple. In the case of the GW potential, the terms proportional to ∇ρ are divergent, and should be regularized at some scale. Comparing with [15], we set |q max | = √ Eµ, and find The first term in the expression above precisely agrees with the result obtained in [15]. In turn, the second term is new and indicates that the higher N contributions to the transverse momentum moments with integer k ≥ 1 are generally dominating as long as the potential integrals are divergent. One should also notice that the second term is independent of w, and can be identified as a purely final state effect (described by the broadening probability P(p)), since in the limit w → 0 the initial distribution is constant in coordinate space and not affected by the shift operatorŜ. In general, non-zero odd moments are the primary effect of the hydrodynamic gradients on the jet broadening, and they can be used to learn about the medium profile with jets. We leave the moments with general real k and their phenomenological implications for future studies.
We now turn back to the partial factorization, and focus on the final state effects. The corresponding portion of (44) is the self-normalized broadening probability P(p), which has multiple phenomenological applications [47,[58][59][60][61]. Its Fourier transform can be obtained from (41) if we set E dN (0) d 2 xdE to be a constant in the transverse directions, decoupling in this way the shift operatorŜ(x): At this point one has to specify the particular model for the medium, fixing the dipole potential, and evaluate P(p) explicitly. We consider the GW model, and in this case the potential reads However, even in this simple case and in the absence of gradients, the Fourier transform cannot be written in a closed form, restricting phenomenological applications.
By this reason, it is instructive to consider P(p) in some limiting regimes. If the gradients are zero, then the probability (53) is a function of two dimensionless combinations of the parameters and momentum: p ⊥ µ and χ ≡ Cg 4 ρ 4πµ 2 L, where χ is the opacity in the GW model. The dipole potential V(x ⊥ ) tends to a constant at large µ x ⊥ , and the large x ⊥ limit of the integration is controlled by the fast oscillation of the Fourier exponential. For sufficiently large momenta such that p ⊥ µ, one can expand the dipole potential in powers of µ x ⊥ , and for the GW model it reads where γ E is the Euler constant. One should notice that the first term in the expansion is common between most of the medium models, since it is fixed by the ultraviolet (UV) Coulomb-like behavior for large momentum exchanges, see e.g. [62][63][64].
If we further require that p 2 ⊥ χ µ 2 , then the exponential of the dipole potential itself can be expanded, and, after expanding the in-medium potentials, we find where we have omitted terms proportional to δ (2) (p) since the argument is away from zero.
The first two terms in this expansion come from the homogeneous case, with the first giving the dominant Coulomb tail expected at large momentum transfers. In turn, the gradient term has an odd (negative) power of p ⊥ and can only be generated in an inhomogeneous medium. In fact, one can compare the asymptotic structure of (56) directly with the first non-trivial moment (52) in the limit of infinitely narrow source (i.e. w → 0). The only contribution to (52) in this limit comes from the N = 2 term, which gets a double Coulomb logarithmic enhancement. Indeed, looking at the last term in (56), one may notice that its dominant contribution to p p 2 ⊥ scales as ∇ρ Eρ dp ⊥ log p ⊥ p ⊥ , and with a similar UV regularization it gives ∇ρ Eρ log E µ 2 in qualitative agreement with the previous result.
Expanding the probability distribution P in two parameters, one has to pay a particular attention to the applicability of the result. Indeed, if the large momentum expansion parameter µ p ⊥ is too small, the gradient contributions could take the leading role. Here we assume that this is not the case, and that the first order gradient correction appearing in the last term of (56) is at least smaller than the dominant Coulomb tail contribution, although it can compete with the second term in the expression.
Another kinematic regime commonly considered in the literature [47,63] corresponds to the intermediate momentum region µ 2 p 2 ⊥ ≤ χ µ 2 . In this case, the exponential of the dipole potential cannot be expanded, while the approximation (55) can still be utilized.
Moreover, the weak logarithmic dependence of the expanded dipole potential on x can be neglected. Indeed, one can introduce a new scale Q 2 µ 2 such that In practice, it can be fixed up to an overall coefficient based on phenomenological arguments, see e.g. [63,65] for further discussion. In this regime, the dipole potential is quadratic, and the leading broadening probability (53) becomes Gaussian. This is a regime of multiple soft interactions, which is widely used in phenomenological models for jet quenching, see e.g. [47,48,58], and often referred to as BDMPS-Z/ASW approximation [22,66]. Taking into account the gradient effects, we find The usual Gaussian result gets an overall p-dependent modulation by a new term containing all the gradient effects. We note that the effect of this contribution should be more important when p 2 ⊥ ∼ χµ 2 log Q 2 µ 2 , roughly at the peak of the p ⊥ -distribution [63]. In this region, the leading behavior of P is fully described neither by a single hard scattering nor by the multiple soft scatterings , but rather there is a competition between these two regimes. Thus, one may expect a more pronounced deviation from the homogeneous solution in this region. We numerically verify this observation below in Fig. 4.
Finally, for µ 2 p 2 ⊥ ≤ χ µ 2 , we can also consider the small logarithmic correction to the dipole potential, neglected in (57), as a perturbation This approach to capture the leading effects beyond the quadratic Gaussian approximation is known as the improved opacity expansion (IOE), see e.g. [62,67,68]. Substituting the potential into the broadening probability, we write it as an expansion P IOE (p) = P 0+1 (p) + P 1+0 (p) + P 2+0 (p) + P 1+1 (p) + P 2+1 (p) + P 3+1 (p) + ... , where the first number in the subscript corresponds to the order in δV, and the second one counts the order in gradients. The leading contribution P 0+1 (p) is given by (58), the term P 1+0 (p) has been discussed in details in [63], and we have to evaluate the mixed contributions. The full derivation of the sub-leading terms goes beyond the scope of this paper, and we focus on their large momentum regime. That will allow us to see that the terms sub-leading in δV give the correct large momentum limit (56) even in the presence of gradients, c.f. [63].
All the x-integrals entering (60) up to the second order can be obtained from three master integrals if we act on them with the appropriate number of momentum derivatives. They read where we omit the explicit form of I 3 for brevity. In the large p ⊥ limit, the momentum derivatives acting on a term reduce its contribution removing powers of p ⊥ . Consequently, one can compare the relative importance of the different contributions in (60) by counting the powers of x entering through V and δV outside the Gaussian exponential. However, since the structure of the gradient terms in (53) involves gradients of the dipole potential, we have to consider higher order terms in δV to be sure that we keep all the relevant contributions.
We first notice that the leading term in the IOE expansion P 0+1 (p) is exponentially gives the leading Coulomb term in this limit. Similarly, the sub-leading non-gradient contribution satisfies P 2+0 (p) ∝ 1 The mixed gradient terms in turn scale as P 1+1 (p) ∝ 1 can check that the next contribution P 3+1 (p) ∝ 1 p 6 ⊥ log 2 p ⊥ is suppressed both by a higher power of inverse momentum as well as by the gradients, and can be omitted. Finally, we which precisely agrees with the large momentum limit (56), showing that the IOE covers both limits considered above. It could be also noticed that the logarithmic terms in P IOE (p) are sensitive to a cancellation of Q between terms of different orders in the IOE, c.f. [65].
Finally, to have a more quantitative understanding of the broadening distribution P(p), we evaluate (53) using the full GW potential given in (54). We also assume that the nuclear matter is near to equilibrium, and its properties are controlled by a single parameter. For instance, in the case of QGP in the large temperature limit, one can use parametric scaling motivated by the equilibrium thermodynamics µ ∝ gT and ρ ∝ T 3 , then The probability P(p) is now a function of θ, the angle between p and ∇T , and a dimensionless quantity c T ≡ ∇T ET 1, see Fig. 3. Performing the remaining angular integration, we can write P as For the particular case of the GW model, V GW (x ⊥ ) asymptotically tends to a constant.
Therefore, the integral in (64) needs to be regulated for values of x ⊥ 1 µ . In particular, one should remove the non-scattering probability term associated with the asymptotic behavior of V GW (x ⊥ ), so that only genuine broadening contributions are taken into account. Thus, following e.g. [63,69], we consider a regularized broadening distribution which differs from (64) by a singular term (delta function) at p = 0.
In Fig. 4 we present numerical evaluations of (65) for several sets of parameters inspired by the particular values used for phenomenological estimates in the literature, see e.g. [17,49,63,[70][71][72][73]. However, we have to additionally stress here that (65) has been evaluated in the simplest limit of a longitudinally uniform slab of nuclear matter with the medium parameters varying in the transverse directions. Thus, any realistic phenomenological applications would require further generalizations of this broadening distribution. Each plot is characterized by a value for the medium opacity χ, the medium length size L, and the Debye mass µ. The values considered for the medium parameters are inspired by the ones typically found in the literature [17,49,63,[70][71][72].
integration up to a scale of the order of 1 µ . This ultraviolet cutoff was systematically varied in order to ensure that the final results do not depend on its particular choice.
At the qualitative level, we observe that the gradient effects are small for lower opacity, since the gradient corrections in (65) are proportional to χ 2 . Comparing all the results shown, we also confirm the claim made above that gradient effects should become more important around the peak of the distribution. Indeed, for the higher opacity scenarios, gradient effects lead to a sizable suppression of the broadening peak (for c T = 0.02 and θ = 0), resulting in an increase of the broadening probability at higher p ⊥ . Conversely, in the case where p is antiparallel to ∇T , we observe an increase of the distribution peak, leading to a depletion of higher momentum modes. Another feature shared by all the setups is the insensitivity of the result to gradient effects at small p ⊥ , due to the fact that the inhomogeneous term always couples to an odd power of p. Finally, it should be also mentioned that the large momentum tail of the distribution is not modified by the gradient corrections, as seen from (62). This is indeed true for the numerically evaluated distributions, but the corresponding large p ⊥ region is not presented in Fig. 4.

V. CONCLUSIONS AND OUTLOOK
In this work, we have derived the parton momentum broadening distribution in an inhomogeneous dense nuclear matter. The final distribution has been obtained using both the GLV and BDMPS-Z formalisms. In the homogeneous case, these two approaches considered above were proven to agree order by order in opacity, see e.g. [27,[74][75][76]. Here we make the next step, showing that the agreement holds even if the matter is not translationallyinvariant in the transverse directions. Importantly, the effects of the medium inhomogeneity only emerge once the medium averaging is performed. As a consequence, all the results for the jet-medium interactions which are independent of the averaging procedure stay unmodified with respect to the homogeneous case. Qualitatively, the modification of the resulting broadening distribution (comparing to the homogeneous baseline) becomes relevant around and above its characteristic peak, introducing a novel angular dependence. As a consequence, the gradient effects can have sizable impact in phenomenological applications, even for momentum averaged quantities.
Our results for the jet broadening can be extended to the case of inelastic energy loss in a dense inhomogeneous medium. Primarily, such an exercise will allow to study how the gradient effects alter the in-medium gluon production rate, giving rise to a non-trivial angular structure of the radiation. One may expect that the effective time factorization of the emission spectrum found in the limit of soft induced radiation in homogeneous media [58,77] would be modified in accordance with (44). Understanding the resulting new factorized form is critical for some aspects of jet quenching phenomenology [47,48,58].
Another important goal we leave for the future studies is to design a jet observable sensitive to the matter effects discussed in this manuscript, c.f. [17]. Particularly, it would be interesting to explore to what extent the jet substructure grooming/tagging techniques, see e.g. [78][79][80][81][82][83] could be used to look for matter inhomogeneity effects on the angular spectrum of the hardest substructures within the jet. Such new observables would provide a window to further probe the medium properties locally along the ideas of jet tomography.
We also note that the presented results are independent of the physical process which generates the background medium and, thus, they are in principle applicable to either the QGP phase or the glasma phase in the HIC context [84,85]. The latter one could be considerably anisotropic, resulting in stronger directional effects. It would be interesting to better understand the relation between the jet quenching formalisms considered here and other approaches which can describe jet evolution in inhomogeneous and evolving media, see e.g. [86][87][88]. In turn, extension to the DIS context would require the present calculation to be updated to the corresponding kinematical regime with particular focus on the factorization of the jet production and its propagation through the medium, see e.g. [29,[89][90][91][92] and references therein.
Finally, it should be mentioned that the dependence of the jet-medium interaction on hydrodynamic gradients could be used to probe the strength of the interactions in the underlying theory. Indeed, comparing the characteristic transport properties and the general jet behavior in an inhomogeneous matter between pQCD and holographic models for jetmedium interactions, see e.g. [48,[93][94][95][96][97][98], one may hope to identify a set of new observables distinguishing the two regimes. We leave this intriguing opportunity for future work. this contribution and averaging over the quantum numbers, we find where the x-dependence has been transformed to momentum space and we assumed that the matter is longitudinally uniform. Notice that the momentum derivatives of the form ∂ ∂(q−q)α , where q is the corresponding momentum in the conjugated amplitude, can act only on the LPM phases at the first order in gradients, since the rest of the integrand is symmetric under exchange of the two momenta in each pair.
Now we turn to the two mixed contributions in the amplitude squared, which involve both SB and DB interactions. It is natural to consider them partially averaged (we choose the contact interaction pair and perform color averaging), then where the subscript indicates that the first and second insertions are SB and DB correspondingly, and At N = 2 this diagram contributes to the amplitude squared multiplied by the diagram with a single SB insertion. Performing the residual averaging, one finds where we have used that the momentum derivatives cancel between the complex conjugated contributions, unless they act on the LMP phases.
Similarly, we can consider the opposite order of the SB and DB insertions, and the corresponding amplitude reads The integral appearing due to the contact interaction is given by exhibiting a non-zero real part even at the leading order in the eikonal expansion. Herẽ µ 2 n ≡ µ 2 +q 2 n , and we keep the subscript i 1 implicit inμ. One can notice that I 21 i 1 results only in q 3,z -poles, leading to exponentially suppressed contributions, and we can easily perform the last q z -integration. Then, (2π) 6 t a proj t a i 3 θ 2,1 θ 1,0 e −i(q 1 +q 2 )·x i 1 e −iq 3 ·x i 3 I 12 where we have used that Re I 21 is sub-eikonal after the q 3,z is substituted (the two integrals I 12 and I 21 coincide in this limit). Turning to the contribution to the squared amplitude, we find where we have adjusted the integration variables to match the form of the other contributions.
Finally, we have to consider two types of diagrams involving two DB interactions. The first of such amplitude reads iM 22 = 1 4 C 2 i 1 ,i 3 d 2 q 1 d 2 q 2 d 2 q 3 d 2 q 4 (2π) 8 θ 2,1 θ 1,0 e −i(q 1 +q 2 )·x i 1 e −i(q 3 +q 4 )·x i 3 × e −i (p−q 1 −q 2 −q 3 −q 4 ) 2 where we have used the properties of the contact q z -integrals from M 12 and M 21 omitting unnecessary details. Combining with the conjugated amplitude, we write where the numeration of momenta has been changed for convenience.
Similarly, one has to take into account the DB amplitude squared. Using the result of [15] we write and notice that the z-integration limits differ from the limits in other diagrams. In fact, in the case of a longitudinally invariant matter, it can be easily seen that this pair of integrals is twice larger than the ordered pair of z-integrals.
Thus, collecting all the contributions, we find the full amplitude squared at N = 2, and it reads in agreement with the general expression (14). Repeating the same calculation for the higher orders in opacity, we find that the structure of (14) is reproduced. In this Appendix, we briefly discuss the imaginary term in the EOMs (31), arising due to the gradient effects. This term results in an imaginary shift of the zeroth order trajectory, which enters for instance into the argument of δ-functions in the last line of (34). Since the δ-functions are originally defined in the real space, this notation is only formal, and the functions should be re-defined.
In the main text, the resulting delta functions are treated as in the case of a real argument.
To justify this treatment, consider the following one-dimensional integral dx e −ipx e ax ≡ (2π) δ(p + ia) .
For a real a, one can understand it as a momentum representation of the shift operator e ax = e ia ∂ ∂p with a purely imaginary argument ia. Thus, we can write dx e −ipx e ax = (2π)e ia ∂ ∂p δ(p) = (2π) δ(p + ia) .
This simple calculation shows that the imaginary term entering the EOM will indeed give rise to a momentum-space shift operator, and its effect can be formally written with the complex argument of a δ-delta. For instance, the one-dimensional analogue to the integral of interest in the main text reduces to dp in dp in δ(p in − p in − ia)J(p in )J † (p in ) = dp in |J(p in )| 2 + ia dp in dp in (2π) δ(p in − p in ) ∂ ∂(p in − p in ) J(p in )J † (p in ) = dp in |J(p in )| 2 + ia 2 dqdQ (2π) δ(q) ∂ ∂q J Q − q 2 J † Q + q 2 = dp in |J(p in )| 2 . (B3)