# Jet broadening in dense inhomogeneous matter

João Barata<EMAIL_ADDRESS>Instituto Galego de Física de Altas
Enerxías, Universidade de Santiago de Compostela, Santiago de Compostela
15782, Spain  Andrey V. Sadofyev<EMAIL_ADDRESS>Instituto Galego de
Física de Altas Enerxías, Universidade de Santiago de Compostela, Santiago de
Compostela 15782, Spain  Institute for Theoretical and Experimental Physics,
Moscow 117218, Russia  Carlos A. Salgado<EMAIL_ADDRESS>Instituto
Galego de Física de Altas Enerxías, Universidade de Santiago de Compostela,
Santiago de Compostela 15782, Spain

###### Abstract

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.

###### Contents

  1. I Introduction
  2. II The GLV Opacity Series
  3. III The BDMPS-Z Formalism
  4. IV The Final State Distribution and Its Properties
  5. V Conclusions and Outlook
  6. A The Gradient Corrections at the $N$th Order in Opacity
  7. B Shift Operator and $\delta$-function

##  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 Bjorken1982EnergyLO , for a more recent review see
Mehtar-Tani:2013pia ; Qin:2015srf ; Blaizot:2015lma ; Sievert:2018imd . 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. Vitev:2002pf ; JET:2013cls ; Betz:2014cza ;
Apolinario:2017sob ; Li:2020zbk ; Li:2020rqj ; Arratia:2020nxw ; He:2020iow ;
Apolinario:2020uvt ; Sadofyev:2021ohn ; Du:2021pqa ; Antiporda:2021hpk  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 Gyulassy:1993hr ;
Baier:1996sk ; Zakharov:1996fv ; Baier:1996kr ; Zakharov:1997uu ; Baier:1998kq
; Gyulassy:1999zd ; Gyulassy:2000fs ; Gyulassy:2000er ; Gyulassy:2002yv ;
Wang:2001ifa ; Zhang:2003wk ; Arnold:2002ja ; Djordjevic:2003zk ; Mehtar-
Tani:2006vpj . 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
Sadofyev:2021ohn . 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 Gyulassy:2000gk ; Gyulassy:2001kr ; Baier:1998yf
, where the medium dilution was considered, and to Armesto:2004pt ;
Armesto:2004vz , 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. Baier:2006pt ; Liu:2006he ; Renk:2006sx . Only more recently, the medium
evolution effects on the jet-medium interactions were formally included into
the Gyulassy-Levai-Vitev (GLV) opacity expansion framework Sadofyev:2021ohn .
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, Sadofyev:2021ohn  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. Lekaveckas:2013lha ; Rajagopal:2015roa ; Sadofyev:2015hxa ;
Li:2016bbh ; Reiten:2019fta ; Arefeva:2020jvo .

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 matter111An extension of (I) to an arbitrarily $z$-dependence of the
medium is given in (III). it is given by (II) or (III), and reads

$\displaystyle(2\pi)^{2}\frac{d\mathcal{N}}{d^{2}\bm{p}dE}$
$\displaystyle=\int
d^{2}{\bm{x}}\,e^{-i{\bm{p}}\cdot{\bm{x}}}e^{-\mathcal{V}(\bm{x})L}\Bigg{\\{}1-i\frac{\mathcal{V}(\bm{x})L^{3}}{6E}\Bigg{(}\frac{\mathcal{V}^{\prime}(\bm{x})}{\mathcal{V}(\bm{x})}\bm{\nabla}\mu^{2}+\frac{1}{\rho}\bm{\nabla}\rho\Bigg{)}\cdot\bm{\nabla}\mathcal{V}(\bm{x})$
$\displaystyle+i\frac{\mathcal{V}(\bm{x})L^{2}}{2E}\left(\frac{\mathcal{V}^{\prime}(\bm{x})}{\mathcal{V}(\bm{x})}\bm{\nabla}\mu^{2}+\frac{1}{\rho}\bm{\nabla}\rho\right)\cdot\bm{\nabla}\Bigg{\\}}\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{x}dE}\,,$
(1)

where the potential $\mathcal{V}(\bm{x})$ describes the interaction of the
parton with the matter, and
$\mathcal{V}^{\prime}(\bm{x})\equiv\frac{\partial}{\partial\mu^{2}}\mathcal{V}(\bm{x})$.
This potential is fixed by the medium model and can be related to an effective
dipole cross-section. The formula (I) 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 $\frac{d{\cal N}^{0}}{d^{2}\bm{p}dE}$,
and propagated through a static longitudinally uniform medium with finite
transverse gradients $\bm{\nabla}\rho$ and $\bm{\nabla}\mu^{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 (I) in two commonly-employed
jet quenching formalisms and discusses its properties. In Section II,
following the previous work done in Sadofyev:2021ohn , we provide a derivation
of (I) 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) Baier:1996kr ; Zakharov:1996fv
. Finally, in Section IV we explore the properties of (I) in a manner suitable
for direct future applications222For a reader interested in immediate
phenomenological applications and not in the technical derivation, (I) (or its
generalization (III)) and Section IV contain all the necessary information. in
commonly used jet quenching models Caucal:2019uvr ; Casalderrey-Solana:2014bpa
; Sievert:2019cwq ; Putschke:2019yrg ; He:2015pra . 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
Sadofyev:2021ohn , we will focus on the spatial gradients of the source
density and Debye mass at zero medium velocity, using scalar QCD333At eikonal
accuracy the spin flips can be ignored, justifying the use of scalar QCD, see
e.g. Sadofyev:2021ohn . to describe the underlying theory. Then, the medium-
induced color field is static and reads

$\displaystyle gA_{\mathrm{ext}}^{\lambda a}(q)=-(2\pi)\,g^{\lambda
0}\,\sum_{j}e^{-i\left(\bm{q}\cdot\bm{x}_{j}+q_{z}z_{j}\right)}\,t_{j}^{a}\,v_{j}(q)\>\delta\left(q^{0}\right)\,,$
(2)

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_{j}^{a}$ and $(\bm{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 Gyulassy:1993hr  for the potential

$\displaystyle
v_{j}(q)\equiv\frac{-g^{2}}{-q_{0}^{2}+q_{\perp}^{2}+q_{z}^{2}+\mu_{j}^{2}-i\epsilon}\,,$
(3)

where $q_{\perp}\equiv|\bm{q}|$ and $\mu_{j}$ is the Debye mass in the HIC
context defined by the local medium properties around the $i$th 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

$\displaystyle
E\,\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{p}dp_{z}}\equiv\frac{1}{2(2\pi)^{3}}\left|J(p)\right|^{2}\,,$
(4)

produced by a hard-scattering event. To form a jet the leading parton has to
be highly energetic $E\simeq p_{z}$, and we will generally work in the eikonal
limit, assuming that any transverse momentum or in-medium characteristic
energy scale $\mu$ 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. Sievert:2018imd .

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 interactions, we
treat the medium color field to be classical and stochastic. We also assume a
color neutrality condition requiring that

$\left\langle
t^{a}_{i}t^{b}_{j}\right\rangle=\frac{1}{d_{\mathrm{tgt}}}\textrm{tr}\left(t^{a}_{i}t^{b}_{j}\right)=\frac{1}{2C_{\bar{R}}}\delta_{ij}\delta^{ab}\,,$

and take into account only pairwise averages, see e.g. Sadofyev:2021ohn . 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 $C_{\bar{R}}$ is the quadratic Casimir of the opposite
representation. The amplitude squared reads

$\displaystyle\left\langle|M|^{2}\right\rangle$
$\displaystyle=\underbrace{\left\langle|M_{0}|^{2}\right\rangle}_{N=0}+\underbrace{\left\langle|M_{1}|^{2}\right\rangle+\left\langle
M_{2}M_{0}^{*}\right\rangle+\left\langle M_{0}M_{2}^{*}\right\rangle}_{N=1}$
$\displaystyle+\underbrace{\left\langle|M_{2}|^{2}\right\rangle+\left\langle
M_{3}M_{1}^{*}\right\rangle+\left\langle
M_{1}M_{3}^{*}\right\rangle+\left\langle
M_{4}M_{0}^{*}\right\rangle+\left\langle
M_{0}M_{4}^{*}\right\rangle}_{N=2}+...\>,$ (5)

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)=\sum\limits_{r}M_{r}$ with $r$ counting the number of the in-
medium field insertions. One should notice that the contributions to
$\left\langle|M|^{2}\right\rangle$ involving an odd number of external fields
average to zero under our assumptions.

Figure 1: The contribution $M_{r}$ with four field insertions ($r=4$) to the
full amplitude.

The contribution $M_{r}$ to the full amplitude with $r$ external field
entries, see Fig. 1, can be written as

$\displaystyle iM_{r}(p_{f})$
$\displaystyle=\prod_{n=1}^{r}\Bigg{[}(-1)\sum_{j_{n}}\int\frac{d^{2}\bm{q}_{n}\,dq_{n,z}}{(2\pi)^{3}}\,t^{a}_{\mathrm{proj}}t_{j_{n}}^{a}$
$\displaystyle\hskip 85.35826pt\times
e^{-i\left(\bm{q}_{n}\cdot\bm{x}_{j_{n}}+q_{n,z}z_{j_{n}}\right)}\,\frac{2E}{p_{n}^{2}+i\epsilon}\,v_{j_{n}}(q_{n})\Bigg{]}J\left(p_{in}\right)\,,$
(6)

where $t^{a}_{proj}$ is the parton ("projectile") color generator,
$p_{n}=p_{f}-\sum\limits_{m=n}^{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\left(p_{in}\right)$ 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 $\mu_{j}$ (varying from point to point), resulting in exponentially
suppressed contributions for a sufficiently large and dilute medium
$\mu_{j_{n}}(z_{j_{n}}-z_{j_{n-1}})\gg 1$. Then,

$\displaystyle iM_{r}$
$\displaystyle=\prod_{n=1}^{r}\Bigg{[}\sum_{j_{n}}\int\frac{d^{2}\bm{q}_{n}}{(2\pi)^{2}}\,it^{a}_{\mathrm{proj}}t_{j_{n}}^{a}\,\theta_{j_{n},j_{n-1}}$
$\displaystyle\hskip 85.35826pt\times
e^{-i\bm{q}_{n}\cdot\bm{x}_{j_{n}}}e^{-iQ_{n}\left(z_{j_{n}}-z_{j_{n-1}}\right)}\,v_{j_{n}}(q_{n})\Bigg{]}J\left(E,\bm{p}_{in}\right)\,,$
(7)

where $J\left(E,\bm{p}_{in}\right)$ is the eikonal limit of the source
function with $p_{in}$ substituted from the corresponding pole,
$\theta_{j_{n},j_{n-1}}\equiv\theta\left(z_{j_{n}}-z_{j_{n-1}}\right)$, the
Landau-Pomeranchuk-Migdal (LPM) phase
$Q_{n}\equiv\frac{p_{n\perp}^{2}-p_{f\perp}^{2}}{2E}$ comes from the small
pole of the propagator $p_{n}^{2}=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.
Sievert:2018imd ; Sadofyev:2021ohn .

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
(II) 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 (II)
involves one or more $\theta(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

$\displaystyle(4E^{2})\int\frac{dq_{n,z}}{2\pi}\frac{dq_{n+1,z}}{2\pi}\frac{v_{j_{n}}(q_{n})v_{j_{n}}(q_{n+1})}{\left(p_{n}^{2}+i\epsilon\right)\left(p_{n+1}^{2}+i\epsilon\right)}e^{-iq_{n,z}\left(z_{j_{n}}-z_{j_{n-1}}\right)}e^{-iq_{n+1,z}\left(z_{j_{n}}-z_{j_{n-1}}\right)}\simeq$
$\displaystyle\hskip
28.45274pt-2iE\,\theta_{j_{n},j_{n-1}}\,e^{-i\left(q^{(p)}_{n,z}+q_{n+1,z}\right)\left(z_{j_{n}}-z_{j_{n-1}}\right)}\,\int\frac{dq_{n+1,z}}{2\pi}\frac{v_{j_{n}}\left(\bm{q}_{n},q^{(p)}_{n,z}\right)v_{j_{n}}(\bm{q}_{n+1},q_{n+1,z})}{p_{n+1}^{2}+i\epsilon}$
$\displaystyle\hskip
56.9055pt\simeq-\frac{1}{2}\,\theta_{j_{n},j_{n-1}}\,v_{j_{n}}(q^{2}_{n\perp})v_{j_{n}}(q^{2}_{n+1\perp})\,e^{-i\left(q^{(p)}_{n,z}+q_{n+1,z}\right)\left(z_{j_{n}}-z_{j_{n-1}}\right)}$
(8)

where $q^{(p)}_{n,z}=Q_{n}-\sum\limits_{m=n+1}^{N}q_{z,m}$ is the small pole
of $p_{n}^{2}=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 $\sum\limits_{m=n+2}^{N}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\left(\bm{q}_{n-1},q^{(p)}_{n-1,z}\right)$ 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.

$\sum_{i}f_{i}=\int d^{2}\bm{x}\,dz\>\rho(\bm{x},z)\>f(\bm{x},z)\,,$ (9)

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 $\bm{x}_{n}$
integral acts only on the corresponding Fourier factor, resulting in

$\int
d^{2}\bm{x}_{n}\,e^{-i(\bm{q}_{n}\pm\overline{\bm{q}}_{n})\cdot\bm{x}_{n}}=(2\pi)^{2}\,\delta^{(2)}(\bm{q}_{n}\pm\overline{\bm{q}}_{n})$

where $\bm{q}_{n}$ and $\overline{\bm{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
Gyulassy:2002yv  for the amplitude squared at $N$th order in opacity:

$\displaystyle\left\langle\left|M\right|^{2}\right\rangle^{(N)}$
$\displaystyle=\prod_{n=1}^{N}\left[(-1)\int\limits_{0}^{z_{n+1}}dz_{n}\int\frac{d^{2}\bm{q}_{n}}{(2\pi)^{2}}\,\mathcal{V}(\bm{q}_{n},z_{n})\right]\>\left|J\left(E,\bm{p}_{in}\right)\right|^{2}\,$
(10)

with

$\mathcal{V}(\bm{q},z)\equiv-\mathcal{C}\,\rho(z)\left(\left|v(q_{\perp}^{2})\right|^{2}-\delta^{(2)}(\bm{q})\int
d^{2}\bm{l}\,\left|v(l_{\perp}^{2})\right|^{2}\right)\,,$

where $\mathcal{V}(\bm{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,
$\mathcal{C}=\frac{\mathcal{C}_{\text{proj}}}{2C_{\bar{R}}}$ is the full color
factor with
$\mathcal{C}_{\text{proj}}\bm{1}=t^{a}_{\text{proj}}t^{a}_{\text{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 $\bm{x}$-space. Indeed, introducing the jet
distribution in the $\bm{x}$-space

$\displaystyle
E\frac{d\mathcal{N}}{d^{2}\bm{x}dE}=\frac{1}{2(2\pi)^{3}}\int\,\frac{d^{2}\bm{p}}{(2\pi)^{2}}\,\left\langle\left|M(\bm{p})\right|^{2}\right\rangle
e^{i\bm{p}\cdot\bm{x}}$ (11)

we readily write

$\displaystyle\frac{d\mathcal{N}}{d^{2}\bm{x}dE}=\sum_{N=0}^{\infty}\int\,\frac{d^{2}\bm{p}\,d^{2}\bm{r}}{(2\pi)^{2}}e^{i\bm{p}\cdot\left(\bm{x}-\bm{r}\right)}\,\frac{(-1)^{N}\left[\mathcal{V}(\bm{r})L\right]^{N}}{N!}\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{r}dE}=e^{-\mathcal{V}(\bm{x})L}\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{x}dE}$
(12)

where for simplicity we set $\rho(z)=const$, and, consequently,
$\mathcal{V}(\bm{q},z)=\mathcal{V}(\bm{q})$.

On the other hand, if the medium is inhomogeneous in the transverse
directions, then the $\bm{x}_{n}$ integrals cannot be simplified without
further assumptions. As in Sadofyev:2021ohn  we focus on the leading
corrections in the case of a slow $\bm{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

$\int
d^{2}\bm{x}_{n}\,x^{\alpha}_{n}\,e^{-i(\bm{q}_{n}\pm\overline{\bm{q}}_{n})\cdot\bm{x}_{n}}=i\,(2\pi)^{2}\,\frac{\partial}{\partial(q_{n}\pm\overline{q}_{n})_{\alpha}}\,\delta^{(2)}(\bm{q}_{n}\pm\overline{\bm{q}}_{n})$

where $\alpha$ is a vector index in the transverse 2D space.

In the absence of medium flow, the two medium parameters of interest are
$\rho$ and $\mu$, and to the leading order in gradients they can be written as

$\displaystyle\rho({\bm{x}},z)\approx\rho(z)+\bm{\nabla}\rho(z)\cdot{\bm{x}}\,,$
$\displaystyle\mu^{2}({\bm{x}},z)\approx\mu^{2}(z)+\bm{\nabla}\mu^{2}(z)\cdot{\bm{x}}\,,$
(13)

where for compactness we use $\rho(z)\equiv\rho(0,z)$. Integrating the
$\delta$-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 Sadofyev:2021ohn  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:

$\displaystyle\left\langle\left|M\right|^{2}\right\rangle^{(N)}$
$\displaystyle=\prod\displaylimits_{n=1}^{N}\left[\int\displaylimits_{0}^{z_{n+1}}\,dz_{n}\int\frac{d^{2}\bm{q}_{n}}{(2\pi)^{2}}\right]\left(1+\frac{1}{E}\sum_{m=1}^{N}(z_{m}-z_{m-1})\bm{p}_{m}\cdot\sum_{k=m}^{N}\hat{\bm{g}}_{k}\right)$
$\displaystyle\hskip
28.45274pt\times(-1)^{N}\mathcal{V}_{1}(\bm{q}_{1})\,...\,\mathcal{V}_{N}(\bm{q}_{N})|J(E,\bm{p}_{in})|^{2}\,,$
(14)

where we have introduced notations $\rho_{k}$, $\mu_{k}$, and
$\mathcal{V}_{k}$ to distinguish different sources of gradients,
$\hat{\bm{g}}_{k}\equiv\left(\bm{\nabla}\rho\frac{\delta}{\delta\rho_{k}}+\bm{\nabla}\mu^{2}\frac{\delta}{\delta\mu^{2}_{k}}\right)$
is an operator generating gradient contributions, and the ordering of the
$z$-integrals is enforced by the $\theta_{n,n-1}$ in (II) and (II). After the
variations with respect to $\rho_{k}$ and $\mu_{k}$ are performed, we again
set the thermodynamic parameters of the same type to be equal and constant in
$z$.

Further simplifying (II) and transforming to the $\bm{x}$-space, we find

$\displaystyle\frac{d\mathcal{N}^{(N)}}{d^{2}\bm{x}dE}$
$\displaystyle=\int\,\frac{d^{2}\bm{p}\,d^{2}\bm{r}}{(2\pi)^{2}}\,e^{i\bm{p}\cdot\left(\bm{x}-\bm{r}\right)}(-1)^{N}\big{[}\mathcal{V}(\bm{r})L\big{]}^{N}\Bigg{\\{}\frac{1}{N!}+\frac{L}{E(N+1)!}$
$\displaystyle\times\sum_{m=1}^{N}\Bigg{[}(N+1-m){\bm{p}}\cdot\left(\frac{\mathcal{V}^{\prime}(\bm{r})}{\mathcal{V}(\bm{r})}\bm{\nabla}\mu^{2}+\frac{1}{\rho}\bm{\nabla}\rho\right)+i(N+1-m)^{2}\frac{\bm{\nabla}\mathcal{V}(\bm{r})}{\rho\,\mathcal{V}(\bm{r})}\cdot\bm{\nabla}\rho$
$\displaystyle+i(N+1-m)\left(\frac{\bm{\nabla}\mathcal{V}^{\prime}(\bm{r})}{\mathcal{V}(\bm{r})}+(N-m)\frac{\mathcal{V}^{\prime}(\bm{r})}{\mathcal{V}(\bm{r})}\frac{\bm{\nabla}\mathcal{V}(\bm{r})}{\mathcal{V}(\bm{r})}\right)\cdot\bm{\nabla}\mu^{2}\Bigg{]}\Bigg{\\}}\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{r}dE}\,,$
(15)

and the opacity series can be again resummed, resulting in

$\displaystyle\frac{d\mathcal{N}}{d^{2}\bm{x}dE}$
$\displaystyle=e^{-\mathcal{V}(\bm{x})L}\Bigg{\\{}\left[1-i\frac{\mathcal{V}(\bm{x})L^{3}}{6E}\Bigg{(}\frac{\mathcal{V}^{\prime}(\bm{x})}{\mathcal{V}(\bm{x})}\bm{\nabla}\mu^{2}+\frac{1}{\rho}\bm{\nabla}\rho\Bigg{)}\cdot\bm{\nabla}\mathcal{V}(\bm{x})\right]\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{x}dE}$
$\displaystyle+i\frac{\mathcal{V}(\bm{x})L^{2}}{2E}\left(\frac{\mathcal{V}^{\prime}(\bm{x})}{\mathcal{V}(\bm{x})}\bm{\nabla}\mu^{2}+\frac{1}{\rho}\bm{\nabla}\rho\right)\cdot\bm{\nabla}\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{x}dE}\Bigg{\\}}\,,$
(16)

where
$\mathcal{V}^{\prime}(\bm{x})=\frac{\partial}{\partial\mu^{2}}\mathcal{V}(\bm{x})$.
Thus, we have derived (I), 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 Mehtar-
Tani:2006vpj ; McLerran:1993ka ; McLerran:1993ni . 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

$\displaystyle gA_{\mathrm{ext}}^{\mu a}(q)=-(2\pi)\,g^{\mu
0}\,v^{a}(q)\>\delta\left(q^{0}\right)\,,$ (17)

where the scattering potential should be set to
$v^{a}(q)=\sum_{i}e^{-i\vv{q}\cdot\vv{x}_{i}}\,t_{i}^{a}\,v_{i}(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 (II) in terms of an effective dressed propagator. For that,
we first write (II) using the in-medium color field (17), then

$\displaystyle
iM_{r}(p_{f})=\prod_{n=1}^{r}\left[(-1)\int\frac{d^{2}\bm{p}_{n}\,dp_{n,z}}{(2\pi)^{3}}\,t^{a}_{\mathrm{proj}}v^{a}(p_{n+1}-p_{n})\,\frac{2E}{p_{n}^{2}+i\epsilon}\right]\,J\left(p_{in}\right)\,,$
(18)

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

$\displaystyle iM_{r}(p_{f})$
$\displaystyle=\prod_{n=1}^{r}\Bigg{[}\int\frac{d^{2}\bm{p}_{n}\,d^{2}\bm{x}_{n}\,dz_{n}}{(2\pi)^{2}}\,\theta_{n,n-1}\,(-i)t^{a}_{\mathrm{proj}}v^{a}(\bm{x}_{n},z_{n})$
$\displaystyle\hskip 56.9055pt\times
e^{-i(\bm{p}_{n+1}-\bm{p}_{n})\cdot\bm{x}_{n}}e^{-i\frac{p_{n\perp}^{2}}{2E}(z_{n}-z_{n-1})}\Bigg{]}\,e^{i\frac{p_{f\perp}^{2}}{2E}z_{r}}\,J\left(E,\bm{p}_{in}\right)\,,$
(19)

where $(\bm{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 ${\bm{p}}_{1}={\bm{p}}_{in}$ and
${\bm{p}}_{r+1}={\bm{p}}_{f}$, we can write the perturbative amplitude as a
convolution between the source and a contribution to an effective single
particle propagator $G_{r}(\bm{p}_{f},L;\bm{p}_{in},0)$:

$\displaystyle
iM_{r}(p)=\int\frac{d^{2}{\bm{p}}_{in}}{(2\pi)^{2}}e^{i\frac{\bm{p}_{f}^{2}}{2E}L}\,G_{r}(\bm{p}_{f},L;\bm{p}_{in},0)\,J\left(E,\bm{p}_{in}\right)\,.$
(20)

Notice that the in-medium potential $v^{a}$ is screened by the Debye mass
$\mu$, and thus has a finite spatial support of size $\sim\frac{1}{\mu}$. 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

$\displaystyle
G(\bm{p}_{f},L;\bm{p}_{in},0)=\sum_{r=0}^{\infty}G_{r}(\bm{p}_{f},L;\bm{p}_{in},0)\,,$
(21)

where it is simple to check that in the case of vacuum propagation it reduces
to the usual Feynman result

$\displaystyle
G_{0}(\bm{p}_{f},L;\bm{p}_{in},0)=(2\pi)^{2}\delta^{(2)}(\bm{p}_{f}-\bm{p}_{in})e^{-i\frac{p_{f\perp}^{2}}{2E}L}\,.$
(22)

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 (III) and (20), which takes the usual Schrödinger-like
form

$\displaystyle\frac{\partial}{\partial L}G(\bm{p}_{f},L;\bm{p}_{in},0)$
$\displaystyle=-i\frac{p_{f\perp}^{2}}{2E}G(\bm{p}_{f},L;\bm{p}_{in},0)$
$\displaystyle+\int\frac{d^{2}\bm{l}d^{2}\bm{x}}{(2\pi)^{2}}\,(-i)t_{\text{proj}}^{a}\,v^{a}({\bm{x}},L)\,e^{-i(\bm{p}_{f}-\bm{l})\cdot{\bm{x}}}\,G(\bm{l},L;\bm{p}_{in},0)\,.$
(23)

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 (III) with the corresponding initial
condition is well known book:Kleinert_path_integrals , and its $\bm{x}$-space
form can be written as a path integral

$\displaystyle G({\bm{x}}_{L},L;{\bm{x}}_{0},0)$
$\displaystyle=\int\limits_{{\bm{x}}_{0}}^{{\bm{x}}_{L}}\mathcal{D}{\bm{r}}\exp\left(\frac{iE}{2}\int\limits_{0}^{L}d\tau\,\dot{{\bm{r}}}^{2}\right)$
$\displaystyle\hskip 28.45274pt\times{\cal
P}\exp\left(-i\int\limits_{0}^{L}d\tau\,t^{a}_{\text{proj}}v^{a}({\bm{r}}(\tau),\tau)\right)\,,$
(24)

where ${\cal P}$ indicates path ordering, and $\bm{x}_{L}\equiv\bm{x}(L)$ and
$\bm{x}_{0}\equiv\bm{x}(0)$ are the boundary conditions for the trajectory.
This effective propagator can be thought of as describing a massive non-
relativistic particle, moving from the initial position ${\bm{x}}_{0}$ at
"time" $\tau=0$ to the final position ${\bm{x}}_{L}$ at "time" $\tau=L$ in a
(random) potential $v^{a}(\bm{r},\tau)$.

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),

$\displaystyle\langle|M|^{2}\rangle=\int\frac{d^{2}{\bm{p}}_{in}d^{2}\overline{{\bm{p}}}_{in}}{(2\pi)^{4}}\langle
G({\bm{p}}_{f},L;{\bm{p}}_{in},0)G^{\dagger}({\bm{p}}_{f},L;\overline{{\bm{p}}}_{in},0)\rangle
J(E,{\bm{p}}_{in})J^{*}(E,\overline{{\bm{p}}}_{in})\,,$ (25)

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 in the BDMPS-Z approach. Since $G({\bm{x}}_{L},L;{\bm{x}}_{0},0)$
is a functional of $v^{a}(\bm{r}(\tau),\tau)$ , one first needs to consider
the average of the in-medium color fields. Since these are assumed to have
Gaussian statistics, only two-point functions of the potentials are non-
trivial. For instance, in the GW model, the corresponding average reads

$\displaystyle\langle
t^{a}_{\text{proj}}v^{a}({\bm{r}},\tau)t^{b}_{\text{proj}}v^{\dagger
b}(\overline{\bm{r}},\overline{\tau})\rangle=\mathcal{C}\,g^{4}\,\int
dz\,d^{2}\bm{x}\,\rho(\bm{x},\,z)$ $\displaystyle\hskip
28.45274pt\times\int\frac{d^{2}{\bm{q}}\,dq_{z}\>d^{2}\overline{\bm{q}}\,d\overline{q}_{z}}{(2\pi)^{6}}\frac{e^{i{\bm{q}}\cdot\left(\bm{r}-\bm{x}\right)}e^{-i\overline{\bm{q}}\cdot\left(\overline{\bm{r}}-\bm{x}\right)}e^{iq_{z}\left(\tau-z\right)}e^{-i\overline{q}_{z}\left(\overline{\tau}-z\right)}}{(q_{\perp}^{2}+q_{z}^{2}+\mu^{2}(\bm{x},z))(\overline{q}_{\perp}^{2}+\overline{q}^{2}_{z}+\mu^{2}(\bm{x},z))}\,,$
(26)

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 (III) when it is squared and
averaged, see e.g. Blaizot:2012fh . In fact, the average in (III) 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 (III) can be neglected without affecting the final result at the accuracy
level being considered. Using this simplification, we find

$\displaystyle\langle
t^{a}_{\text{proj}}v^{a}({\bm{r}},\tau)t^{b}_{\text{proj}}v^{\dagger
b}(\overline{\bm{r}},\overline{\tau})\rangle$
$\displaystyle\simeq\left(1+\frac{{\bm{r}}(\tau)+\overline{\bm{r}}(\tau)}{2}\cdot\hat{\bm{g}}\right)$
$\displaystyle\hskip
28.45274pt\times\mathcal{C}\,\delta(\tau-\overline{\tau})\,\rho\,g^{4}\int\frac{d^{2}\bm{q}}{(2\pi)^{2}}\frac{e^{i\bm{q}\cdot({\bm{r}}-\overline{\bm{r}})}}{(\bm{q}^{2}+\mu^{2})^{2}}\,,$
(27)

where $\bm{x}$ in $\rho(\bm{x},\,z)$ and $\mu^{2}(\bm{x},\,z)$ has been
replaced with $-i\frac{\partial}{\partial(\bm{q}-\overline{\bm{q}})}$ acting
on everything but the delta function $\delta^{(2)}(\bm{q}-\overline{\bm{q}})$,
and we assume that $\rho$ and $\mu^{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 $|{\bm{r}}-\overline{\bm{r}}|$, but also on its
center of mass transverse position $\frac{{\bm{r}}+\overline{\bm{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.

Figure 2: Diagrammatic illustration of (III). Here, the amplitude level
process is depicted in the upper half of the plot while the conjugate
amplitude part is depicted in the bottom one, with the dashed cut in between.

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 (III), we can write the average of two Wilson lines as444The
exponentiation follows from the fact that the averaging statistics is still
Gaussian.

$\displaystyle\left\langle{\cal
P}\exp\left(-i\int_{0}^{L}d\tau\,t_{\text{proj}}^{a}v^{a}({\bm{r}}(\tau),\tau)\right){\cal
P}\exp\left(i\int_{0}^{L}d\overline{\tau}\,t_{\text{proj}}^{b}v^{b}(\overline{\bm{r}}(\overline{\tau}),\overline{\tau})\right)\right\rangle$
$\displaystyle\hskip
28.45274pt=\exp\left\\{-\int\limits_{0}^{L}\,d\tau\,\left[1+\frac{\bm{r}(\tau)+\overline{\bm{r}}(\tau)}{2}\cdot\hat{\bm{g}}\right]\mathcal{V}\left(\bm{r}(\tau)-\overline{\bm{r}}(\tau)\right)\right\\}\,,$
(28)

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

$\displaystyle\left\langle
G({\bm{x}}_{L},L;{\bm{x}}_{0},0)G^{\dagger}(\overline{{\bm{x}}}_{L},L;\overline{{\bm{x}}}_{0},0)\right\rangle\equiv\left\langle
G({\bm{x}}_{L};{\bm{x}}_{0})G^{\dagger}(\overline{{\bm{x}}}_{L};\overline{{\bm{x}}}_{0})\right\rangle$
$\displaystyle\hskip
28.45274pt=\int\limits_{{\bm{x}}_{0}}^{{\bm{x}}_{L}}\mathcal{D}{\bm{r}}\int\limits_{\overline{{\bm{x}}}_{0}}^{\overline{{\bm{x}}}_{L}}\mathcal{D}\overline{\bm{r}}\,\exp\left\\{\frac{iE}{2}\int_{0}^{L}d\tau\,\left[\dot{{\bm{r}}}^{2}-\dot{\overline{\bm{r}}}^{2}\right]\right\\}$
$\displaystyle\hskip
56.9055pt\times\exp\left\\{-\int\limits_{0}^{L}\,d\tau\,\left[1+\frac{{\bm{r}}(\tau)+\overline{\bm{r}}(\tau)}{2}\cdot\hat{\bm{g}}\right]\mathcal{V}\left({\bm{r}}(\tau)-\overline{\bm{r}}(\tau)\right)\right\\}$
$\displaystyle\hskip
28.45274pt=\int\limits_{\bm{u}_{0}}^{\bm{u}_{L}}\mathcal{D}{\bm{u}}\int\limits_{\bm{w}_{0}}^{\bm{w}_{L}}\mathcal{D}\bm{w}\,\exp\left\\{\int_{0}^{L}d\tau\,\left[iE\,\dot{\bm{u}}\cdot\dot{\bm{w}}-\left(1+\bm{w}\cdot\hat{\bm{g}}\right)\mathcal{V}\left(\bm{u}(\tau)\right)\right]\right\\}\,,$
(29)

where $\bm{u}\equiv\bm{r}-\overline{\bm{r}}$ and
$\bm{w}\equiv\frac{\bm{r}+\overline{\bm{r}}}{2}$. The path integral is simple
enough to be evaluated analytically, since the $\bm{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 motion555One should notice that this EOM
involves a complex potential term, which in fact should be understood as a
perturbation around the trivial trajectory, see Appendix B. (EOM). Thus, using
the standard methods book:Kleinert_path_integrals ; Blaizot:2012fh ;
Apolinario:2014csa , we can reduce this expression to

$\displaystyle\left\langle
G({\bm{x}}_{L};{\bm{x}}_{0})G^{\dagger}(\overline{{\bm{x}}}_{L};\overline{{\bm{x}}}_{0})\right\rangle=\left(\frac{E}{2\pi
L}\right)^{2}\,\frac{\exp\left\\{iE\left(\bm{w}\cdot\dot{{\bm{u}}}_{c}\right)\Big{|}_{0}^{L}-\int\limits_{0}^{L}\,d\tau\,\mathcal{V}\left(\bm{u}_{c}(\tau)\right)\right\\}}{1+\frac{i}{EL}\hat{\bm{g}}\cdot\int\limits^{L}_{0}d\zeta\int\limits^{\zeta}_{0}d\xi\,\xi\,\bm{\nabla}\mathcal{V}\left({\bm{u}}_{c}(\xi)\right)}\,,$
(30)

where $\bm{u}_{c}(\tau)$ is the classical solution satisfying

$\displaystyle E\ddot{\bm{u}}=i\hat{\bm{g}}\,\mathcal{V}({\bm{u}}(\tau))\,,$
(31)

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 general, (31) does not admit an analytical
solution. However, since we are only interested in the leading gradient terms,
it can be solved in an expansion
${\bm{u}}_{c}={\bm{u}}_{c}^{(0)}+{\bm{u}}_{c}^{(1)}$, where
${\bm{u}}_{c}^{(0)}$ has no gradient dependence and ${\bm{u}}_{c}^{(1)}$ is
linear in gradients666Notice that the imaginary part of the classical EOM is
proportional to the gradient terms. Thus, substituting $\bm{u}_{c}(\tau)$ 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..

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"

$\displaystyle{\bm{u}}^{(0)}_{c}(\tau)=\frac{{\bm{u}}_{L}-{\bm{u}}_{0}}{L}\tau+{\bm{u}}_{0}\,,$
(32)

leading to the uniform broadening result in (12), see e.g. Blaizot:2012fh ;
Apolinario:2014csa . In turn, the leading gradient correction to the
trajectory is purely imaginary

$\displaystyle{\bm{u}}^{(1)}_{c}(\tau)$
$\displaystyle=\frac{i}{E}\,\hat{\bm{g}}\left\\{\int\limits^{\tau}_{0}\,d\zeta\,\int\limits_{0}^{\zeta}\,d\xi\,\mathcal{V}\left({\bm{u}}^{(0)}_{c}(\xi)\right)-\frac{\tau}{L}\int\limits^{L}_{0}\,d\zeta\,\int\limits_{0}^{\zeta}\,d\xi\,\mathcal{V}\left({\bm{u}}^{(0)}_{c}(\xi)\right)\right\\}\,,$
(33)

and satisfies the trivial boundary conditions
${\bm{u}}^{(1)}_{c}(L)={\bm{u}}^{(1)}_{c}(0)={\bm{0}}$.

The momentum scale ${\bm{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 (30) we find

$\displaystyle\langle
G(\bm{p}_{f},L;\bm{p}_{in},0)G^{\dagger}(\bm{p}_{f},L;\overline{\bm{p}}_{in},0)\rangle$
$\displaystyle\hskip 28.45274pt=\left(\frac{E}{2\pi
L}\right)^{2}\int\,d^{2}{\bm{u}}_{0}\,d^{2}\bm{w}_{0}\,d^{2}{\bm{u}}_{L}\,d^{2}\bm{w}_{L}\,e^{-i\bm{p}_{f}\cdot{\bm{u}}_{L}}e^{i{\bm{u}}_{0}\cdot\frac{\bm{p}_{in}+\overline{\bm{p}}_{in}}{2}+i\bm{w}_{0}\cdot\left(\bm{p}_{in}-\overline{\bm{p}}_{in}\right)}$
$\displaystyle\hskip
56.9055pt\times\frac{\exp\left\\{iE\left(\bm{w}\cdot\dot{{\bm{u}}}_{c}\right)\Big{|}_{0}^{L}-\int\limits_{0}^{L}\,d\tau\,\mathcal{V}\left(\bm{u}_{c}(\tau)\right)\right\\}}{1+\frac{i}{EL}\hat{\bm{g}}\cdot\int\limits^{L}_{0}d\zeta\int\limits^{\zeta}_{0}d\xi\,\xi\,\bm{\nabla}\mathcal{V}\left({\bm{u}}_{c}(\xi)\right)}$
$\displaystyle\hskip
28.45274pt=\frac{(2\pi)^{2}}{L^{2}}\int\,d^{2}{\bm{u}}_{0}\,d^{2}{\bm{u}}_{L}\,e^{-i\bm{p}_{f}\cdot{\bm{u}}_{L}}e^{i{\bm{u}}_{0}\cdot\frac{\bm{p}_{in}+\overline{\bm{p}}_{in}}{2}}$
$\displaystyle\hskip
56.9055pt\times\frac{\delta^{(2)}\left(\dot{{\bm{u}}}_{c}(L)\right)\delta^{(2)}\left(\bm{p}_{in}-\overline{\bm{p}}_{in}-E\dot{{\bm{u}}}_{c}(0)\right)}{1+\frac{i}{EL}\hat{\bm{g}}\cdot\int\limits^{L}_{0}d\zeta\int\limits^{\zeta}_{0}d\xi\,\xi\,\bm{\nabla}\mathcal{V}\left({\bm{u}}_{c}(\xi)\right)}\exp\left\\{-\int\limits_{0}^{L}\,d\tau\,\mathcal{V}\left(\bm{u}_{c}(\tau)\right)\right\\}\,.$
(34)

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

$\displaystyle\langle|M|^{2}\rangle\simeq\frac{1}{(2\pi
L)^{2}}\int\,d^{2}\bm{P}_{in}\,d^{2}{\bm{u}}_{0}\,d^{2}{\bm{u}}_{L}\,e^{-i\bm{p}_{f}\cdot{\bm{u}}_{L}}e^{i\bm{P}_{in}\cdot{\bm{u}}_{0}}$
$\displaystyle\hskip
113.81102pt\times\frac{\delta^{(2)}\left(\dot{{\bm{u}}}_{c}(L)\right)\exp\left\\{-\int\limits_{0}^{L}\,d\tau\,\mathcal{V}\left(\bm{u}_{c}(\tau)\right)\right\\}}{1+\frac{i}{EL}\hat{\bm{g}}\cdot\int\limits^{L}_{0}d\zeta\int\limits^{\zeta}_{0}d\xi\,\xi\,\bm{\nabla}\mathcal{V}\left({\bm{u}}_{c}(\xi)\right)}\left|J\left(E,\bm{P}_{in}\right)\right|^{2}\,,$
(35)

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
$\bm{P}_{in}=\frac{\bm{p}_{in}+\overline{\bm{p}}_{in}}{2}$. We have also used
that the leading shift of $\bm{p}_{in}-\overline{\bm{p}}_{in}$ in
$J\left(\bm{p}_{in}\right)J\left(\overline{\bm{p}}_{in}\right)$ is zero at the
first order in gradients, since the first derivative of this symmetric
function with respect to $\bm{p}_{in}-\overline{\bm{p}}_{in}$ is zero at
$\bm{p}_{in}=\overline{\bm{p}}_{in}$.

Thus, we can write the corresponding Fourier transformed distribution as

$\displaystyle\frac{d\mathcal{N}}{d^{2}\bm{x}dE}\simeq\frac{1}{L^{2}}\int\,d^{2}{\bm{u}}_{0}\,d^{2}{\bm{u}}_{L}\,\delta^{(2)}(\bm{x}-\bm{u}_{L})$
$\displaystyle\hskip
56.9055pt\times\frac{\delta^{(2)}\left(\dot{{\bm{u}}}_{c}(L)\right)\exp\left\\{-\int\limits_{0}^{L}\,d\tau\,\mathcal{V}\left(\bm{u}_{c}(\tau)\right)\right\\}}{1+\frac{i}{EL}\hat{\bm{g}}\cdot\int\limits^{L}_{0}d\zeta\int\limits^{\zeta}_{0}d\xi\,\xi\,\bm{\nabla}\mathcal{V}\left({\bm{u}}_{c}(\xi)\right)}\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{u}_{0}dE}\,.$
(36)

The velocity constraint at $\tau=L$ allows to simplify the form of
${\bm{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 Blaizot:2015lma . In the present case, this constraint
implies

$\displaystyle\frac{{\bm{u}}_{L}-{\bm{u}}_{0}}{L}+\frac{i}{E}\,\hat{\bm{g}}\Bigg{\\{}\int\limits_{0}^{L}\,d\xi\,\mathcal{V}\left({\bm{u}}^{(0)}_{c}(\xi)\right)-\frac{1}{L}\int\limits^{L}_{0}\,d\zeta\,\int\limits_{0}^{\zeta}\,d\xi\,\mathcal{V}\left({\bm{u}}^{(0)}_{c}(\xi)\right)\Bigg{\\}}=0\,,$
(37)

and one can use it to relate ${\bm{u}}_{0}$ and ${\bm{u}}_{L}$.

At the leading order in gradients we find

$\displaystyle{\bm{u}}_{L}-{\bm{u}}_{0}+\frac{iL^{2}}{2E}\,\hat{\bm{g}}\,\mathcal{V}\left(\bm{u}_{L}\right)\simeq
0\,.$ (38)

and, consequently, ${\bm{u}}_{c}$ takes a remarkably simple form

$\displaystyle{\bm{u}}_{c}(\tau)={\bm{u}}_{L}+\frac{i}{E}\,\hat{\bm{g}}\,\mathcal{V}\left({\bm{u}}_{L}\right)\left\\{\frac{(\tau-L)^{2}}{2}\right\\}\,.$
(39)

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 $\mathcal{V}(\bm{x})$. One should also notice that the
delta function constraining the velocity at $\tau=L$ has a non-trivial
dependence on $\bm{u}_{0}$ which enters the argument also through
$\mathcal{V}\left({\bm{u}}^{(0)}_{c}(\xi)\right)$, leading to an additional
multiple:

$\displaystyle
L^{2}\left[1-\frac{iL}{E}\hat{\bm{g}}\cdot\Bigg{\\{}\int\limits_{0}^{L}\,d\xi-\frac{1}{L}\int\limits^{L}_{0}\,d\zeta\,\int\limits_{0}^{\zeta}\,d\xi\Bigg{\\}}\,\left(1-\frac{\xi}{L}\right)\,\bm{\nabla}\mathcal{V}\left({\bm{u}}^{(0)}_{c}(\xi)\right)\right]^{-1}\,.$
(40)

Combining all the previous results and expanding in the smallness of
gradients, we conclude that the $\bm{x}$-space form of the jet distribution
reads

$\displaystyle\frac{d\mathcal{N}}{d^{2}\bm{x}dE}$
$\displaystyle\simeq\exp\left\\{-\mathcal{V}\left(\bm{x}\right)L\right\\}\Bigg{\\{}\left[1-\frac{iL^{3}}{6E}\bm{\nabla}\mathcal{V}\left(\bm{x}\right)\cdot\hat{\bm{g}}\,\mathcal{V}\left(\bm{x}\right)\right]\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{x}dE}$
$\displaystyle\hskip
56.9055pt+\frac{iL^{2}}{2E}\,\hat{\bm{g}}\,\mathcal{V}\left(\bm{x}\right)\cdot\bm{\nabla}\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{x}dE}\Bigg{\\}}\,,$
(41)

and one can easily see that it exactly coincides with (II), 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 $\bm{w}$ dependence in (III) would no longer be
linear).

Before we turn to the properties of the distribution (III), 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

$\displaystyle\frac{d\mathcal{N}}{d^{2}\bm{x}dE}\simeq\exp\left\\{-\int_{0}^{L}\,d\tau\,\mathcal{V}\left(\bm{x},\tau\right)\right\\}$
$\displaystyle\hskip
14.22636pt\times\left\\{\rule{0.0pt}{28.45274pt}\right.\left[1-\frac{i}{E}\int\limits_{0}^{L}\,d\tau\,\bm{\nabla}\mathcal{V}\left(\bm{x},\tau\right)\cdot\left(\int\limits_{L}^{\tau}\,d\zeta\,\int\limits_{0}^{\zeta}\,d\xi+(L-\tau)\int\limits^{L}_{0}\,d\xi\right)\hat{\bm{g}}(\xi)\mathcal{V}\left(\bm{x},\xi\right)\right]$
$\displaystyle\hskip
28.45274pt+\frac{i}{E}\int\limits_{0}^{L}d\zeta\int\limits^{\zeta}_{0}\,d\xi\,\hat{\bm{g}}\left(\xi\right)\mathcal{V}\left(\bm{x},\xi\right)\cdot\bm{\nabla}\left.\rule{0.0pt}{28.45274pt}\right\\}\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{x}dE}\,,$
(42)

where
$\hat{\bm{g}}\left(\tau\right)=\left(\bm{\nabla}\rho(\tau)\frac{\delta}{\delta\rho}+\bm{\nabla}\mu^{2}(\tau)\frac{\delta}{\delta\mu^{2}}\right)$.

##  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 (II) (or equivalently (III)). 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
non-zero 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 (II) and illustrate our results with simple numerical
estimates for the broadening probability, using phenomenologically relevant
parameters.

Figure 3: An illustration of a single parton evolving in the presence of an
inhomogeneous static slab of matter. The initial jet direction gets modified
due to the presence of the temperature gradient. This modification depends on
the angle $\theta$ between $\bm{\nabla}T$ and the transverse momentum
accumulated during the evolution, as shown by the dashed shape.

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
$\bm{x}$-space distribution at $\bm{x}={\bm{0}}$, which corresponds to the
volume integral of the momentum-space distribution. At this point, the
$\bm{x}$-space dipole potential is zero $\mathcal{V}(\bm{0})=0$ as well as its
transverse gradient $\bm{\nabla}\mathcal{V}(\bm{0})={\bm{0}}$, and the full
distribution (I) is equal to
$\int\,d^{2}\bm{p}\,\frac{dN^{(0)}}{d^{2}\bm{p}\,dE}$. Thus, the unitarity
condition is satisfied by the broadened jet momentum distribution even in an
inhomogeneous medium

$\displaystyle\int\,d^{2}\bm{p}\,\frac{d\mathcal{N}}{d^{2}\bm{p}dE}=\frac{d\mathcal{N}}{d^{2}\bm{x}dE}\Bigg{|}_{{\bm{x}}=0}=\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{x}dE}\Bigg{|}_{{\bm{x}}=0}=\int\,d^{2}\bm{p}\,\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{p}dE}\,.$
(43)

One could also notice that the distribution (III) is not completely factorized
anymore, as it is in the uniform case. The gradient effects enter not only
through the probe-medium interactions, but also through a shift of the
argument in the initial hard distribution $\frac{dN^{(0)}}{d^{2}\bm{p}\,dE}$.
Thus, the relation between the final and initial distributions involves an
operator

$\displaystyle\frac{d\mathcal{N}}{d^{2}\bm{x}dE}$ $\displaystyle={\cal
P}\left(\bm{x}\right)\hat{\mathsf{S}}\left(\bm{x}\right)\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{x}dE}\,,$
(44)

where ${\cal P}({\bm{x}})$ is a transform of ${\cal P}({\bm{p}})$, the
(broadening) probability for a particle of energy $E$ to acquire transverse
momentum ${\bm{p}}$ due to propagation in the medium for a distance $L$, and
$\hat{\mathsf{S}}\left(\bm{x}\right)$ is a shift operator being an identity
operator in the absence of gradients. One can immediately see that ${\cal
P}({\bm{p}})$ is self-normalized, and that
$\hat{\mathsf{S}}\left(\bm{x}\right)$ 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:

$\displaystyle\left\langle
F(\bm{p})\right\rangle=\frac{\int\,d^{2}\bm{p}\,F(\bm{p})\,\frac{d\mathcal{N}}{d^{2}\bm{p}dE}}{\int\,d^{2}\bm{p}\,\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{p}dE}}\,,$
(45)

where $F(\bm{p})$ is an arbitrary function. Following Sadofyev:2021ohn , we
focus on a Gaussian initial distribution

$\displaystyle E\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{p}\,dE}$
$\displaystyle=\frac{f(E)}{2\pi w^{2}}e^{-\frac{p_{\perp}^{2}}{2w^{2}}}\,,$
(46)

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 $\langle
p_{\perp}^{2k}\rangle=\langle
p_{\perp}^{2k}\rangle_{\nabla\rho=\nabla\mu^{2}={\bm{0}}}$, and we may focus
on the odd moments $\langle{\bm{p}}\>p_{\perp}^{2k}\rangle$ sensitive to the
anisotropy777For an earlier attempt to study the matter anisotropy with
directional jet quenching, see Majumder:2006wi ., see e.g. Sadofyev:2021ohn ;
Antiporda:2021hpk .

Starting with the simplest case of the averaged momentum
$\left\langle\bm{p}\right\rangle$ at the given energy, we find it convenient
to use the opacity expanded form of the distribution (II) rather than its
resummed form. Doing so, one finds that the moment is controlled by the
following integral

$\displaystyle\frac{1}{2\pi
w^{2}}\int\,d^{2}\bm{p}\,p^{\alpha}\,p_{m}^{\beta}\,e^{-\frac{p_{in\perp}^{2}}{2w^{2}}}$
$\displaystyle=\delta^{\alpha\beta}w^{2}+\sum_{i=1}^{N}q^{\alpha}_{i}\sum_{j=1}^{m-1}q^{\beta}_{j}\,,$
(47)

where $\bm{p}_{n}=\bm{p}_{f}-\sum\limits_{m=n}^{N}\bm{q}_{m}$,
$\bm{p}_{in}=\bm{p}_{1}$, and $\alpha$ and $\beta$ 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 (II) gives zero. Thus, for the
contribution at the $N$th order in opacity, we find

$\displaystyle\left\langle p^{\alpha}\right\rangle^{(N)}$
$\displaystyle=\frac{L^{N+1}}{E\,(N+1)!}\prod\displaylimits_{n=1}^{N}\left[\int\frac{d^{2}q_{n}}{(2\pi)^{2}}\right]\sum_{m=1}^{N}\left[\left(\sum_{i=1}^{N}q^{\alpha}_{i}\sum_{j=1}^{m-1}q^{\beta}_{j}\right)\sum_{k=m}^{N}\hat{g}^{\beta}_{k}\right]$
$\displaystyle\hskip
56.9055pt\times(-1)^{N}\mathcal{V}_{1}(\bm{q}_{1})\,...\,\mathcal{V}_{N}(\bm{q}_{N})\,,$
(48)

where we have additionally used that
$\int\,d^{2}\bm{q}_{n}\,\mathcal{V}_{n}(\bm{q}_{n})=0$ removing the $w^{2}$
term. For any $N>1$ this expression is zero due to the same property of
$\mathcal{V}_{n}(\bm{q}_{n})$, since the integrand involves only two powers of
$\bm{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 $\left\langle\bm{p}\right\rangle={\bm{0}}$ to all orders in
opacity.

The first non-zero moment of the momentum broadening distribution at the first
order in opacity is $\left\langle\bm{p}\>p_{\perp}^{2}\right\rangle$
Sadofyev:2021ohn . We can again perform the final momentum averaging, and find

$\displaystyle\frac{1}{2\pi
w^{2}}\int\,d^{2}\bm{p}\,p^{\alpha}\,p_{m}^{\beta}\,p_{\perp}^{2}\,e^{-\frac{\bm{p}_{in}^{2}}{2w^{2}}}$
$\displaystyle=\left[4w^{4}+w^{2}\left(\sum_{i=1}^{N}\bm{q}_{i}\right)^{2}\right]\delta^{\alpha\beta}+4w^{2}\left(\sum_{i=1}^{N}q^{\alpha}_{i}\right)\left(\sum_{j=1}^{m-1}q^{\beta}_{j}\right)$
$\displaystyle+2w^{2}\left(\sum_{i=1}^{N}q^{\alpha}_{i}\right)\left(\sum_{j=1}^{N}q^{\beta}_{j}\right)+\left(\sum_{i=1}^{N}q^{\alpha}_{i}\right)\left(\sum_{j=1}^{m-1}q^{\beta}_{j}\right)\left(\sum_{l=1}^{N}\bm{q}_{l}\right)^{2}\,.$
(49)

Thus, one can see that the only non-zero contributions may come at $N=1$

$\displaystyle\left\langle p^{\alpha}\>p_{\perp}^{2}\right\rangle^{(1)}$
$\displaystyle=-\frac{w^{2}L^{2}}{E}\,\hat{g}^{\alpha}\int\frac{d^{2}\bm{q}}{(2\pi)^{2}}\,q_{\perp}^{2}\,\mathcal{V}(\bm{q})\,$
(50)

and at $N=2$

$\displaystyle\left\langle p^{\alpha}\>p_{\perp}^{2}\right\rangle^{(2)}$
$\displaystyle=\frac{L^{3}}{6E}\int\frac{d^{2}\bm{q}_{1}}{(2\pi)^{2}}\,q_{1\perp}^{2}\,\mathcal{V}(\bm{q}_{1})\,\hat{g}^{\alpha}\int\frac{d^{2}\bm{q}_{2}}{(2\pi)^{2}}\,q_{2\perp}^{2}\,\mathcal{V}(\bm{q}_{2})\,,$
(51)

while the higher orders in opacity decouple. In the case of the GW potential,
the terms proportional to $\bm{\nabla}\rho$ are divergent, and should be
regularized at some scale. Comparing with Sadofyev:2021ohn , we set
$|\bm{q}_{\text{max}}|=\sqrt{E\mu}$, and find

$\displaystyle\left\langle p^{\alpha}\>p_{\perp}^{2}\right\rangle$
$\displaystyle=\frac{w^{2}L^{2}\mu^{2}}{E\,\lambda}\frac{\nabla^{\alpha}\rho}{\rho}\ln\frac{E}{\mu}+\frac{L^{3}\mu^{4}}{6E\,\lambda^{2}}\frac{\nabla^{\alpha}\rho}{\rho}\left(\ln\frac{E}{\mu}\right)^{2}\,,$
(52)

where $\lambda=\frac{1}{\rho\sigma_{0}}$ with
$\sigma_{0}=\int\frac{d^{2}\bm{q}}{(2\pi)^{2}}\,\mathcal{C}\,|v(q_{\perp}^{2})|^{2}$.
The first term in the expression above precisely agrees with the result
obtained in Sadofyev:2021ohn . In turn, the second term is new and indicates
that the higher $N$ contributions to the transverse momentum moments with
integer $k\geq 1$ are generally dominating as long as the potential integrals
are divergent. Such 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 ${\cal P}(\bm{p})$, which has multiple phenomenological
applications Caucal:2019uvr ; Blaizot:2013vha ; Kutak:2018dim ; Barata:2021byj
; Soudi:2021aar . Its Fourier transform can be obtained from (III) if we set
$E\frac{d\mathcal{N}^{(0)}}{d^{2}\bm{x}dE}$ to be a constant in the transverse
directions, decoupling in this way the shift operator
$\hat{\mathsf{S}}(\bm{x})$:

$\displaystyle{\cal P}(\bm{p})$
$\displaystyle=\int\,d^{2}\bm{x}\,e^{-i\bm{p}\cdot\bm{x}}e^{-\mathcal{V}\left(\bm{x}\right)L}\left[1-\frac{iL^{3}}{6E}\bm{\nabla}\mathcal{V}\left(\bm{x}\right)\cdot\hat{\bm{g}}\mathcal{V}\left(\bm{x}\right)\right]\,.$
(53)

At this point one has to specify the particular model for the medium, fixing
the dipole potential, and evaluate ${\cal P}(\bm{p})$ explicitly. We consider
the GW model, and in this case the potential reads

$\mathcal{V}^{\rm
GW}(x_{\perp})=\frac{\mathcal{C}g^{4}\rho}{4\pi\mu^{2}}\left(1-\mu\,x_{\perp}K_{1}(\mu\,x_{\perp})\right)\,,$
(54)

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 ${\cal P}(\bm{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:
$\frac{p_{\perp}}{\mu}$ and
$\chi\equiv\frac{\mathcal{C}g^{4}\rho}{4\pi\mu^{2}}L$, where $\chi$ is the
opacity in the GW model. The dipole potential $\mathcal{V}(x_{\perp})$ tends
to a constant at large $\mu\,x_{\perp}$, and the large $x_{\perp}$ limit of
the integration is controlled by the fast oscillation of the Fourier
exponential. For sufficiently large momenta such that $p_{\perp}\gg\mu$, one
can expand the dipole potential in powers of $\mu\,x_{\perp}$, and for the GW
model it reads

$\displaystyle\frac{4L}{\chi}\mathcal{V}^{\rm
GW}(x_{\perp})=\mu^{2}x_{\perp}^{2}\log\frac{4e^{1-2\gamma_{E}}}{\mu^{2}x_{\perp}^{2}}+\mathcal{O}\left(\mu^{4}x^{4}_{\perp}\right)\,,$
(55)

where $\gamma_{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. Barata:2020sav ; Barata:2020rdn ; Djordjevic:2009cr .

If we further require that $p_{\perp}^{2}\gg\chi\,\mu^{2}$, then the
exponential of the dipole potential itself can be expanded, and, after
expanding the in-medium potentials, we find

$\displaystyle{\cal P}(\bm{p})$
$\displaystyle\simeq\int\,d^{2}\bm{x}\,e^{-i\bm{p}\cdot\bm{x}}\Bigg{[}1-\mathcal{V}\left(x_{\perp}\right)L+\frac{1}{2}\mathcal{V}^{2}\left(x_{\perp}\right)L^{2}-\frac{iL^{3}}{6E}\bm{\nabla}\mathcal{V}\left(x_{\perp}\right)\cdot\hat{\bm{g}}\mathcal{V}\left(x_{\perp}\right)\Bigg{]}$
$\displaystyle\simeq\frac{4\pi\mu^{2}\chi}{p_{\perp}^{4}}+\frac{16\pi\mu^{4}\chi^{2}}{p_{\perp}^{6}}\left(\log\frac{p_{\perp}^{2}}{\mu^{2}}-2\right)$
$\displaystyle\hskip
113.81102pt+\frac{4\pi\mu^{4}\chi^{2}L}{3E}\left[\frac{\bm{\nabla}\rho}{\rho}\left(\log\frac{p_{\perp}^{4}}{\mu^{4}}-4\right)-\frac{\bm{\nabla}\mu^{2}}{\mu^{2}}\right]\cdot\frac{\bm{p}}{p_{\perp}^{6}}\,,$
(56)

where we have omitted terms proportional to $\delta^{(2)}(\bm{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_{\perp}$ and can only be generated in an inhomogeneous medium. In
fact, one can compare the asymptotic structure of (IV) directly with the first
non-trivial moment (52) in the limit of infinitely narrow source (i.e. $w\to
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 (IV), one may notice that its dominant contribution to
$\langle{\bm{p}}\>p_{\perp}^{2}\rangle$ scales as
$\left|\frac{\bm{\nabla}\rho}{E\rho}\right|\int
dp_{\perp}\frac{\log\,p_{\perp}}{p_{\perp}}$, and with a similar UV
regularization it gives
$\left|\frac{\bm{\nabla}\rho}{E\rho}\right|\left(\log\frac{E}{\mu}\right)^{2}$
in qualitative agreement with the previous result.

Expanding the probability distribution ${\cal P}$ in two parameters, one has
to pay a particular attention to the applicability of the result. Indeed, if
the large momentum expansion parameter $\frac{\mu}{p_{\perp}}$ 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 (IV) 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 Barata:2020rdn
; Caucal:2019uvr  corresponds to the intermediate momentum region $\mu^{2}\ll
p_{\perp}^{2}\leq\chi\,\mu^{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 $\bm{x}$ can be neglected. Indeed, one can introduce a new scale
$\mathcal{Q}^{2}\gg\mu^{2}$ such that

$\displaystyle\frac{4L}{\chi}\mathcal{V}^{\rm
GW}(x_{\perp})\simeq\mu^{2}x_{\perp}^{2}\left(\log\frac{\mathcal{Q}^{2}}{\mu^{2}}+\log\frac{4e^{1-2\gamma_{E}}}{\mathcal{Q}^{2}x_{\perp}^{2}}\right)\simeq\mu^{2}x_{\perp}^{2}\log\frac{\mathcal{Q}^{2}}{\mu^{2}}\,.$
(57)

In practice, it can be fixed up to an overall coefficient based on
phenomenological arguments, see e.g. Barata:2020rdn ; Barata:2021wuf  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. Caucal:2019uvr ; Casalderrey-Solana:2014bpa ;
Blaizot:2013vha , and often referred to as BDMPS-Z/ASW approximation
Zakharov:1997uu ; Armesto:2003jh . Taking into account the gradient effects,
we find

$\displaystyle{\cal P}(\bm{p})$
$\displaystyle\simeq\int\,d^{2}\bm{x}\,e^{-i\bm{p}\cdot\bm{x}}e^{-\frac{1}{4}\chi\mu^{2}x_{\perp}^{2}\log\frac{\mathcal{Q}^{2}}{\mu^{2}}}\left[1-\frac{i\mu^{4}\chi^{2}L}{48E}\log\frac{\mathcal{Q}^{2}}{\mu^{2}}\,x_{\perp}^{2}\,\bm{x}\cdot\left(\frac{\bm{\nabla}\rho}{\rho}\log\frac{\mathcal{Q}^{2}}{\mu^{2}}-\frac{\bm{\nabla}\mu^{2}}{\mu^{2}}\right)\right]$
$\displaystyle=\frac{4\pi}{\chi\mu^{2}\log\frac{\mathcal{Q}^{2}}{\mu^{2}}}\Bigg{[}1+\frac{L}{6E}\,\frac{p_{\perp}^{2}-2\chi\mu^{2}\log\frac{\mathcal{Q}^{2}}{\mu^{2}}}{\chi\mu^{2}\log\frac{\mathcal{Q}^{2}}{\mu^{2}}}\,\left(\frac{\bm{\nabla}\rho}{\rho}-\frac{1}{\log\frac{\mathcal{Q}^{2}}{\mu^{2}}}\frac{\bm{\nabla}\mu^{2}}{\mu^{2}}\right)\cdot\bm{p}\Bigg{]}e^{-\frac{p_{\perp}^{2}}{\chi\mu^{2}\log\frac{\mathcal{Q}^{2}}{\mu^{2}}}}\,.$
(58)

The usual Gaussian result gets an overall ${\bm{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}_{\perp}\sim\chi\mu^{2}\log\frac{\mathcal{Q}^{2}}{\mu^{2}}$, roughly at
the peak of the $p_{\perp}$-distribution Barata:2020rdn . In this region, the
leading behavior of ${\cal 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 $\mu^{2}\ll p_{\perp}^{2}\leq\chi\,\mu^{2}$, we can also consider
the small logarithmic correction to the dipole potential, neglected in (57),
as a perturbation

$\displaystyle\mathcal{V}^{\rm
GW}(x_{\perp})\simeq\frac{\chi\mu^{2}}{4L}x_{\perp}^{2}\left(\log\frac{\mathcal{Q}^{2}}{\mu^{2}}+\log\frac{4e^{1-2\gamma_{E}}}{\mathcal{Q}^{2}x_{\perp}^{2}}\right)=\mathcal{V}_{0}(x_{\perp})+\delta\mathcal{V}(x_{\perp})\,.$
(59)

This approach to capture the leading effects beyond the quadratic Gaussian
approximation is known as the improved opacity expansion (IOE), see e.g.
Mehtar-Tani:2019ygg ; Mehtar-Tani:2019tvy ; Barata:2020sav . Substituting the
potential into the broadening probability, we write it as an expansion

$\displaystyle{\cal P}_{\rm IOE}(\bm{p})$ $\displaystyle={\cal
P}_{0+1}(\bm{p})+{\cal P}_{1+0}(\bm{p})+{\cal P}_{2+0}(\bm{p})+{\cal
P}_{1+1}(\bm{p})+{\cal P}_{2+1}(\bm{p})+{\cal P}_{3+1}(\bm{p})+...\,,$ (60)

where the first number in the subscript corresponds to the order in
$\delta\mathcal{V}$, and the second one counts the order in gradients. The
leading contribution ${\cal P}_{0+1}(\bm{p})$ is given by (IV), the term
${\cal P}_{1+0}(\bm{p})$ has been discussed in details in Barata:2020rdn , 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
$\delta\mathcal{V}$ give the correct large momentum limit (IV) even in the
presence of gradients, c.f. Barata:2020rdn .

All the $\bm{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

$\displaystyle I_{1}(\bm{p},a)$
$\displaystyle=\int\,d^{2}\bm{x}\,e^{-i\bm{p}\cdot\bm{x}}e^{-ax_{\perp}^{2}}=\frac{\pi}{a}e^{-\frac{p_{\perp}^{2}}{4a}}\,,$
$\displaystyle I_{2}(\bm{p},a,b)$
$\displaystyle=\int\,d^{2}\bm{x}\,e^{-i\bm{p}\cdot\bm{x}}e^{-ax_{\perp}^{2}}\,\log\frac{1}{b^{2}x_{\perp}^{2}}=\frac{\pi}{a}e^{-\frac{p_{\perp}^{2}}{4a}}\left[Ei\left(\frac{p_{\perp}^{2}}{4a}\right)-\log\frac{b^{2}p_{\perp}^{2}}{4a^{2}}\right]$
$\displaystyle\simeq\frac{4\pi}{p_{\perp}^{2}}+\frac{16\pi
a}{p_{\perp}^{4}}+\mathcal{O}\left(\frac{a^{2}}{p_{\perp}^{6}}\right)\,,$
$\displaystyle I_{3}(\bm{p},a,b)$
$\displaystyle=\int\,d^{2}\bm{x}\,e^{-i\bm{p}\cdot\bm{x}}e^{-ax_{\perp}^{2}}\,\log^{2}\frac{1}{b^{2}_{\perp}x_{\perp}^{2}}$
$\displaystyle\simeq\frac{16\pi}{p_{\perp}^{2}}\left(\gamma_{E}+\log\frac{p_{\perp}}{2b}\right)+\frac{64\pi
a}{p_{\perp}^{4}}\left(\gamma_{E}-1+\log\frac{p_{\perp}}{2b}\right)+\mathcal{O}\left(\frac{a^{2}}{p_{\perp}^{6}}\log\frac{p_{\perp}}{2b}\right)\,,$
(61)

where we omit the explicit form of $I_{3}$ for brevity. In the large
$p_{\perp}$ limit, the momentum derivatives acting on a term reduce its
contribution removing powers of $p_{\perp}$. Consequently, one can compare the
relative importance of the different contributions in (60) by counting the
powers of $\bm{x}$ entering through $\mathcal{V}$ and $\delta\mathcal{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 $\delta\mathcal{V}$ to be sure that we keep all the
relevant contributions.

We first notice that the leading term in the IOE expansion ${\cal
P}_{0+1}(\bm{p})$ is exponentially suppressed for
$p_{\perp}^{2}\gg\chi\mu^{2}$, while ${\cal
P}_{1+0}(\bm{p})\propto\frac{1}{p_{\perp}^{4}}$ gives the leading Coulomb term
in this limit. Similarly, the sub-leading non-gradient contribution satisfies
${\cal P}_{2+0}(\bm{p})\propto\frac{1}{p_{\perp}^{6}}\log p_{\perp}$. The
mixed gradient terms in turn scale as ${\cal
P}_{1+1}(\bm{p})\propto\frac{1}{p_{\perp}^{5}}$ and ${\cal
P}_{2+1}(\bm{p})\propto\frac{1}{p_{\perp}^{5}}\log p_{\perp}$, while one can
check that the next contribution ${\cal
P}_{3+1}(\bm{p})\propto\frac{1}{p_{\perp}^{6}}\log^{2}p_{\perp}$ is suppressed
both by a higher power of inverse momentum as well as by the gradients, and
can be omitted. Finally, we find

$\displaystyle{\cal P}_{\rm IOE}(\bm{p})$
$\displaystyle\simeq\frac{4\pi\mu^{2}\chi}{p_{\perp}^{4}}+\frac{16\pi\mu^{4}\chi^{2}}{p_{\perp}^{6}}\left(\log\frac{p_{\perp}^{2}}{\mu^{2}}-2\right)$
$\displaystyle\hskip
28.45274pt+\frac{4\pi\mu^{4}\chi^{2}L}{3E}\left[\frac{\bm{\nabla}\rho}{\rho}\left(\log\frac{p_{\perp}^{4}}{\mu^{4}}-4\right)-\frac{\bm{\nabla}\mu^{2}}{\mu^{2}}\right]\cdot\frac{\bm{p}}{p_{\perp}^{6}}\,,$
(62)

which precisely agrees with the large momentum limit (IV), showing that the
IOE covers both limits considered above. It could be also noticed that the
logarithmic terms in ${\cal P}_{\rm IOE}(\bm{p})$ are sensitive to a
cancellation of $\mathcal{Q}$ between terms of different orders in the IOE,
c.f. Barata:2021wuf .

Finally, to have a more quantitative understanding of the broadening
distribution ${\cal P}({\bm{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 $\mu\propto gT$
and $\rho\propto T^{3}$, then

$\displaystyle\frac{\bm{\nabla}\rho}{\rho}=3\frac{\bm{\nabla}T}{T}\,,\quad\frac{\bm{\nabla}\mu^{2}}{\mu^{2}}=2\frac{\bm{\nabla}T}{T}\,.$
(63)

The probability ${\cal P}({\bm{p}})$ is now a function of $\theta$, the angle
between ${\bm{p}}$ and $\bm{\nabla}T$, and a dimensionless quantity
$c_{T}\equiv\left|\frac{\bm{\nabla}T}{ET}\right|\ll 1$, see Fig. 3. Performing
the remaining angular integration, we can write ${\cal P}$ as

$\displaystyle{\cal P}({\bm{p}})$
$\displaystyle=2\pi\int\limits^{\infty}_{0}dx_{\perp}\,x_{\perp}\,e^{-\mathcal{V}^{\rm
GW}(x_{\perp})L}\Bigg{\\{}J_{0}(p_{\perp}x_{\perp})-\frac{\chi^{2}\mu^{2}L}{6}\,c_{T}\,x_{\perp}\,K_{0}(\mu\,x_{\perp})\,J_{1}(p_{\perp}x_{\perp})$
$\displaystyle\hskip
28.45274pt\times\left[1-3\mu\,x_{\perp}K_{1}(\mu\,x_{\perp})+\mu^{2}x_{\perp}^{2}K_{2}(\mu\,x_{\perp})\right]\cos{\theta}\Bigg{\\}}\,.$
(64)

For the particular case of the GW model, $\mathcal{V}^{\rm GW}(x_{\perp})$
asymptotically tends to a constant. Therefore, the integral in (IV) needs to
be regulated for values of $x_{\perp}\gg\frac{1}{\mu}$. In particular, one
should remove the non-scattering probability term associated with the
asymptotic behavior of $\mathcal{V}^{\rm GW}(x_{\perp})$, so that only genuine
broadening contributions are taken into account. Thus, following e.g.
Feal:2018jbm ; Barata:2020rdn , we consider a regularized broadening
distribution

$\displaystyle{\cal P}_{r}({\bm{p}})$ $\displaystyle={\cal
P}({\bm{p}})-2\pi\int\limits^{\infty}_{0}dx_{\perp}\,x_{\perp}\,e^{-\mathcal{V}^{\rm
GW}(\infty)L}\,J_{0}(p_{\perp}x_{\perp})\,,$ (65)

which differs from (IV) by a singular term (delta function) at
${\bm{p}}={\bm{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. Vitev:2008vk ; Xu:2014ica ;
Sievert:2019cwq ; Feal:2019xfl ; Barata:2020rdn ; Andres:2020vxs ;
Antiporda:2021hpk . 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 $\chi$, the medium length size
$L$, and the Debye mass $\mu$. Gradient effects are controlled by $\theta$ and
$c_{T}$. We take the medium size to be $L=5\,{\rm fm}$, the typical size of
nuclear matter, for all the plots. To explore the sensitivity of ${\cal
P}_{r}({\bm{p}})$ to the matter gradients, we chose to either take
$c_{T}=0.02$ or $c_{T}=0.002$ with the temperature gradient either parallel
($\theta=0$) or antiparallel ($\theta=\pi$) to $\bm{p}$. The particular values
of $c_{T}$ should be understood in the following way. In the hydrodynamic
regime the temperature gradient is expected not to be large, i.e.
$|\bm{\nabla}T|<T^{2}$, while the characteristic jet energy $E$ should be much
larger than the temperature $T$. If one takes
$\left|\frac{\bm{\nabla}T}{T^{2}}\right|\sim 1$ and
$\left|\frac{\bm{\nabla}T}{T^{2}}\right|\sim 0.1$, then the $c_{T}$ values
correspond to $E\sim 50\,T$, and one can further scale the value of the
gradient and jet energy together. Finally, when numerically evaluating (65),
we perform the remaining integration up to a scale of the order of
$\frac{1}{\mu}$. This ultraviolet cutoff was systematically varied in order to
ensure that the final results do not depend on its particular choice.

Figure 4: Numerical evaluation of (65) for six different medium
parametrizations. Each setup is characterized by its values for the opacity
$\chi$, the Debye mass $\mu$, and a medium length $L=5\,{\rm fm}$. The values
considered for the medium parameters are inspired by the ones typically found
in the literature Vitev:2008vk ; Xu:2014ica ; Sievert:2019cwq ; Feal:2019xfl ;
Barata:2020rdn ; Antiporda:2021hpk .

At the qualitative level, we observe that the gradient effects are small for
lower opacity, since the gradient corrections in (65) are proportional to
$\chi^{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
$\theta=0$), resulting in an increase of the broadening probability at higher
$p_{\perp}$. Conversely, in the case where ${\bm{p}}$ is antiparallel to
$\bm{\nabla}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_{\perp}$,
due to the fact that the inhomogeneous term always couples to an odd power of
${\bm{p}}$.

##  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. Wiedemann:2000ez ; Wiedemann:2000za ; Gyulassy:2002yv ;
Qiu:2003pm . Here we make the next step, showing that the agreement holds even
if the matter is not translationally-invariant 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 Mehtar-Tani:2012mfa ;
Blaizot:2013vha  would be modified in accordance with (44). Understanding the
resulting new factorized form is critical for some aspects of jet quenching
phenomenology Blaizot:2013vha ; Caucal:2019uvr ; Casalderrey-Solana:2014bpa .

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.
Antiporda:2021hpk . Particularly, it would be interesting to explore to what
extent the jet substructure grooming/tagging techniques, see e.g.
Larkoski:2014wba ; Chien:2016led ; Larkoski:2017jix ; Apolinario:2017qay ;
Caucal:2021cfb ; Karlberg:2021kwr  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
Ipp:2020mjc ; Carrington:2021dvw . 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. Agostini:2019avp ;
Hauksson:2021okc ; Altinoluk:2021lvu . 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. Armesto:2013fca ;
Wang:2001ifa .

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 jet-medium
interactions, see e.g. Casalderrey-Solana:2014bpa ; Casalderrey-Solana:2015vaa
; Casalderrey-Solana:2016jvj ; Rajagopal:2016uip ; Brewer:2017fqy ;
Brewer:2018mpk , one may hope to identify a set of new observables
distinguishing the two regimes. We leave this intriguing opportunity for
future work.

## Acknowledgements

The authors would like to thank C. Andres, N. Armesto, F. Dominguez, M.
Sievert, I. Vitev, and B. Wu for discussions and comments on this work. This
work is supported by European Research Council project ERC-2018-ADG-835105
YoctoLHC; by Maria de Maetzu excellence program under project MDM-2016-0692
and CEX2020-001035-M; by Spanish Research State Agency under project
PID2020-119632GB-I00; and by Xunta de Galicia (Centro singular de
investigación de Galicia accreditation 2019-2022), by European Union ERDF. The
work of A.S. is also supported by the Marie Sklodowska-Curie Individual
Fellowship under JetT project (project reference 101032858);

##  Appendix A The Gradient Corrections at the $N$th Order in Opacity

In this Appendix, we present the jet broadening distribution in an
inhomogeneous matter at leading orders in opacity, and discuss the general
structure of the opacity series. While $N=0$ case is trivial, the leading
gradient corrections at $N=1$ have been obtained in Sadofyev:2021ohn , and we
only quote the result here

$\displaystyle\left\langle\left|M\right|^{2}\right\rangle^{(1)}$
$\displaystyle=-\int\displaylimits_{0}^{L}\,dz\int\frac{d^{2}\bm{q}}{(2\pi)^{2}}\left\\{1+\left[\frac{(\bm{p}-\bm{q})}{E}z\right]\cdot\hat{\bm{g}}\right\\}\mathcal{V}_{1}(\bm{q})|J(E,\bm{p}-\bm{q})|^{2}\,,$
(66)

writing it in the notations of this work.

Turning to the $N=2$ case, we notice that there are five types of
contributions, which we have to consider separately. We start with the
simplest diagram involving no contact interactions, and the corresponding
amplitude reads

$\displaystyle iM_{11}$
$\displaystyle=\sum_{i_{1},i_{2}}\int\frac{d^{2}\bm{q}_{1}\,d^{2}\bm{q}_{2}}{(2\pi)^{4}}\,\left(it^{a}_{\mathrm{proj}}t_{i_{1}}^{a}\right)\left(it^{b}_{\mathrm{proj}}t_{i_{2}}^{b}\right)\,\theta_{2,1}\theta_{1,0}\,e^{-i\bm{q}_{1}\cdot\bm{x}_{i_{1}}}e^{-i\bm{q}_{2}\cdot\bm{x}_{i_{2}}}v_{i_{2}}(q_{2})v_{i_{1}}(q_{1})$
$\displaystyle\hskip 28.45274pt\times
e^{-i\frac{(p-q_{1}-q_{2})_{\perp}^{2}-p_{\perp}^{2}}{2E}z_{i_{1}}}e^{-i\frac{(p-q_{2})_{\perp}^{2}-p_{\perp}^{2}}{2E}\left(z_{i_{2}}-z_{i_{1}}\right)}\,J\left(E,\bm{p}-\bm{q}_{1}-\bm{q}_{2}\right)\,,$
(67)

where the subscript of $M$ indicates that there are two SB insertions. This
amplitude has the same form as in the well known homogeneous case, and we omit
its derivation. Squaring this contribution and averaging over the quantum
numbers, we find

$\displaystyle\left\langle\left|M\right|^{2}\right\rangle^{(2)}_{11,11}$
$\displaystyle=\mathcal{C}^{2}\int\limits^{L}_{0}\,dz_{2}\int\limits^{z_{2}}_{0}\,dz_{1}\,\int\frac{d^{2}\bm{q}_{1}\,d^{2}\bm{q}_{2}}{(2\pi)^{4}}\Bigg{\\{}1+\frac{\bm{p}-\bm{q}_{1}-\bm{q}_{2}}{E}z_{1}\cdot\hat{\bm{g}}_{1}+\frac{\bm{p}-\bm{q}_{1}-\bm{q}_{2}}{E}z_{1}\cdot\hat{\bm{g}}_{2}$
$\displaystyle+\frac{\bm{p}-\bm{q}_{2}}{E}(z_{2}-z_{1})\cdot\hat{\bm{g}}_{2}\Bigg{\\}}\,\rho_{2}\,\rho_{1}\,v_{2}^{2}(q_{2\perp}^{2})v_{1}^{2}(q_{1\perp}^{2})|J(E,\bm{p}-\bm{q}_{1}-\bm{q}_{2})|^{2}$
(68)

where the $\bm{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 $\frac{\partial}{\partial(q-\overline{q})_{\alpha}}$,
where $\overline{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

$\displaystyle iM_{12}$
$\displaystyle=\mathcal{C}\sum_{i_{1},i_{2}}\int\frac{d^{2}\bm{q}_{1}\,d^{2}\bm{q}_{2}\,d^{2}\bm{q}_{3}}{(2\pi)^{6}}\,\left(t^{a}_{\mathrm{proj}}t_{i_{1}}^{a}\right)\,\theta_{2,1}\theta_{1,0}\,e^{-i\bm{q}_{1}\cdot\bm{x}_{i_{1}}}e^{-i\left(\bm{q}_{2}+\bm{q}_{3}\right)\cdot\bm{x}_{i_{2}}}v_{i_{1}}(\bm{q}_{1})\mathcal{I}^{12}_{i_{2}}\left(\bm{q}_{2},\bm{q}_{3}\right)$
$\displaystyle\times
e^{-i\frac{(p-q_{1}-q_{2}-q_{3})_{\perp}^{2}-p_{\perp}^{2}}{2E}z_{i_{1}}}e^{-i\frac{(p-q_{2}-q_{3})_{\perp}^{2}-p_{\perp}^{2}}{2E}\left(z_{i_{2}}-z_{i_{1}}\right)}\,J\left(E,\bm{p}-\bm{q}_{1}-\bm{q}_{2}-\bm{q}_{3}\right)\,,$
(69)

where the subscript indicates that the first and second insertions are SB and
DB correspondingly, and

$\displaystyle\mathcal{I}^{12}_{i_{2}}$
$\displaystyle=2E\int\frac{dq_{3,z}}{(2\pi)}\frac{v_{i_{2}}\left(\bm{q}_{2},\frac{(p-q_{2}-q_{3})_{\perp}^{2}-p_{\perp}^{2}}{2E}-q_{3,z}\right)v_{i_{2}}\left(\bm{q}_{3},q_{3,z}\right)}{\left(p-q_{3}\right)^{2}+i\epsilon}$
$\displaystyle\hskip
28.45274pt\simeq-\frac{i}{2}v_{i_{2}}\left(q_{2\perp}^{2}\right)v_{i_{2}}\left(q_{3\perp}^{2}\right)\,.$
(70)

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

$\displaystyle\left\langle M_{12}M^{*}_{1}\right\rangle+c.c.$
$\displaystyle=-\mathcal{C}^{2}\int\limits^{L}_{0}\,dz_{2}\int\limits^{z_{2}}_{0}\,dz_{1}\,\int\frac{d^{2}\bm{q}_{1}\,d^{2}\bm{q}_{2}}{(2\pi)^{4}}\Bigg{\\{}1+\frac{\bm{p}-\bm{q}_{1}}{E}z_{1}\cdot\hat{\bm{g}}_{1}+\frac{\bm{p}-\bm{q}_{1}}{E}z_{1}\cdot\hat{\bm{g}}_{2}$
$\displaystyle+\frac{\bm{p}}{E}(z_{2}-z_{1})\cdot\hat{\bm{g}}_{2}\Bigg{\\}}\,\rho_{2}\,\rho_{1}\,v_{2}^{2}(q_{2\perp}^{2})v_{1}^{2}(q_{1\perp}^{2})|J(E,\bm{p}-\bm{q}_{1})|^{2}\,.$
(71)

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

$\displaystyle iM_{21}$
$\displaystyle=i\mathcal{C}\sum_{i_{1},i_{3}}\int\frac{d^{2}\bm{q}_{1}\,d^{2}\bm{q}_{2}\,d^{2}\bm{q}_{3}\,dq_{3,z}}{(2\pi)^{7}}\,\left(t^{a}_{\mathrm{proj}}t_{i_{3}}^{a}\right)\,\theta_{1,0}\,e^{-i(\bm{q}_{1}+\bm{q}_{2})\cdot\bm{x}_{i_{1}}}e^{-i\bm{q}_{3}\cdot\bm{x}_{i_{3}}}\mathcal{I}^{21}_{i_{1}}\left(\bm{q}_{1},\bm{q}_{2},\bm{q}_{3},q_{3,z}\right)$
$\displaystyle\times
e^{-i\frac{(p-q_{1}-q_{2}-q_{3})_{\perp}^{2}-p_{\perp}^{2}}{2E}z_{i_{1}}}\frac{2E\,v_{i_{3}}(q_{3})e^{-iq_{3,z}\left(z_{i_{3}}-z_{i_{1}}\right)}}{(p-q_{3})^{2}+i\epsilon}\,J\left(E,\bm{p}-\bm{q}_{1}-\bm{q}_{2}-\bm{q}_{3}\right)\,.$
(72)

The integral appearing due to the contact interaction is given by

$\displaystyle\mathcal{I}^{21}_{i_{1}}$
$\displaystyle=2E\int\frac{dq_{2,z}}{(2\pi)}\frac{v_{i_{1}}\left(\bm{q}_{1},\frac{(p-q_{1}-q_{2}-q_{3})_{\perp}^{2}-p_{\perp}^{2}}{2E}-q_{2,z}-q_{3,z}\right)v_{i_{1}}\left(\bm{q}_{2},q_{2,z}\right)}{\left(p-q_{2}-q_{3}\right)^{2}+i\epsilon}$
$\displaystyle\simeq-\frac{i}{2}v_{i_{1}}\left(q_{1\perp}^{2}\right)v_{i_{1}}\left(q_{2\perp}^{2}+q_{3,z}^{2}\right)+\frac{\tilde{\mu}_{1}+2\tilde{\mu}_{2}}{\tilde{\mu}_{1}\tilde{\mu}_{2}(\tilde{\mu}_{2}^{2}+q_{3,z}^{2})\left((\tilde{\mu}_{1}+\tilde{\mu}_{2})^{2}+q_{3,z}^{2}\right)}q_{3,z}\,,$
(73)

exhibiting a non-zero real part even at the leading order in the eikonal
expansion. Here $\tilde{\mu}_{n}^{2}\equiv\mu^{2}+\bm{q}_{n}^{2}$, and we keep
the subscript $i_{1}$ implicit in $\tilde{\mu}$. One can notice that
$\mathcal{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,

$\displaystyle iM_{21}$
$\displaystyle=\mathcal{C}\sum_{i_{1},i_{3}}\int\frac{d^{2}\bm{q}_{1}\,d^{2}\bm{q}_{2}\,d^{2}\bm{q}_{3}}{(2\pi)^{6}}\,\left(t^{a}_{\mathrm{proj}}t_{i_{3}}^{a}\right)\,\theta_{2,1}\theta_{1,0}\,e^{-i(\bm{q}_{1}+\bm{q}_{2})\cdot\bm{x}_{i_{1}}}e^{-i\bm{q}_{3}\cdot\bm{x}_{i_{3}}}\mathcal{I}^{12}_{i_{1}}\left(\bm{q}_{1},\bm{q}_{2}\right)$
$\displaystyle\hskip 28.45274pt\times
e^{-i\frac{(p-q_{1}-q_{2}-q_{3})_{\perp}^{2}-p_{\perp}^{2}}{2E}z_{i_{1}}}e^{-i\frac{(p-q_{3})_{\perp}^{2}-p_{\perp}^{2}}{2E}\left(z_{i_{3}}-z_{i_{1}}\right)}\,J\left(E,\bm{p}-\bm{q}_{1}-\bm{q}_{2}-\bm{q}_{3}\right)\,,$
(74)

where we have used that $\text{Re}\,\mathcal{I}^{21}$ is sub-eikonal after the
$q_{3,z}$ is substituted (the two integrals $\mathcal{I}^{12}$ and
$\mathcal{I}^{21}$ coincide in this limit). Turning to the contribution to the
squared amplitude, we find

$\displaystyle\left\langle M_{21}M^{*}_{1}\right\rangle+c.c.$
$\displaystyle=-\mathcal{C}^{2}\int\limits^{L}_{0}\,dz_{2}\int\limits^{z_{2}}_{0}\,dz_{1}\,\int\frac{d^{2}\bm{q}_{1}\,d^{2}\bm{q}_{2}}{(2\pi)^{4}}\Bigg{\\{}1+\frac{\bm{p}-\bm{q}_{2}}{E}z_{1}\cdot\hat{\bm{g}}_{1}+\frac{\bm{p}-\bm{q}_{2}}{E}z_{1}\cdot\hat{\bm{g}}_{2}$
$\displaystyle+\frac{\bm{p}-\bm{q}_{2}}{E}(z_{2}-z_{1})\cdot\hat{\bm{g}}_{2}\Bigg{\\}}\,\rho_{2}\,\rho_{1}\,v_{2}^{2}(q_{2\perp}^{2})v_{1}^{2}(q_{1\perp}^{2})|J(E,\bm{p}-\bm{q}_{2})|^{2}\,,$
(75)

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

$\displaystyle\left\langle iM_{22}\right\rangle$
$\displaystyle=\frac{1}{4}\mathcal{C}^{2}\sum_{i_{1},i_{3}}\int\frac{d^{2}\bm{q}_{1}\,d^{2}\bm{q}_{2}\,d^{2}\bm{q}_{3}\,d^{2}\bm{q}_{4}}{(2\pi)^{8}}\,\theta_{2,1}\theta_{1,0}\,e^{-i\left(\bm{q}_{1}+\bm{q}_{2}\right)\cdot\bm{x}_{i_{1}}}e^{-i\left(\bm{q}_{3}+\bm{q}_{4}\right)\cdot\bm{x}_{i_{3}}}$
$\displaystyle\hskip
14.22636pt\times\,e^{-i\frac{(p-q_{1}-q_{2}-q_{3}-q_{4})_{\perp}^{2}-p_{\perp}^{2}}{2E}z_{i_{1}}}e^{-i\frac{(p-q_{3}-q_{4})_{\perp}^{2}-p_{\perp}^{2}}{2E}\left(z_{i_{3}}-z_{i_{1}}\right)}$
$\displaystyle\hskip 28.45274pt\times
v_{i_{1}}(q_{1\perp}^{2})v_{i_{1}}(q_{2\perp}^{2})v_{i_{3}}(q_{3\perp}^{2})v_{i_{3}}(q_{4\perp}^{2})\,J\left(E,\bm{p}-\bm{q}_{1}-\bm{q}_{2}-\bm{q}_{3}-\bm{q}_{4}\right)\,,$
(76)

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

$\displaystyle\left\langle M_{22}M^{*}_{0}\right\rangle+c.c.$
$\displaystyle=\frac{1}{2}\mathcal{C}^{2}\int\limits^{L}_{0}\,dz_{2}\int\limits^{z_{2}}_{0}\,dz_{1}\,\int\frac{d^{2}\bm{q}_{1}\,d^{2}\bm{q}_{2}}{(2\pi)^{4}}\Bigg{\\{}1+\frac{\bm{p}}{E}z_{1}\cdot\hat{\bm{g}}_{1}+\frac{\bm{p}}{E}z_{1}\cdot\hat{\bm{g}}_{2}$
$\displaystyle\hskip
28.45274pt+\frac{\bm{p}}{E}(z_{2}-z_{1})\cdot\hat{\bm{g}}_{2}\Bigg{\\}}\,\rho_{2}\,\rho_{1}\,v_{2}^{2}(q_{2\perp}^{2})v_{1}^{2}(q_{1\perp}^{2})|J(E,\bm{p})|^{2}\,,$
(77)

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 Sadofyev:2021ohn  we write

$\displaystyle\left\langle M_{2}M^{*}_{2}\right\rangle$
$\displaystyle=\frac{1}{4}\mathcal{C}^{2}\int\limits^{L}_{0}\,dz_{2}\int\limits^{L}_{0}\,dz_{1}\,\int\frac{d^{2}\bm{q}_{1}\,d^{2}\bm{q}_{2}}{(2\pi)^{4}}\Bigg{\\{}1+\frac{\bm{p}}{E}z_{1}\cdot\hat{\bm{g}}_{1}+\frac{\bm{p}}{E}z_{2}\cdot\hat{\bm{g}}_{2}\Bigg{\\}}$
$\displaystyle\hskip
28.45274pt\times\rho_{2}\,\rho_{1}\,v_{2}^{2}(q_{2\perp}^{2})v_{1}^{2}(q_{1\perp}^{2})|J(E,\bm{p})|^{2}\,,$
(78)

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

$\displaystyle\left\langle\left|M\right|^{2}\right\rangle^{(2)}$
$\displaystyle=\int\displaylimits_{0}^{L}\,dz_{2}\int\displaylimits_{0}^{z_{2}}\,dz_{1}\int\frac{d^{2}\bm{q}}{(2\pi)^{2}}\Bigg{\\{}1+\frac{\bm{p}-\bm{q}_{1}-\bm{q}_{2}}{E}z_{1}\cdot\hat{\bm{g}}_{1}+\frac{\bm{p}-\bm{q}_{1}-\bm{q}_{2}}{E}z_{1}\cdot\hat{\bm{g}}_{2}$
$\displaystyle+\frac{\bm{p}-\bm{q}_{2}}{E}(z_{2}-z_{1})\cdot\hat{\bm{g}}_{2}\Bigg{\\}}\mathcal{V}_{1}(\bm{q}_{1})\mathcal{V}_{2}(\bm{q}_{2})|J(E,\bm{p}-\bm{q}_{1}-\bm{q}_{2})|^{2}\,,$
(79)

in agreement with the general expression (II). Repeating the same calculation
for the higher orders in opacity, we find that the structure of (II) is
reproduced.

To understand the origin of the all order pattern, one may consider the
different classes of the contributions to the amplitude squared, distinguished
by the color pairing. For the given pairing, one can consider three higher
order diagrams obtained by attaching two extra in-medium field insertions in
all possible ways – one extra SB interaction on each of the sides of the cut
and one extra DB interaction on either one of the sides of the cut. Thus, for
each color pairing at $N=n$ we obtain a set of contributions at $N=n+1$.
Comparing to the $N=2$ case, we can see that each of these contributions has
the same structure as in the case of the already studied diagrams. Paying a
particular attention to the LPM phases, we can see that at the first order in
gradients, all the phases may depend only on the momenta of the two additional
interactions which are responsible for the gradient contributions arising at
$N=n+1$.

##  Appendix B Shift Operator and $\delta$-function

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
$\delta$-functions in the last line of (III). Since the $\delta$-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

$\displaystyle\int dx\,e^{-ipx}e^{ax}\equiv(2\pi)\,\delta(p+ia)\,.$ (80)

For a real $a$, one can understand it as a momentum representation of the
shift operator $e^{ax}=e^{ia\frac{\partial}{\partial p}}$ with a purely
imaginary argument $ia$. Thus, we can write

$\displaystyle\int dx\,e^{-ipx}\,e^{ax}=(2\pi)e^{ia\frac{\partial}{\partial
p}}\delta(p)=(2\pi)\,\delta(p+ia)\,\,.$ (81)

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$-delta. For instance,
the one-dimensional analogue to the integral of interest in the main text
reduces to

$\displaystyle\int
dp_{in}d\bar{p}_{in}\,\delta(\bar{p}_{in}-p_{in}-ia)J(p_{in})J^{\dagger}(\bar{p}_{in})$
$\displaystyle=\int dp_{in}\,|J(p_{in})|^{2}+ia\int
dp_{in}d\bar{p}_{in}\,(2\pi)\,\delta(\bar{p}_{in}-p_{in})\frac{\partial}{\partial(\bar{p}_{in}-p_{in})}\left[J(p_{in})J^{\dagger}(\bar{p}_{in})\right]$
$\displaystyle=\int dp_{in}\,|J(p_{in})|^{2}+\frac{ia}{2}\int
dqdQ\,(2\pi)\,\delta(q)\frac{\partial}{\partial
q}\left[J\left(\frac{Q-q}{2}\right)J^{\dagger}\left(\frac{Q+q}{2}\right)\right]$
$\displaystyle=\int dp_{in}\,|J(p_{in})|^{2}\,.$ (82)

## References

  * (1) J. D. Bjorken, “Energy Loss of Energetic Partons in Quark - Gluon Plasma: Possible Extinction of High p(t) Jets in Hadron - Hadron Collisions”.
  * (2) Y. Mehtar-Tani, J. G. Milhano & K. Tywoniuk, “Jet physics in heavy-ion collisions”, Int. J. Mod. Phys. A 28, 1340013 (2013), arXiv:1302.2579.
  * (3) G.-Y. Qin & X.-N. Wang, “Jet quenching in high-energy heavy-ion collisions”, Int. J. Mod. Phys. E 24, 1530014 (2015), arXiv:1511.00790.
  * (4) J.-P. Blaizot & Y. Mehtar-Tani, “Jet Structure in Heavy Ion Collisions”, Int. J. Mod. Phys. E 24, 1530012 (2015), arXiv:1503.05958.
  * (5) M. D. Sievert & I. Vitev, “Quark branching in QCD matter to any order in opacity beyond the soft gluon emission limit”, Phys. Rev. D98, 094010 (2018), arXiv:1807.03799.
  * (6) I. Vitev & M. Gyulassy, “High $p_{T}$ tomography of $d$ \+ Au and Au+Au at SPS, RHIC, and LHC”, Phys. Rev. Lett. 89, 252301 (2002), hep-ph/0209161.
  * (7) JET Collaboration, K. M. Burke et al., “Extracting the jet transport coefficient from jet quenching in high-energy heavy-ion collisions”, Phys. Rev. C 90, 014909 (2014), arXiv:1312.5003.
  * (8) B. Betz & M. Gyulassy, “Constraints on the Path-Length Dependence of Jet Quenching in Nuclear Collisions at RHIC and LHC”, JHEP 1408, 090 (2014), arXiv:1404.6378, [Erratum: JHEP 10, 043 (2014)].
  * (9) L. Apolinário, J. G. Milhano, G. P. Salam & C. A. Salgado, “Probing the time structure of the quark-gluon plasma with top quarks”, Phys. Rev. Lett. 120, 232301 (2018), arXiv:1711.03105.
  * (10) H. T. Li, Z. L. Liu & I. Vitev, “Heavy meson tomography of cold nuclear matter at the electron-ion collider”, Phys. Lett. B 816, 136261 (2021), arXiv:2007.10994.
  * (11) H. T. Li & I. Vitev, “Nuclear Matter Effects on Jet Production at Electron-Ion Colliders”, Phys. Rev. Lett. 126, 252001 (2021), arXiv:2010.05912.
  * (12) M. Arratia, Z.-B. Kang, A. Prokudin & F. Ringer, “Jet-based measurements of Sivers and Collins asymmetries at the future electron-ion collider”, Phys. Rev. D 102, 074015 (2020), arXiv:2007.07281.
  * (13) Y. He, L.-G. Pang & X.-N. Wang, “Gradient Tomography of Jet Quenching in Heavy-Ion Collisions”, Phys. Rev. Lett. 125, 122301 (2020), arXiv:2001.08273.
  * (14) L. Apolinário, A. Cordeiro & K. Zapp, “Time reclustering for jet quenching studies”, Eur. Phys. J. C 81, 561 (2021), arXiv:2012.02199.
  * (15) A. V. Sadofyev, M. D. Sievert & I. Vitev, “Ab initio coupling of jets to collective flow in the opacity expansion approach”, Phys. Rev. D 104, 094044 (2021), arXiv:2104.09513.
  * (16) Y.-L. Du, D. Pablos & K. Tywoniuk, “Jet tomography in heavy ion collisions with deep learning”, arXiv:2106.11271.
  * (17) L. Antiporda, J. Bahder, H. Rahman & M. D. Sievert, “Jet Drift and Collective Flow in Heavy-Ion Collisions”, arXiv:2110.03590.
  * (18) M. Gyulassy & X.-n. Wang, “Multiple collisions and induced gluon Bremsstrahlung in QCD”, Nucl. Phys. B420, 583 (1994), nucl-th/9306003.
  * (19) R. Baier, Y. L. Dokshitzer, A. H. Mueller, S. Peigne & D. Schiff, “Radiative energy loss and p(T) broadening of high-energy partons in nuclei”, Nucl. Phys. B484, 265 (1997), hep-ph/9608322.
  * (20) B. G. Zakharov, “Fully quantum treatment of the Landau-Pomeranchuk-Migdal effect in QED and QCD”, JETP Lett. 63, 952 (1996), hep-ph/9607440.
  * (21) R. Baier, Y. L. Dokshitzer, A. H. Mueller, S. Peigne & D. Schiff, “Radiative energy loss of high-energy quarks and gluons in a finite volume quark - gluon plasma”, Nucl. Phys. B 483, 291 (1997), hep-ph/9607355.
  * (22) B. G. Zakharov, “Radiative energy loss of high-energy quarks in finite size nuclear matter and quark - gluon plasma”, JETP Lett. 65, 615 (1997), hep-ph/9704255.
  * (23) R. Baier, Y. L. Dokshitzer, A. H. Mueller & D. Schiff, “Medium induced radiative energy loss: Equivalence between the BDMPS and Zakharov formalisms”, Nucl. Phys. B 531, 403 (1998), hep-ph/9804212.
  * (24) M. Gyulassy, P. Levai & I. Vitev, “Jet quenching in thin quark gluon plasmas. 1. Formalism”, Nucl. Phys. B571, 197 (2000), hep-ph/9907461.
  * (25) M. Gyulassy, P. Levai & I. Vitev, “NonAbelian energy loss at finite opacity”, Phys. Rev. Lett. 85, 5535 (2000), nucl-th/0005032.
  * (26) M. Gyulassy, P. Levai & I. Vitev, “Reaction operator approach to nonAbelian energy loss”, Nucl. Phys. B594, 371 (2001), nucl-th/0006010.
  * (27) M. Gyulassy, P. Levai & I. Vitev, “Reaction operator approach to multiple elastic scatterings”, Phys. Rev. D 66, 014005 (2002), nucl-th/0201078.
  * (28) X.-N. Wang & X.-f. Guo, “Multiple parton scattering in nuclei: Parton energy loss”, Nucl. Phys. A696, 788 (2001), hep-ph/0102230.
  * (29) B.-W. Zhang, E. Wang & X.-N. Wang, “Heavy quark energy loss in nuclear medium”, Phys. Rev. Lett. 93, 072301 (2004), nucl-th/0309040.
  * (30) P. B. Arnold, G. D. Moore & L. G. Yaffe, “Photon and gluon emission in relativistic plasmas”, JHEP 0206, 030 (2002), hep-ph/0204343.
  * (31) M. Djordjevic & M. Gyulassy, “Heavy quark radiative energy loss in QCD matter”, Nucl. Phys. A733, 265 (2004), nucl-th/0310076.
  * (32) Y. Mehtar-Tani, “Relating the description of gluon production in pA collisions and parton energy loss in AA collisions”, Phys. Rev. C 75, 034908 (2007), hep-ph/0606236.
  * (33) M. Gyulassy, I. Vitev & X. Wang, “High p(T) azimuthal asymmetry in noncentral A+A at RHIC”, Phys. Rev. Lett. 86, 2537 (2001), nucl-th/0012092.
  * (34) M. Gyulassy, I. Vitev, X.-N. Wang & P. Huovinen, “Transverse expansion and high p(T) azimuthal asymmetry at RHIC”, Phys. Lett. B 526, 301 (2002), nucl-th/0109063.
  * (35) R. Baier, Y. L. Dokshitzer, A. H. Mueller & D. Schiff, “Radiative energy loss of high-energy partons traversing an expanding QCD plasma”, Phys. Rev. C 58, 1706 (1998), hep-ph/9803473.
  * (36) N. Armesto, C. A. Salgado & U. A. Wiedemann, “Measuring the collective flow with jets”, Phys. Rev. Lett. 93, 242301 (2004), hep-ph/0405301.
  * (37) N. Armesto, C. A. Salgado & U. A. Wiedemann, “Low-p(T) collective flow induces high-p(T) jet quenching”, Phys. Rev. C 72, 064910 (2005), hep-ph/0411341.
  * (38) R. Baier, A. H. Mueller & D. Schiff, “How does transverse (hydrodynamic) flow affect jet-broadening and jet-quenching ?”, Phys. Lett. B649, 147 (2007), nucl-th/0612068.
  * (39) H. Liu, K. Rajagopal & U. A. Wiedemann, “Wilson loops in heavy ion collisions and their calculation in AdS/CFT”, JHEP 0703, 066 (2007), hep-ph/0612168.
  * (40) T. Renk, J. Ruppert, C. Nonaka & S. A. Bass, “Jet-quenching in a 3D hydrodynamic medium”, Phys. Rev. C75, 031902 (2007), nucl-th/0611027.
  * (41) M. Lekaveckas & K. Rajagopal, “Effects of Fluid Velocity Gradients on Heavy Quark Energy Loss”, JHEP 1402, 068 (2014), arXiv:1311.5577.
  * (42) K. Rajagopal & A. V. Sadofyev, “Chiral drag force”, JHEP 1510, 018 (2015), arXiv:1505.07379.
  * (43) A. V. Sadofyev & Y. Yin, “The charmonium dissociation in an “anomalous wind””, JHEP 1601, 052 (2016), arXiv:1510.06760.
  * (44) S. Li, K. A. Mamo & H.-U. Yee, “Jet quenching parameter of the quark-gluon plasma in a strong magnetic field: Perturbative QCD and AdS/CFT correspondence”, Phys. Rev. D94, 085016 (2016), arXiv:1605.00188.
  * (45) J. Reiten & A. V. Sadofyev, “Drag force to all orders in gradients”, JHEP 2007, 146 (2020), arXiv:1912.08816.
  * (46) I. Y. Aref’eva, A. A. Golubtsova & E. Gourgoulhon, “Holographic drag force in 5d Kerr-AdS black hole”, arXiv:2004.12984.
  * (47) P. Caucal, E. Iancu & G. Soyez, “Deciphering the $z_{g}$ distribution in ultrarelativistic heavy ion collisions”, JHEP 1910, 273 (2019), arXiv:1907.04866.
  * (48) J. Casalderrey-Solana, D. C. Gulhan, J. G. Milhano, D. Pablos & K. Rajagopal, “A Hybrid Strong/Weak Coupling Approach to Jet Quenching”, JHEP 1410, 019 (2014), arXiv:1405.3864, [Erratum: JHEP 09, 175 (2015)].
  * (49) M. D. Sievert, I. Vitev & B. Yoon, “A complete set of in-medium splitting functions to any order in opacity”, Phys. Lett. B 795, 502 (2019), arXiv:1903.06170.
  * (50) J. H. Putschke et al., “The JETSCAPE framework”, arXiv:1903.07706.
  * (51) Y. He, T. Luo, X.-N. Wang & Y. Zhu, “Linear Boltzmann Transport for Jet Propagation in the Quark-Gluon Plasma: Elastic Processes and Medium Recoil”, Phys. Rev. C 91, 054908 (2015), arXiv:1503.03313, [Erratum: Phys.Rev.C 97, 019902 (2018)].
  * (52) L. D. McLerran & R. Venugopalan, “Gluon distribution functions for very large nuclei at small transverse momentum”, Phys. Rev. D 49, 3352 (1994), hep-ph/9311205.
  * (53) L. D. McLerran & R. Venugopalan, “Computing quark and gluon distribution functions for very large nuclei”, Phys. Rev. D 49, 2233 (1994), hep-ph/9309289.
  * (54) H. Kleinert, “Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets”.
  * (55) J.-P. Blaizot, F. Dominguez, E. Iancu & Y. Mehtar-Tani, “Medium-induced gluon branching”, JHEP 1301, 143 (2013), arXiv:1209.4585.
  * (56) L. Apolinário, N. Armesto, J. G. Milhano & C. A. Salgado, “Medium-induced gluon radiation and colour decoherence beyond the soft approximation”, JHEP 1502, 119 (2015), arXiv:1407.0599.
  * (57) A. Majumder, B. Muller & S. A. Bass, “Longitudinal Broadening of Quenched Jets in Turbulent Color Fields”, Phys. Rev. Lett. 99, 042301 (2007), hep-ph/0611135.
  * (58) J.-P. Blaizot, F. Dominguez, E. Iancu & Y. Mehtar-Tani, “Probabilistic picture for medium-induced jet evolution”, JHEP 1406, 075 (2014), arXiv:1311.5823.
  * (59) K. Kutak, W. Placzek & R. Straka, “Solutions of evolution equations for medium-induced QCD cascades”, Eur. Phys. J. C 79, 317 (2019), arXiv:1811.06390.
  * (60) J. Barata, F. Domίnguez, C. A. Salgado & V. Vila, “A modified in-medium evolution equation with color coherence”, JHEP 2105, 148 (2021), arXiv:2101.12135.
  * (61) I. Soudi, “Energy loss and equilibration of a highly energetic parton in QCD plasmas”.
  * (62) J. Barata & Y. Mehtar-Tani, “Improved opacity expansion at NNLO for medium induced gluon radiation”, JHEP 2010, 176 (2020), arXiv:2004.02323.
  * (63) J. Barata, Y. Mehtar-Tani, A. Soto-Ontoso & K. Tywoniuk, “Revisiting transverse momentum broadening in dense QCD media”, Phys. Rev. D 104, 054047 (2021), arXiv:2009.13667.
  * (64) M. Djordjevic, “Theoretical formalism of radiative jet energy loss in a finite size dynamical QCD medium”, Phys. Rev. C 80, 064909 (2009), arXiv:0903.4591.
  * (65) J. Barata, Y. Mehtar-Tani, A. Soto-Ontoso & K. Tywoniuk, “Medium-induced radiative kernel with the Improved Opacity Expansion”, arXiv:2106.07402.
  * (66) N. Armesto, C. A. Salgado & U. A. Wiedemann, “Medium induced gluon radiation off massive quarks fills the dead cone”, Phys. Rev. D 69, 114003 (2004), hep-ph/0312106.
  * (67) Y. Mehtar-Tani & K. Tywoniuk, “Improved opacity expansion for medium-induced parton splitting”, JHEP 2006, 187 (2020), arXiv:1910.02032.
  * (68) Y. Mehtar-Tani, “Gluon bremsstrahlung in finite media beyond multiple soft scattering approximation”, JHEP 1907, 057 (2019), arXiv:1903.00506.
  * (69) X. G. Feal, “High energy scattering and emission in QED&QCD media”, arXiv:1812.06903.
  * (70) I. Vitev & B.-W. Zhang, “A Systematic study of direct photon production in heavy ion collisions”, Phys. Lett. B 669, 337 (2008), arXiv:0804.3805.
  * (71) J. Xu, A. Buzzatti & M. Gyulassy, “Azimuthal jet flavor tomography with CUJET2.0 of nuclear collisions at RHIC and LHC”, JHEP 1408, 063 (2014), arXiv:1402.2956.
  * (72) X. Feal, C. A. Salgado & R. A. Vazquez, “Jet quenching test of the QCD matter created at RHIC and the LHC needs opacity-resummed medium induced radiation”, Phys. Lett. B 816, 136251 (2021), arXiv:1911.01309.
  * (73) C. Andres, L. Apolinário & F. Dominguez, “Medium-induced gluon radiation with full resummation of multiple scatterings for realistic parton-medium interactions”, JHEP 2007, 114 (2020), arXiv:2002.01517.
  * (74) U. A. Wiedemann, “Transverse dynamics of hard partons in nuclear media and the QCD dipole”, Nucl. Phys. B 582, 409 (2000), hep-ph/0003021.
  * (75) U. A. Wiedemann, “Gluon radiation off hard quarks in a nuclear environment: Opacity expansion”, Nucl. Phys. B 588, 303 (2000), hep-ph/0005129.
  * (76) J.-w. Qiu & I. Vitev, “Transverse momentum diffusion and broadening of the back-to-back dihadron correlation function”, Phys. Lett. B 570, 161 (2003), nucl-th/0306039.
  * (77) Y. Mehtar-Tani, C. A. Salgado & K. Tywoniuk, “The Radiation pattern of a QCD antenna in a dense medium”, JHEP 1210, 197 (2012), arXiv:1205.5739.
  * (78) A. J. Larkoski, S. Marzani, G. Soyez & J. Thaler, “Soft Drop”, JHEP 1405, 146 (2014), arXiv:1402.2657.
  * (79) Y.-T. Chien & I. Vitev, “Probing the Hardest Branching within Jets in Heavy-Ion Collisions”, Phys. Rev. Lett. 119, 112301 (2017), arXiv:1608.07283.
  * (80) A. J. Larkoski, I. Moult & B. Nachman, “Jet Substructure at the Large Hadron Collider: A Review of Recent Advances in Theory and Machine Learning”, Phys. Rept. 841, 1 (2020), arXiv:1709.04464.
  * (81) L. Apolinário, J. G. Milhano, M. Ploskon & X. Zhang, “Novel subjet observables for jet quenching in heavy-ion collisions”, Eur. Phys. J. C 78, 529 (2018), arXiv:1710.07607.
  * (82) P. Caucal, A. Soto-Ontoso & A. Takacs, “Dynamically groomed jet radius in heavy-ion collisions”, arXiv:2111.14768.
  * (83) A. Karlberg, G. P. Salam, L. Scyboz & R. Verheyen, “Spin correlations in final-state parton showers and jet observables”, Eur. Phys. J. C 81, 681 (2021), arXiv:2103.16526.
  * (84) A. Ipp, D. I. Müller & D. Schuh, “Anisotropic momentum broadening in the 2+1D Glasma: analytic weak field approximation and lattice simulations”, Phys. Rev. D 102, 074001 (2020), arXiv:2001.10001.
  * (85) M. E. Carrington, A. Czajka & S. Mrowczynski, “Jet quenching in glasma”, arXiv:2112.06812.
  * (86) P. Agostini, T. Altinoluk & N. Armesto, “Non-eikonal corrections to multi-particle production in the Color Glass Condensate”, Eur. Phys. J. C 79, 600 (2019), arXiv:1902.04483.
  * (87) S. Hauksson, S. Jeon & C. Gale, “The momentum broadening of energetic partons in an anisotropic plasma”, arXiv:2109.04575.
  * (88) T. Altinoluk & G. Beuf, “Quark and scalar propagators at next-to-eikonal accuracy in the CGC through a dynamical background gluon field”, arXiv:2109.01620.
  * (89) N. Armesto, H. Ma, M. Martinez, Y. Mehtar-Tani & C. A. Salgado, “Coherence Phenomena between Initial and Final State Radiation in a Dense QCD Medium”, JHEP 1312, 052 (2013), arXiv:1308.2186.
  * (90) J. Casalderrey-Solana, D. C. Gulhan, J. G. Milhano, D. Pablos & K. Rajagopal, “Predictions for Boson-Jet Observables and Fragmentation Function Ratios from a Hybrid Strong/Weak Coupling Model for Jet Quenching”, JHEP 1603, 053 (2016), arXiv:1508.00815.
  * (91) J. Casalderrey-Solana, D. Gulhan, G. Milhano, D. Pablos & K. Rajagopal, “Angular Structure of Jet Quenching Within a Hybrid Strong/Weak Coupling Model”, JHEP 1703, 135 (2017), arXiv:1609.05842.
  * (92) K. Rajagopal, A. V. Sadofyev & W. van der Schee, “Evolution of the jet opening angle distribution in holographic plasma”, Phys. Rev. Lett. 116, 211603 (2016), arXiv:1602.04187.
  * (93) J. Brewer, K. Rajagopal, A. Sadofyev & W. Van Der Schee, “Evolution of the Mean Jet Shape and Dijet Asymmetry Distribution of an Ensemble of Holographic Jets in Strongly Coupled Plasma”, JHEP 1802, 015 (2018), arXiv:1710.03237.
  * (94) J. Brewer, A. Sadofyev & W. van der Schee, “Jet shape modifications in holographic dijet systems”, arXiv:1809.10695.