aainstitutetext: Physics Department, Brookhaven National Laboratory, Upton, NY
11973, USAbbinstitutetext: Instituto Galego de Física de Altas Enerxías IGFAE,
Universidade de Santiago de Compostela, E- 15782 Galicia-Spainccinstitutetext:
Axencia Galega de Innovación (GAIN), Xunta de Galicia, Galicia-
Spainddinstitutetext: Laboratório de Instrumentação e Física Experimental de
Partículas (LIP), Av. Prof. Gama Pinto, 2, 1649-003 Lisbon,
Portugaleeinstitutetext: Departamento de Física, Instituto Superior Técnico
(IST), Universidade de Lisboa, Av. Rovisco Pais 1, 1049-001 Lisbon, Portugal

# Gluon to $q\bar{q}$ antenna in anisotropic QCD matter: spin-polarized and
azimuthal jet observables

João Barata b,c Carlos A. Salgado b,d,e João M. Silva
jlourenco@bnl.gov,joao.m.da.silva<EMAIL_ADDRESS>
###### Abstract

We study the production of a quark-antiquark antenna in the presence of a
dense and anisotropic QCD medium. We assume the antenna to originate from an
unpolarized gluon state, and consider both massless and massive final states.
The medium anisotropy is captured by allowing the jet quenching coefficient to
take different magnitudes in orthogonal directions with respect to the jet
axis. We find that the final particle distribution is sensitive to the medium
anisotropy, and more importantly, that this effect couples directly to the
helicity/spin of the final states. We propose to look into these effects by
performing a Fourier decomposition of the particle distribution inside the
jet. In our medium model, we find that the spin independent terms contribute
to the even harmonics of the cosine series. The helicity/spin dependence
enters only through the sine Fourier series. We further explore the spin
dependence by examining the degree of polarization of the final states in
different directions. Our results indicate that the anisotropies present in
the QCD matter produced in heavy ion collisions can be probed by studying
azimuthal and spin observables inside jets.

##  1 Introduction

The quark-gluon plasma (QGP) is a hot and dense QCD state of matter produced
in the early epochs of the Universe, where quarks and gluons are not confined
inside hadrons. On Earth, the QGP has been experimentally produced in the
aftermath of ultrarelativistic heavy ion collisions (URHICs), at the
Relativistic Heavy Ion Collider (RHIC)PHENIX:2001hpc ; STAR:2005gfr ;
STAR:2020xiv  and the Large Hadron Collider (LHC) Muller:2012zq ;
ATLAS:2018gwx ; CMS:2021vui . These experiments have allowed to test the many-
body properties of QCD in the most extreme regimes, and lead to the remarkable
discovery that the QGP is a nearly ideal liquid, see e.g. Busza:2018rrf  for a
review on URHICs.

The QGP produced at RHIC and the LHC appears as an intermediate phase of a
complex evolution, bridging a sequence of different QCD states of matter.
Generally, it is believed that in the first $1$ fm/c (or $3$ ys) after the
ions collide, the produced matter is in a highly anisotropic and out-of-
equilibrium state, dominated by strong polarized chromoelectric and magnetic
fields, which eventually evolve into an equilibrium state preceding the QGP
phase Schlichting:2019abc ; Berges:2020fwq . Afterwards, the QGP expands
hydrodynamicaly, with spatial gradients driving the evolution, until the
plasma temperature is sufficiently low, at which point the confinement
mechanism leads to the production of a gas of hadrons, later detected by the
experiment Jaiswal:2016hex ; Busza:2018rrf .111Despite the existence of
several sophisticated models describing the many parts of this evolution,
there is no overall widely accepted theoretical picture, nor are the current
measurements capable of probing all the aspects characterizing the different
stages of URHICs uniquely. To probe the QCD matter produced in URHICs, and
ultimately extract the properties of the QGP, it is common to use energetic
probes generated in the same event, which are sensitive to the underlying
medium as they evolve through it. One of the most popular and powerful ones
are hadronic jets, i.e. collimated QCD cascades originating from highly
boosted partons, whose fragmentation occurs simultaneously with the matter
evolution. The presence of a background medium induces modifications to the
final particle distribution inside jets, with part of the energy thermalizing
to the medium. The study of these medium induced effects on jets is typically
referred to as jet quenching Qin:2015srf ; Cunqueiro:2021wls ;
Apolinario:2022vzg . Since jets’ evolution covers a wide range of energy
scales, they offer an ideal probe of the spacetime structure of the QGP and
its transport properties.

So far, most theoretical jet quenching studies have focused on understanding
how mechanisms such as parton momentum broadening or stimulated radiative
energy loss affect jet observables and can be related to fundamental
properties of the plasma. Such studies are typically performed using
simplified medium models, mainly focused on capturing the QGP phase, commonly
assuming the medium to be static, isotropic, and homogeneous. However, as
mentioned above, the states of matter produced in real URHICs can be highly
dynamical and anisotropic, either due to the presence of polarized fields or
due to spatial gradients, for example. If such type of effects are
sufficiently strong, one should expect the final particle distribution inside
jets to reflect an angular modulation, which can be correlated to the
structure of the medium.222We note nonetheless that since the hydrodynamic
expansion tends to isotropize the system, many of these effects can be
ultimately washed out. The observation of such an effect is highly non-trivial
from the experimental point of view, since it requires correlating the jet
measurement to soft sector physics, but it would allow to use jets as detailed
probes of the structure of matter in URHICs. From the theory side, including
such contributions into account is technically challenging, and requires
further expanding the current theoretical toolkit Wang:2002ri ; Vitev:2004bh ;
Armesto:2004pt .

In this work, we take a step in this direction, by studying how a quark-
antiquark state branching from an initial unpolarized gluon couples to the
anisotropies of an underlying QCD medium, similarly to what was done recently
for pure gluonic systems Hauksson:2023tze . Assuming the medium to be
sufficiently long, we compute the leading order cross-sections for this
process in the large $N_{c}$ limit, and show that the final distributions are
sensitive to the matter structure. In particular, we find that helicity/spin
states couple distinctly to the medium anisotropies, and thus they offer a way
to extract the details of the underlying matter, not accessible with standard
jet observables. To this end, we propose extracting the azimuthal harmonic
inside jets, similar to what is done in the soft sector of URHICs, see also
Baty:2021ugw ; Zhao:2024wqs . We show below that, for the case where the jet
transport coefficient is distinct in two orthogonal directions transverse to
the jet axis, the leading terms in the harmonic decomposition of the particle
distribution inside a jet, for small anisotropy, read

$\displaystyle\frac{2\pi}{dN^{h}/dz}\frac{dN^{h}}{dzd\phi}=1+v_{2}\cos\left(2\phi\right)+w_{2}^{(h)}\sin\left(2\phi\right)\,,$
(1)

where $dN^{h}$ denotes the helicity/spin dependent final particle
distribution, $\phi$ is the azimuthal angle around the jet axis, $z$ is the
energy sharing fraction of the antenna, $v_{2},w_{2}$ vanish in isotropic
matter, and all the helicity/spin dependence lies in $w^{(h)}_{2}$.

This work is structured as follows. In Section 2, we provide the $g\to
q\bar{q}$ cross-section in the presence of an anisotropic medium, for both
massless and massive final states. We discuss jet observables capable of
probing and extracting the features of the medium anisotropy in Section 3,
followed by a summary and discussion of our results in Section 4. Several
technical details of the calculation are shown in Appendices A and B.

##  2 Quark pair production in a dense anisotropic QCD medium

We start by deriving the double differential cross-section for the process
presented in Fig. 1. This consists in the branching of a gluon to a $q\bar{q}$
pair (i.e. a quark antenna) inside a medium of length $L^{+}$ and density
$n(x^{+})$, where $x^{+}$ denotes the light-cone time. We shall assume that
the medium length is sufficiently long, such that we only take into
consideration processes where the branching takes place inside the medium.
Nonetheless, the calculation shown below can be easily extended to
fragmentation in finite media. We compute the quark antenna cross-section for
both massless and massive final states.

The most important element in our calculation is that we assume the medium to
be anisotropic. In particular, following Hauksson:2023tze , we incorporate
this picture by allowing the jet quenching parameter $\hat{q}$ to take
different values along two orthogonal directions in the transverse plane to
the jet axis, which we shall denote as $x$ and $y$. The case when
$\hat{q}_{x}=\hat{q}_{y}$ corresponds to the limit of an isotropic background.
As we show below, even when considering an unpolarized initial state, the
final particle distribution has a non-trivial azimuthal profile which depends
on the final states’ polarization. This suggests that studying the
helicity/spin dependent final state distributions inside jets can give non-
trivial information regarding the structure of the underlying medium; we
further explore this idea in the next section. Nonetheless, we note that our
matter model is rather simplistic, and many effects that tend to isotropize
the system, could qualitative change the results. We will disregard such
contributions in this work. Finally, in what follows we shall not discuss the
physics from which the anisotropic $\hat{q}$ emerges. We note that equivalent
forms for this parameter have been found in the early stages of heavy ions
collisions Boguslavski:2023waw ; Boguslavski:2023alu ; Ipp:2020mjc ;
Ipp:2020nfu ; Avramescu:2023qvv ; Carrington:2021dvw ; Carrington:2022bnv ;
Barata:2024xwy , and in the QGP epoch Kuzmin:2023hko ; Barata:2023zqg ;
Barata:2023qds ; Barata:2022utc ; Andres:2022ndd ; Barata:2022krd ;
Sadofyev:2021ohn ; He:2020iow ; Xiao:2024ffk ; Fu:2022idl ; He:2022evt ;
Kuzmin:2024smy .333See also Hauksson:2020wsm ; Hauksson:2021okc  for related
work. How all these different sources give rise to an anisotropic transport
coefficient and its impact on the evolution of the partonic cascade in the
medium are active topics of research. The observables proposed in this work
would contribute to pin down the relevant dynamics associated to an
anisotropic transport coefficient.

Figure 1: Leading order diagrams contributing to the antenna production. We
indicate the most relevant spacetime, polarization and color indices entering
the main text calculation, representing the (conjugate) amplitude on top
(bottom). The production point and initial momentum of the gluon is assumed to
be the same on both sides of the cut. Notice that although the diagram is the
same as in vacuum propagation, here one should associate a medium modified
propagator to each partonic leg.

We consider first the case of massless quark states. Since we are interested
in the limit of a dense QCD background medium, when computing the diagrams
shown in Fig. 1, one has to resum contributions describing the scattering of
the jet partons with the medium constituents. This can be efficiently
performed in a path integral formalism, as originally discussed by BDMPS-Z
Zakharov:1996fv ; Baier:1996kr , see also Attems:2022ubu  for a recent related
calculation. Very briefly,444See e.g. Mehtar-Tani:2013pia ; Blaizot:2015lma
for recent reviews. in this approach, to each line in Fig. 1 one associates an
effective quantum mechanical propagator, which for massless quarks takes the
form

$\mathcal{G}^{ij}(y^{+},\boldsymbol{y};x^{+},\boldsymbol{x}|p^{+})=\int^{\boldsymbol{r}(y^{+})=\boldsymbol{y}}_{\boldsymbol{r}(x^{+})=\boldsymbol{x}}\mathcal{D}\boldsymbol{r}(\xi)\exp{\left\\{i\frac{p^{+}}{2}\int_{x^{+}}^{y^{+}}ds^{+}\,\dot{\boldsymbol{r}}^{2}\right\\}}U^{ij}(x^{+},y^{+},\boldsymbol{r}(\xi))\,,$
(2)

with the endpoints denoting the starting and final transverse positions of the
quark evolution in $x^{+}$ at fixed light-cone energy $p^{+}$. We assume
$p^{+}\gg|{\boldsymbol{p}}|$ to be the largest scale in the problem, such that
the evolution is highly collimated and all the dynamics are constrained to the
plane transverse to the jet axis.555Of course, particle branching leads to
energy sharing between the daughter partons. The successive scatterings with
the medium lead to color precession, which is captured by the path-ordered
fundamental Wilson line

$U^{ij}(x^{+},y^{+},\boldsymbol{r})=\mathcal{P}\exp{\left\\{ig\int_{x^{+}}^{y^{+}}ds^{+}\,\mathcal{A}_{a}^{-}(s^{+},\boldsymbol{r}(s^{+}))t^{a}_{ij}\right\\}}\,.$
(3)

Here $\mathcal{A}$ denote the stochastic background field representing the
underlying QCD medium. When computing any observable one should average over
field configurations, which we assume follow a Gaussian distribution

$\langle\mathcal{A}_{a}^{-}(x^{+},x^{-},\boldsymbol{x})\mathcal{A}_{b}^{-}(y^{+},y^{-},\boldsymbol{y})\rangle=\delta_{ab}n(x^{+})\delta(x^{+}-y^{+})\gamma(\boldsymbol{x-y})\,,$
(4)

with all correlations being local in color, time, and transverse position.
Also, $\gamma({\boldsymbol{x}})$ is directly related to the in-medium elastic
scattering rate. Combined with the effective Feynman rules detailed in
Appendix A, one can compute processes such as the one in Fig. 1 resumming all
possible gluon exchanges with the medium at leading order in eikonality.

The amplitude describing the top diagram in Fig. 1 reads

$\displaystyle i\mathcal{M}^{hh^{\prime}}$
$\displaystyle=\frac{1}{2q_{0}^{+}}e^{i\frac{\boldsymbol{p_{1}}^{2}}{2p_{1}^{+}}L^{+}}e^{i\frac{\boldsymbol{p_{2}}^{2}}{2p_{2}^{+}}L^{+}}\sum_{\lambda=\pm
1}\int_{x_{v}^{+}}\int_{\boldsymbol{x_{v},{\boldsymbol{x}}_{q},{\boldsymbol{x}}_{\bar{q}},{\boldsymbol{x}}_{g},{\boldsymbol{x}},{\boldsymbol{y}}}}\int_{x_{g}^{+}}\mathcal{M}_{0}^{\lambda,a}(q_{0}^{+},x_{g}^{+},{\boldsymbol{x}}_{g})e^{-i\boldsymbol{p_{1}\cdot
x_{q}}}e^{-i\boldsymbol{p_{2}\cdot x_{\bar{q}}}}$
$\displaystyle\times\overline{\mathcal{G}}^{lj}(L^{+},\boldsymbol{x_{\bar{q}}};x_{v}^{+},\boldsymbol{x_{v}}|p_{2}^{+})\left(igt_{kl}^{b}V^{\lambda
hh^{\prime}}(z,\boldsymbol{{\boldsymbol{x}}},\boldsymbol{{\boldsymbol{y}}})\right)\mathcal{G}^{ba}_{A}(x_{v}^{+},\boldsymbol{x_{v}-{\boldsymbol{x}}};x_{g}^{+},\boldsymbol{x_{g}}|q_{0}^{+})$
$\displaystyle\times\mathcal{G}^{ik}(L^{+},\boldsymbol{{\boldsymbol{x}}_{q}};x_{v}^{+},\boldsymbol{x_{v}-{\boldsymbol{y}}}|p_{1}^{+})\,,$
(5)

where

$\displaystyle V^{\lambda
hh^{\prime}}(z,{\boldsymbol{x}},{\boldsymbol{y}})=2\gamma^{\lambda
h}(z)\delta_{h,-h^{\prime}}\boldsymbol{\epsilon}_{\lambda}\cdot
i(\delta({\boldsymbol{x}})\delta^{{}^{\prime}}({\boldsymbol{y}})+z\delta^{{}^{\prime}}({\boldsymbol{x}})\delta({\boldsymbol{y}}))\,,$
$\displaystyle\mathcal{M}_{0}^{\lambda,a}(q_{0}^{+},x_{g}^{+},{\boldsymbol{x}}_{g})=\int_{p^{-},{\boldsymbol{p}}}e^{-ip^{-}x_{g}^{+}}e^{i{\boldsymbol{p}}\cdot{\boldsymbol{x}}_{g}}\mathcal{M}_{0}^{\mu,a}(p)\epsilon^{\ast}_{\lambda,\mu}(p)\,,$
(6)

and $\mathcal{M}_{0}^{\mu,a}(p)$ is the momentum space current that produces
the gluon. We impose a fixed initial energy for the gluon $p^{+}=q_{0}^{+}$,
and we shall denote the light cone energy fraction carried by the quark by
$z\equiv p_{1}^{+}/q_{0}^{+}$. The notation used for the integrals is, for
example:

$\int_{x_{v}^{+}}\int_{\boldsymbol{x_{v}}}\int_{{\boldsymbol{p}}}=\int_{x_{0}^{+}}^{L^{+}}dx_{v}^{+}\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}d^{2}\boldsymbol{x_{v}}\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\frac{d^{2}{\boldsymbol{p}}}{(2\pi)^{2}}\,.$
(7)

In what follows we shall take $x_{0}^{+}=0$.

The vertex $V^{\lambda hh^{\prime}}$ function in Eq. (2) is written in the
circular polarization basis and as a function of the transverse polarization
of the initial gluon $\lambda=\pm 1$ and of the helicities of the quark
($h=\pm 1$) and anti-quark ($h^{\prime}=\mp 1$) (cf. Appendix A). The basis
for polarization states is, however, not relevant since, as discussed below,
we will average over them. As for the fermions’ helicity states, these will be
discussed further ahead and we keep them explicit since we aim at computing
the helicity dependent cross-section. To obtain this usual form for
$V^{\lambda hh^{\prime}}$, we applied the completeness relation
$\sum_{\lambda=\pm
1}\epsilon_{\lambda,i}({\boldsymbol{p}})\epsilon_{\lambda,\mu}^{\ast}=d_{i\mu}({\boldsymbol{p}})$
to one of the two transverse projectors making up the Lorentz structure of the
gluon propagator Blaizot:2012fh . To obtain a more familiar expression for the
vertex function involving transverse space derivatives, one should recall the
property

$\displaystyle\int_{-\infty}^{\infty}dy\hskip
2.84526pt\delta^{\prime}(y)f(x-y)=f^{\prime}(x)\,,$ (8)

which is valid for an arbitrary smooth test function $f$. In the coordinate-
space notation used in Eq. (2), the transverse coordinate displacements are
associated with the gluon ($\boldsymbol{{\boldsymbol{x}}}$) and the quark
($\boldsymbol{{\boldsymbol{y}}}$) propagators. The derivatives of
$\delta$-functions appearing in $V^{\lambda hh^{\prime}}$ are transverse
vectors which are contracted with the gluon’s polarization vector
$\boldsymbol{\epsilon}_{\lambda}=\frac{1}{\sqrt{2}}(1,\lambda i)$.

In order to factorize the initial current in light-cone time we assume

$\int_{p^{-}}e^{-ip^{-}x_{g}^{+}}\mathcal{M}_{0}^{\mu,a}(p)=\mathcal{M}_{0}^{\mu,a}(q_{0}^{+},x_{g}^{+},{\boldsymbol{p}})\delta(x_{g}^{+}-x_{0}^{+})\,,$
(9)

which means the gluon is produced at the initial time of the medium
$x_{0}^{+}$ with a fixed energy $q_{0}^{+}$. We will stop explicitly writing
the initial current’s dependence on $q_{0}^{+}$ and $x_{0}^{+}$ from now on.
Also, when squaring the amplitude we take a gluon with a randomised
polarization and color (which is fixed to be the same in the conjugate
amplitude)

$\begin{split}&\sum_{\lambda,\bar{\lambda}=\pm
1}\mathcal{M}_{0}^{\mu,a}({\boldsymbol{p}})\epsilon_{\lambda,\mu}^{\ast}\left(\mathcal{M}_{0}^{\bar{\mu},\bar{a}}(\bar{\boldsymbol{p}})\epsilon_{\bar{\lambda},\bar{\mu}}^{\ast}\right)^{\dagger}\rightarrow\sum_{\lambda,\bar{\lambda}=\pm
1}\frac{\delta^{a\bar{a}}}{N_{c}^{2}-1}\frac{\delta^{\lambda\bar{\lambda}}}{2}\mathcal{M}_{0}({\boldsymbol{p}})\mathcal{M}_{0}^{\ast}(\bar{\boldsymbol{p}})\,,\end{split}$
(10)

where the correction factors below each Kronecker-$\delta$ account for the
average over the fixed choice of quantum numbers. Hence, under these
assumptions that impose the production at fixed time of an initial gluon with
fixed light-cone energy and quantum numbers, we can work with the following
amplitude

$\displaystyle i\mathcal{M}^{\lambda h}$
$\displaystyle=\frac{1}{2q_{0}^{+}}e^{i\frac{\boldsymbol{p_{1}}^{2}}{2p_{1}^{+}}L^{+}}e^{i\frac{\boldsymbol{p_{2}}^{2}}{2p_{2}^{+}}L^{+}}\int_{{\boldsymbol{p}}}\mathcal{M}_{0}({\boldsymbol{p}})\int_{x_{v}^{+}}\int_{\boldsymbol{x_{v},{\boldsymbol{x}}_{q},{\boldsymbol{x}}_{\bar{q}},{\boldsymbol{x}}_{g},{\boldsymbol{x}},{\boldsymbol{y}}}}e^{i\boldsymbol{p\cdot
x_{g}}}e^{-i\boldsymbol{p_{1}\cdot x_{q}}}e^{-i\boldsymbol{p_{2}\cdot
x_{\bar{q}}}}$
$\displaystyle\times\overline{\mathcal{G}}^{lj}(L^{+},\boldsymbol{x_{\bar{q}}};x_{v}^{+},\boldsymbol{x_{v}}|p_{2}^{+})\left(igt_{kl}^{b}V^{\lambda
h}(z,\boldsymbol{{\boldsymbol{x}}},\boldsymbol{{\boldsymbol{y}}})\right)\mathcal{G}^{ba}_{A}(x_{v}^{+},\boldsymbol{x_{v}-{\boldsymbol{x}}};x_{0}^{+},\boldsymbol{x_{g}}|q_{0}^{+})$
$\displaystyle\times\mathcal{G}^{ik}(L^{+},\boldsymbol{{\boldsymbol{x}}_{q}};x_{v}^{+},\boldsymbol{x_{v}-{\boldsymbol{y}}}|p_{1}^{+})\,,$
(11)

as long as we average over the initial gluon’s quantum numbers when squaring.
Since the vacuum QCD vertex imposes $h^{\prime}=-h$, we dropped the explicit
dependence on the anti-quark’s helicity $h^{\prime}$.

The color and medium averaged double differential cross-section can be written
as

$4z(1-z)(2\pi)^{5}\sigma_{0}\frac{dN^{\lambda
h}}{dzd^{2}\boldsymbol{p_{1}}d^{2}\boldsymbol{p_{2}}}=2\text{Re}\,\left[\frac{1}{N_{c}^{2}-1}\sum_{\rm
colors}\langle\mathcal{M}\mathcal{M}^{\dagger}\rangle\right]_{\overline{x}_{v}^{+}>x_{v}^{+}}\,,$
(12)

where $\sigma_{0}$ denotes the initial cross-section to produce the gluon:

$\sigma_{0}=\int_{{\boldsymbol{p}}}\left|\mathcal{M}_{0}({\boldsymbol{p}})\right|^{2}\,.$
(13)

We note that with this choice of normalization, the vacuum limit for this
process reads666Note here the extra $1/2$ factor compared to the usual
normalization for the vacuum piece; this can be reabsorbed in the definition
of $\sigma_{0}$.

$\displaystyle\frac{dN^{h}}{dzd\theta_{q\bar{q}}d\phi}=\frac{\alpha_{s}}{(2\pi)^{2}}\frac{P_{qg}(z)}{\theta_{q\bar{q}}}\,.$
(14)

The local form of the medium average in the $+$ component in Eq. (4) greatly
simplifies the structure of the calculation, since it allows to split the
overall average into the product of shorter time averages. Using the
composition rule for the in-medium propagators

$\mathcal{G}^{ik}(L^{+},\boldsymbol{x_{q}};x_{v}^{+},\boldsymbol{x_{v}-y}|p_{1}^{+})=\int_{\boldsymbol{z_{3}}}\mathcal{G}^{ii_{1}}(L^{+},\boldsymbol{x_{q}};\overline{x}_{v}^{+},\boldsymbol{z_{3}}|p_{1}^{+})\mathcal{G}^{i_{1}k}(\overline{x}_{v}^{+},\boldsymbol{z_{3}};x_{v}^{+},\boldsymbol{x_{v}-y}|p_{1}^{+})\,,$
(15)

allowing to directly split the medium averages in Eq. (12) into three regions:
$(x_{0}^{+},x_{v}^{+})$, $(x_{v}^{+},\overline{x}_{v}^{+})$ and
$(\overline{x}_{v}^{+},L^{+})$. This results in the expression for the cross-
section:

$\begin{split}&4z(1-z)(2\pi)^{5}\frac{dN^{\lambda
h}}{dzd^{2}\boldsymbol{p_{1}}d^{2}\boldsymbol{p_{2}}}=\frac{2g^{2}}{(2q_{0}^{+})^{2}}\text{Re}\,\int_{{\boldsymbol{p}},\bar{\boldsymbol{p}}}\mathcal{M}_{0}({\boldsymbol{p}})\mathcal{M}_{0}(\bar{\boldsymbol{p}})^{\ast}\\\
&\times\int_{X_{v}^{+}}\int_{\boldsymbol{X}}e^{i\boldsymbol{p\cdot
x_{g}}}e^{-i\boldsymbol{\bar{p}\cdot\overline{x}_{g}}}e^{-i\boldsymbol{p_{1}\cdot(x_{q}-\overline{x}_{q})}}e^{-i\boldsymbol{p_{2}\cdot(x_{\bar{q}}-\overline{x}_{\bar{q}})}}V^{\lambda
h}(z,\boldsymbol{x},\boldsymbol{y})\left(V^{\lambda
h}(z,\boldsymbol{\overline{x}},\boldsymbol{\overline{y}})\right)^{\ast}\\\\[8.61108pt]
&\times\frac{1}{N_{c}^{2}-1}\left\langle\mathcal{G}^{ba}_{A}(x_{v}^{+},\boldsymbol{x_{v}-x};x_{0}^{+},\boldsymbol{x_{g}}|q_{0}^{+})(\mathcal{G}^{\dagger})^{a\overline{b}_{1}}_{A}(x_{v}^{+},\boldsymbol{z_{1}};x_{0}^{+},\boldsymbol{\overline{x}_{g}}|q_{0}^{+})\right\rangle\\\
&\times
t_{kl}^{b}t_{\overline{l}\overline{k}}^{\overline{b}}\left\langle(\mathcal{G}^{\dagger})^{\overline{b}_{1}\overline{b}}_{A}(\overline{x}_{v}^{+},\boldsymbol{\overline{x}_{v}-\overline{x}};x_{v}^{+},\boldsymbol{z_{1}}|q_{0}^{+})\overline{\mathcal{G}}^{lj_{1}}(\overline{x}_{v}^{+},\boldsymbol{z_{2}};x_{v}^{+},\boldsymbol{x_{v}}|p_{2}^{+})\mathcal{G}^{i_{1}k}(\overline{x}_{v}^{+},\boldsymbol{z_{3}};x_{v}^{+},\boldsymbol{x_{v}-y}|p_{1}+)\right\rangle\\\
&\times\left\langle\mathcal{G}^{ii_{1}}(L^{+},\boldsymbol{x_{q}};\overline{x}_{v}^{+},\boldsymbol{z_{3}}|p_{1}^{+})(\mathcal{G}^{\dagger})^{\overline{k}i}(L^{+},\boldsymbol{\overline{x}_{q}};\overline{x}_{v}^{+},\boldsymbol{\overline{x}_{v}-\overline{y}}|p_{1}^{+})\right.\\\
&\times\left.\overline{\mathcal{G}}^{j_{1}j}(L^{+},\boldsymbol{x_{\bar{q}}};\overline{x}_{v}^{+},\boldsymbol{z_{2}}|p_{2}^{+})(\overline{\mathcal{G}}^{\dagger})^{j\overline{l}}(L^{+},\boldsymbol{\overline{x}_{\bar{q}}};\overline{x}_{v}^{+},\boldsymbol{\overline{x}_{v}}|p_{2}^{+})\right\rangle\,,\\\
\end{split}$ (16)

where, for conciseness, we write
$\boldsymbol{X}=\\{\boldsymbol{x_{v},\overline{x}_{v},x_{q},\overline{x}_{q},x_{\bar{q}},\overline{x}_{\bar{q}},x_{g},\overline{x}_{g},x,\overline{x},y,\overline{y},z_{1},z_{2},z_{3}}\\}$
and $X_{v}^{+}=\\{x_{v}^{+},\overline{x}_{v}^{+}>x_{v}^{+}\\}$.

To further simplify Eq. (16), we perform the remaining two, three and four
body averages. To that end, one should project the above object in both the
singlet and octet color subspaces. However, since we do not consider the
presence of coherent color fields in the medium, the octet contribution should
be suppressed and we leave its study to future work, see e.g. Barata:2023uoi ;
Blaizot:2017ypk ; Blaizot:1996az ; Braaten:1990it  for related discussions.
With this in mind, in each time region the system must be in a color singlet
state, since Eq. (4) is local in color space. As a consequence of this
statement, we can factorize the color structure in each time region and
contract them as detailed in Appendix B.1. Using this, the four last lines in
Eq. (16) can be written as a product
$\mathcal{S}^{(2)}\mathcal{S}^{(3)}\mathcal{S}^{(4)}$, corresponding to the
same two, three and four body averages but now projected onto a color singlet
state and fully contracted in color. The object $\mathcal{S}^{(2)}$ describes
the gluon broadening in the first region and $\mathcal{S}^{(3)}$ the
production of the $q\bar{q}$ pair antenna. The color structure in
$\mathcal{S}^{(4)}$ is further simplified by applying the large-$N_{c}$ limit,
in which case it describes the independent broadening of the quark and anti-
quark. In Appendix B.2 and B.3 we further simplify these objects by
analytically solving the resulting path integrals either fully
($\mathcal{S}^{(2)}$ and $\mathcal{S}^{(4)}$) or partially
($\mathcal{S}^{(3)}$). After carrying out these calculations, an emission
kernel naturally appears (cf. Appendix B.3)

$\begin{split}\mathcal{K}(\boldsymbol{u}_{1}^{i},&\boldsymbol{u}_{1}^{f},\boldsymbol{u}_{2}^{i},\boldsymbol{u}_{2}^{f})\equiv\int\mathcal{D}\boldsymbol{u}_{1}\exp{\left\\{\frac{i\tilde{q}_{0}^{+}}{2}\int
ds^{+}\dot{\boldsymbol{u}}_{1}^{2}\right\\}}\\\
&\times\exp{\left\\{-\frac{1}{2}\int
ds^{+}\,n(s^{+})(N_{c}(\sigma(\boldsymbol{u}_{2}^{0}+z\boldsymbol{u}_{1})+\sigma(\boldsymbol{u}_{2}^{0}-(1-z)\boldsymbol{u}_{1}))-\frac{1}{N_{c}}\sigma(-\boldsymbol{u}_{1}))\right\\}}\\\
&\equiv\int_{\boldsymbol{t}_{1},\boldsymbol{t}_{2},\boldsymbol{t}_{3},\boldsymbol{t}_{4}}e^{i\boldsymbol{t}_{1}\cdot\boldsymbol{u}_{1}^{i}}e^{i\boldsymbol{t}_{2}\cdot\boldsymbol{u}_{1}^{f}}e^{i\boldsymbol{t}_{3}\cdot\boldsymbol{u}_{2}^{i}}e^{i\boldsymbol{t}_{4}\cdot\boldsymbol{u}_{2}^{f}}\tilde{\mathcal{K}}(\boldsymbol{t}_{1},\boldsymbol{t}_{2},\boldsymbol{t}_{3},\boldsymbol{t}_{4})\,,\end{split}$
(17)

where $\tilde{q}_{0}^{+}=z(1-z)q_{0}^{+}$ and the full cross-section then
reads

$\begin{split}4z(1-z)(2\pi)^{5}&\frac{dN^{\lambda
h}}{dzd^{2}\boldsymbol{p_{1}}d^{2}\boldsymbol{p_{2}}}=\frac{g^{2}}{(2q_{0}^{+})^{2}}\text{Re}\,\int_{{\boldsymbol{p}}}|\mathcal{M}_{0}({\boldsymbol{p}})|^{2}\int_{X_{v}^{+}}\int_{\boldsymbol{v_{1},v_{2},v_{3}}}\\\
&\times
e^{i{\boldsymbol{p}}\cdot\boldsymbol{v_{1}}}\exp{\left\\{-\frac{C_{A}}{2}\int_{x_{0}^{+}}^{x_{v}^{+}}ds^{+}\,n(s^{+})\sigma(\boldsymbol{v_{1}})\right\\}}\\\
&\times e^{-i\boldsymbol{p_{1}\cdot v_{3}}}e^{-i\boldsymbol{p_{2}\cdot
v_{2}}}\exp{\left\\{-\frac{C_{F}}{2}\int_{\bar{x}_{v}^{+}}^{L^{+}}ds^{+}\,n(s^{+})(\sigma(\boldsymbol{v_{3}})+\sigma(\boldsymbol{-v_{2}}))\right\\}}\\\
&\times\int_{\boldsymbol{t}_{1},\boldsymbol{t}_{2},\boldsymbol{t}_{3},\boldsymbol{t}_{4}}\tilde{\mathcal{K}}(\boldsymbol{t}_{1},\boldsymbol{t}_{2},\boldsymbol{t}_{3},\boldsymbol{t}_{4})\int_{\boldsymbol{y,\overline{y},x,\overline{x},\Delta_{v}}}V^{\lambda
h}(z,\boldsymbol{x},\boldsymbol{y})\left(V^{\lambda
h}(z,\boldsymbol{\overline{x}},\boldsymbol{\overline{y}})\right)^{\ast}\\\
&\times\left(\frac{q_{0}^{+}}{2\pi\Delta
t}\right)^{2}\exp{\left\\{\frac{iq_{0}^{+}}{2}\frac{\Delta\boldsymbol{u}_{2}}{\Delta
t}\left(2\Delta\boldsymbol{r}_{0}+\Delta\boldsymbol{u}_{2}\right)\right\\}}e^{i\boldsymbol{t}_{1}\cdot\boldsymbol{u}_{1}^{i}}e^{i\boldsymbol{t}_{2}\cdot\boldsymbol{u}_{1}^{f}}e^{i\boldsymbol{t}_{3}\cdot\boldsymbol{u}_{2}^{i}}e^{i\boldsymbol{t}_{4}\cdot\boldsymbol{u}_{2}^{f}}\,.\end{split}$
(18)

Here we used the notation
$\Delta\boldsymbol{X}=\boldsymbol{X^{f}}-\boldsymbol{X^{i}}$ and, for
completeness, we have
$\Delta{\boldsymbol{r}}_{0}={\boldsymbol{x}}+\boldsymbol{v_{1}}-\bar{\boldsymbol{x}}+\boldsymbol{\Delta_{v}}$
and

$\begin{split}&\boldsymbol{u_{1}^{i}}=\boldsymbol{y}\\\
&\boldsymbol{u_{2}^{i}}=\boldsymbol{v_{1}}+{\boldsymbol{x}}-z\boldsymbol{y}\end{split}\qquad\qquad\begin{split}&\boldsymbol{u_{1}^{f}}=\boldsymbol{v_{2}}-\boldsymbol{v_{3}}+\boldsymbol{\overline{y}}\\\
&\boldsymbol{u_{2}^{f}}=\boldsymbol{v_{2}}+\overline{{\boldsymbol{x}}}-z(\boldsymbol{v_{2}}-\boldsymbol{v_{3}}+\boldsymbol{\overline{y}})\,.\end{split}$
(19)

Note that we have changed variables from $\bar{\boldsymbol{x}}_{v}$ to
$\boldsymbol{\Delta_{v}}=\bar{\boldsymbol{x}}_{v}-{\boldsymbol{x}}_{v}$ and
integrated over ${\boldsymbol{x}}_{v}$ and $\bar{\boldsymbol{p}}$, as
explained at the end of Appendix B.3. This sets
${\boldsymbol{p}}=\bar{\boldsymbol{p}}$, i.e., the same transverse momentum
flows out of the initial current in amplitude and its complex conjugate.

Let us take note on an important aspect regarding phase space dependence. It
is an intuitive idea that one should keep the spectrum differential in some
directional quantity in order to be able to probe anisotropy induced
modifications. However, as we will see, integrating over all momenta still
leaves trace of the medium anisotropy. Despite this, we choose to keep the
spectrum differential in some angle(s) which we will define ahead, since this
provides an additional handle on probing anisotropic effects. In fact, keeping
some differential dependence on an angle in the transverse plane is necessary
to keep spin/helicity dependence. To that end, we switch our phase space
variables from $\\{\boldsymbol{p_{1}},\boldsymbol{p_{2}}\\}$ to
$\\{\boldsymbol{P}^{CM},\boldsymbol{P}^{rel}\\}$, where the center of
transverse momentum and the relative transverse momentum are, respectively,
given by

$\displaystyle\boldsymbol{P}^{CM}=\boldsymbol{p_{1}}+\boldsymbol{p_{2}}\,,\quad\boldsymbol{P}^{rel}=(1-z)\boldsymbol{p_{1}}-z\boldsymbol{p_{2}}\,,$
(20)

which is a transformation with unit jacobian. $\boldsymbol{P}^{rel}$
($-\boldsymbol{P}^{rel}$) is the momentum of the quark (anti-quark) in the
center of momentum frame.

We now carry out the derivatives in the vertex, integrate out the center of
momentum information $\boldsymbol{P}^{CM}$ and at the same time average over
gluon polarizations $\lambda$. By using the massless, helicity-dependent
vertex in Eq. (2), this leads us to:

$\begin{split}&4z(1-z)(2\pi)^{3}\frac{dN^{h}}{dzd^{2}\boldsymbol{P}^{rel}}=-\frac{2g^{2}}{(2q_{0}^{+})^{2}z(1-z)}\text{Re}\,\int_{X_{v}^{+}}\int_{\boldsymbol{v}}\int_{\boldsymbol{t}_{1},\boldsymbol{t}_{2}}\\\
&\times
e^{-i(\boldsymbol{P}^{rel}+\boldsymbol{t}_{2})\cdot\boldsymbol{v}}\exp{\left\\{-\frac{C_{F}}{2}\int_{\bar{x}_{v}^{+}}^{L^{+}}ds^{+}\,n(s^{+})(\sigma((1-z)\boldsymbol{v})+\sigma(z\boldsymbol{v}))\right\\}}\\\
&\times\left[P_{qg}^{vac}(z)\left(\boldsymbol{t}_{1}\cdot\boldsymbol{t}_{2}\right)+\frac{ih}{2}(1-2z)\left(\boldsymbol{t}_{1}\times\boldsymbol{t}_{2}\right)_{z}\right]\tilde{\mathcal{K}}(\boldsymbol{t}_{1},\boldsymbol{t}_{2})\,,\end{split}$
(21)

where

$\begin{split}&\tilde{\mathcal{K}}(\boldsymbol{t}_{1},\boldsymbol{t}_{2})=\int_{\boldsymbol{a_{1}},\boldsymbol{a_{2}}}e^{-i\boldsymbol{a_{1}}\cdot\boldsymbol{t}_{1}}e^{-i\boldsymbol{a_{2}}\cdot\boldsymbol{t}_{2}}\int_{\boldsymbol{a_{1}}}^{\boldsymbol{a_{2}}}\mathcal{D}\boldsymbol{u}\exp{\left\\{\frac{i\tilde{q}_{0}^{+}}{2}\int_{x_{v}^{+}}^{\bar{x}_{v}^{+}}ds^{+}\dot{\boldsymbol{u}}^{2}\right\\}}\\\
&\times\exp{\left\\{-\frac{1}{2}\int_{x_{v}^{+}}^{\bar{x}_{v}^{+}}ds^{+}\,n(s^{+})(N_{c}(\sigma(z\boldsymbol{u})+\sigma(-(1-z)\boldsymbol{u}))-\frac{1}{N_{c}}\sigma(-\boldsymbol{u}))\right\\}}\,,\end{split}$
(22)

and
$\left(\boldsymbol{t}_{1}\times\boldsymbol{t}_{2}\right)_{z}=\boldsymbol{t}_{1x}\boldsymbol{t}_{2y}-\boldsymbol{t}_{1y}\boldsymbol{t}_{2x}$.
Note that this integration in $\boldsymbol{P}^{CM}$ allowed us to factorize
the integral over ${\boldsymbol{p}}$ of the initial current squared
($\sigma_{0}$ in Eq. (13)) from the remaining of the amplitude, so we have
dropped it. This integration also resulted in the gluon broadening being
completely integrated out.

At this point we have to provide an explicit form for the dipole cross-section
$\sigma$. Here we use the harmonic approximation Zakharov:1996fv , whereby the
logarithmic dependence on the transverse position, in the small distance
limit, is dropped in $\sigma$, see e.g. Mehtar-Tani:2019tvy ; Barata:2020sav
for further discussion on the limitations of this approximation. Additionally
assuming that the medium is static, we can write

$N_{c}\int_{x^{+}}^{y^{+}}ds^{+}\,n(s^{+})\sigma(\boldsymbol{r})=N_{c}(y^{+}-x^{+})n\sigma(\boldsymbol{r})=\frac{(y^{+}-x^{+})}{2}\left(\hat{q}_{x}r_{x}^{2}+\hat{q}_{y}r_{y}^{2}\right)\,,$
(23)

where we have written the components of the transverse vector as
${\boldsymbol{r}}=(r_{x},r_{y})$. Furthermore, we have naturally introduced a
form of medium anisotropy by considering different values for the jet
quenching parameter $\hat{q}$ in different directions in the transverse plane
Hauksson:2023tze . Combining this with Eq. (21), we obtain

$\begin{split}&4z(1-z)(2\pi)^{3}\frac{dN^{h}}{dzd^{2}\boldsymbol{P}^{rel}}=\frac{2g^{2}}{(2q_{0}^{+})^{2}z(1-z)}\\\
&\times\text{Re}\,\int_{X_{v}^{+}}\int_{\boldsymbol{v}}e^{-i\boldsymbol{P}^{rel}\cdot\boldsymbol{v}}\exp{\left\\{-\frac{C_{F}}{2N_{c}}\Delta_{L}^{+}P_{qg}^{vac}(z)\left(\hat{q}_{x}v_{x}^{2}+\hat{q}_{y}v_{y}^{2}\right)\right\\}}\\\\[4.30554pt]
&\times\left[P_{qg}^{vac}(z)\left(\boldsymbol{\nabla_{a_{1}}}\cdot\boldsymbol{\nabla_{a_{2}}}\right)+\frac{ih}{2}(1-2z)\left(\boldsymbol{\nabla_{a_{1}}}\times\boldsymbol{\nabla_{a_{2}}}\right)_{z}\right]\mathcal{K}(\boldsymbol{a_{1}},\boldsymbol{a_{2}})\Big{|}_{\boldsymbol{a_{1}}=0,\boldsymbol{a_{2}}=-\boldsymbol{v}}\,,\end{split}$
(24)

where $\Delta_{L}^{+}=L^{+}-\overline{x}_{v}^{+}$. Importantly, note that if
the medium was isotropic there would be no helicity ($h$) dependence attached
to medium parameters $\hat{q}_{x},\hat{q}_{y}$. The path integral
$\mathcal{K}(\boldsymbol{a_{1}},\boldsymbol{a_{2}})$ has a well known solution
in the literature Apolinario:2014csa ; Kleinert:2004ev

$\begin{split}\mathcal{K}(\boldsymbol{a_{1}},\boldsymbol{a_{2}})&=\sqrt{\frac{\tilde{q}_{0}^{+}\Omega_{x}}{2\pi
i\sin{\Omega_{x}\Delta
t}}}\exp{\left\\{i\frac{\tilde{q}_{0}^{+}\Omega_{x}}{2\sin{\Omega_{x}\Delta
t}}\left[\left(a_{1x}^{2}+a_{2x}^{2}\right)\cos{\Omega_{x}\Delta
t}-2a_{1x}a_{2x}\right]\right\\}}\\\
&\times\sqrt{\frac{\tilde{q}_{0}^{+}\Omega_{y}}{2\pi i\sin{\Omega_{y}\Delta
t}}}\exp{\left\\{i\frac{\tilde{q}_{0}^{+}\Omega_{y}}{2\sin{\Omega_{y}\Delta
t}}\left[\left(a_{1y}^{2}+a_{2y}^{2}\right)\cos{\Omega_{y}\Delta
t}-2a_{1y}a_{2y}\right]\right\\}}\,,\end{split}$ (25)

with

$\Omega_{i}=\left(\frac{1-i}{\sqrt{2}}\right)\omega_{i},\qquad\omega_{i}\approx\sqrt{\frac{\hat{q}_{i}P_{qg}^{vac}(z)}{\tilde{q}_{0}^{+}}},\qquad\Delta
t=\overline{x}_{v}^{+}-x_{v}^{+}\,,$ (26)

where $P_{qg}^{vac}(z)=T_{R}(z^{2}+(1-z)^{2})$ is the vacuum splitting
function, we took $C_{F}/N_{c}\approx 1/2$ and we ignored a $1/N_{c}^{2}$ term
in $\omega_{i}$.

We now choose to write the result in terms of two angles: the angle between
the 3-momenta $\vec{p}_{1}$ and $\vec{p_{2}}$ and the azimuthal angle of the
relative transverse momentum $\boldsymbol{P}^{rel}$. These are given by

$\displaystyle|\boldsymbol{P}^{rel}|=\sin\left(\frac{\theta_{q\bar{q}}}{2}\right)\tilde{q}_{0}^{+}\sqrt{2}\approx\frac{\theta_{q\bar{q}}}{\sqrt{2}}\tilde{q}_{0}^{+}\,,\quad\phi=\arctan\frac{{\boldsymbol{P}}^{rel}_{y}}{{\boldsymbol{P}}^{rel}_{x}}\,,$
(27)

where the expression for $\theta_{q\overline{q}}$ assumes
$(|{\boldsymbol{p}}_{i}|/p_{i}^{+})^{2}\ll 1$. Using the results in Eq. (24)
and (25), one can then write the differential cross-section as a function of
these two angles as

$\displaystyle\frac{dN^{h}}{dzd\theta_{q\bar{q}}d\phi}=\frac{\alpha_{s}\theta_{q\bar{q}}}{8\pi^{2}}\text{Im}\int_{X_{v}^{+}}\frac{\sqrt{c_{1x}}\sqrt{c_{1y}}}{\sqrt{c_{3x}}\sqrt{c_{3y}}}\exp{\left\\{i\frac{\tilde{q}_{0}^{+}\theta^{2}_{q\bar{q}}}{4}\left(\frac{\cos^{2}\phi}{c_{3x}}+\frac{\sin^{2}\phi}{c_{3y}}\right)\right\\}}$
$\displaystyle\times\Bigg{[}P_{qg}^{vac}(z)\Bigg{(}\tilde{q}_{0}^{+}c_{4}+\tilde{q}_{0}^{+}\left(\frac{c_{1x}c_{2x}}{c_{3x}}+\frac{c_{1y}c_{2y}}{c_{3y}}\right)+i\frac{(\tilde{q}_{0}^{+})^{2}\theta^{2}_{q\bar{q}}}{2}\left(\frac{c_{1x}c_{2x}}{c_{3x}^{2}}\cos^{2}\phi+\frac{c_{1y}c_{2y}}{c_{3y}^{2}}\sin^{2}\phi\right)\Bigg{)}$
$\displaystyle+\frac{h(1-2z)}{4}(\tilde{q}_{0}^{+})^{2}\theta^{2}_{q\bar{q}}\frac{c_{5}\cos\phi\sin\phi}{c_{3x}c_{3y}}\Bigg{]}\,,$
(28)

where we have introduced the compact notation

$\begin{split}&c_{1i}=\frac{\Omega_{i}}{2i\sin{\Omega_{i}\Delta
t}},\qquad\qquad c_{2i}=\frac{\Omega_{i}}{\tan{\Omega_{i}\Delta t}},\qquad
c_{3i}=\frac{\Delta_{L}^{+}\Omega_{i}^{2}}{2}-c_{2i},\\\
&c_{4}=\left(c_{1x}+c_{1y}\right),\hskip 28.45274pt\qquad
c_{5}=\left(c_{1y}c_{2x}-c_{1x}c_{2y}\right)\,.\end{split}$ (29)

We should emphasize that if we were to integrate over $\boldsymbol{P}^{rel}$
in Eq. (24) the only term left would be the one proportional to $c_{4}$, which
is nevertheless sensitive to an anisotropic $\hat{q}$. That term would be
similar to the one obtained for the decay rate of a $y$-polarized gluon to two
$z-$polarized gluons in Hauksson:2023tze , in which case a (linear)
polarization dependence exists even after integrating over all momenta.
Finally, the purely azimuthal distribution is given by777It is worth pointing
out that $\text{Im}(c_{3i})<0$ for all $x_{v}^{+}$ and $\bar{x}_{v}^{+}$
guarantees the convergence of both integrals in $\boldsymbol{v}$ and
$|\boldsymbol{P}^{rel}|$ and is responsible for the sign attached to the
second term in Eq. (2).

$\begin{split}&\frac{dN^{h}}{dzd\phi}=\frac{\alpha_{s}}{(2\pi)^{2}}\text{Re}\,\int_{X_{v}^{+}}\frac{\sqrt{c_{1x}}\sqrt{c_{1y}}}{\sqrt{c_{3x}}\sqrt{c_{3y}}}\left(P_{qg}^{vac}(z)\Bigg{(}\frac{(c_{1y}c_{2y}c_{3x}-c_{1x}c_{2x}c_{3y})(c_{3y}\cos^{2}{\phi}-c_{3x}\sin^{2}{\phi})}{\left(c_{3y}\cos^{2}{\phi}+c_{3x}\sin^{2}{\phi}\right)^{2}}\right.+\\\
&\left.+\frac{c_{3x}c_{3y}c_{4}}{c_{3y}\cos^{2}{\phi}+c_{3x}\sin^{2}{\phi}}\Bigg{)}+ih(1-2z)\frac{c_{3x}c_{3y}c_{5}\sin{\phi}\cos{\phi}}{\left(c_{3y}\cos^{2}{\phi}+c_{3x}\sin^{2}{\phi}\right)^{2}}\right)\,.\end{split}$
(30)

Note that the distribution in Eq. (30) has a linear divergence $1/\Delta t$ as
$\Delta t\rightarrow 0$, as noted in Hauksson:2023tze , which is removed by
subtracting a contribution corresponding to $\hat{q}_{i}\rightarrow 0$ to all
our results.

A similar calculation for the case of massive final states can also be
performed, which we here largely omit and refer the reader to Appendix A. The
most relevant modification when introducing masses, while still working in the
high energy limit, is the appearance of a mass dependent term where the spins
of the quark ($r$) and anti-quark ($s$) can be anti-parallel in the vertex
(see Eq. (51)). An additional phase factor proportional to $m^{2}$ (cf. Eq.
(43)) is also generated by the mass term. The polarization averaged result is
given by:

$\begin{split}&\frac{dN^{rs}}{dzd\phi}=\frac{\alpha_{s}}{(2\pi)^{2}}\text{Re}\,\int_{X_{v}^{+}}\frac{\sqrt{c_{1x}}\sqrt{c_{1y}}}{\sqrt{c_{3x}}\sqrt{c_{3y}}}e^{-i\frac{m^{2}}{2\tilde{q}_{0}^{+}}\Delta
t}\\\
&\times\Bigg{\\{}\delta_{rs}\left(P_{qg}^{vac}(z)\left[\frac{(c_{1y}c_{2y}c_{3x}-c_{1x}c_{2x}c_{3y})(c_{3y}\cos^{2}{\phi}-c_{3x}\sin^{2}{\phi})}{\left(c_{3y}\cos^{2}{\phi}+c_{3x}\sin^{2}{\phi}\right)^{2}}\right.\right.+\\\
&\left.\left.+\frac{c_{3x}c_{3y}c_{4}}{c_{3y}\cos^{2}{\phi}+c_{3x}\sin^{2}{\phi}}\right]+ir(1-2z)\frac{c_{3x}c_{3y}c_{5}\sin{\phi}\cos{\phi}}{\left(c_{3y}\cos^{2}{\phi}+c_{3x}\sin^{2}{\phi}\right)^{2}}\right)\\\
&\delta_{r,-s}\left(\frac{m^{2}}{4\tilde{q}_{0}^{+}}\frac{c_{3x}c_{3y}}{c_{3y}\cos^{2}{\phi}+c_{3x}\sin^{2}{\phi}}\right)\Bigg{\\}}\,,\end{split}$
(31)

where now we use the quark’s spin-projection along $z$ rather than the
helicity $h$, since we are considering finite masses.

##  3 Jet observables

In the previous section, we provided the cross-section to produce a quark
antenna in a dense anisotropic medium, for both massless and massive final
states. Using these results, we now study the impact of the medium anisotropy,
encapsulated into $\hat{q}$, on jet observables. We will first study the total
particle distributions. This is complemented by a study of the obtained
observables in terms of their Fourier decomposition, similar to the extraction
of flow harmonics commonly performed in soft physics studies in URHICs, see
e.g. Snellings:2011sz . Finally we discuss how measurements of spin along
different directions in the transverse plane can give information about jet
polarization.

###  3.1 Azimuthal particle distribution

We first consider the azimuthal particle distribution given by Eqs. (30) and
(31). The massless case only has one possible helicity configuration
($\delta_{h^{\prime},-h}$) as we have seen and for the massive case we will
focus on the same spin configuration ($\delta^{rs}$), given by the first term
in Eq. (31). The numerical results are presented in terms of the dimensionless
quantities

$\zeta=\frac{\sqrt{\hat{q}_{y}}-\sqrt{\hat{q}_{x}}}{\sqrt{\hat{q}_{y}}+\sqrt{\hat{q}_{x}}}\,,\quad
r=L^{+}\frac{\left(\sqrt{\hat{q}_{y}}+\sqrt{\hat{q}_{x}}\right)}{2\sqrt{q_{0}^{+}}}\,,\qquad\mu=\frac{\sqrt{2}m^{2}}{\left(\sqrt{\hat{q}_{y}}+\sqrt{\hat{q}_{x}}\right)\sqrt{q_{0}^{+}}}\,,$
(32)

which give a measure of the anisotropy $\zeta$, a rescaled medium length $r$
and a rescaled quark mass $\mu$. For the numerical evaluations we estimate the
values of these parameters by considering $L^{+}/\sqrt{2}=10$ fm,
$q_{0}^{+}/\sqrt{2}=100$ GeV, and $\hat{q}_{x}=1\,{\rm GeV}^{2}{\rm fm}^{-1}$.
Writing $\hat{q}_{y}=c\hat{q}_{x}$ and using $m=4$ GeV/c2, we have that

$\displaystyle c=2:\,\zeta=0.17\,,r=3.20\,,\mu=1.76\,,$ $\displaystyle
c=20:\,\zeta=0.63\,,r=7.28\,,\mu=0.78\,.$ (33)

Thus, in what follows we take $0.10\leq\zeta\leq 0.95$, $r=5$, and
$\mu=0.1,0.4$.888To better illustrate the anisotropic effects, in the
numerical study we used slightly smaller values for $\mu$, since the mass
suppresses the medium signature. This can be cleanly seen in the results
shown.

In Fig. 2 we show the evolution of the particle distribution for massless
quark as a function of $\phi$ for several values of the anisotropic
coefficient $\zeta$. The result for massive quarks is shown in Fig. 2. The
left columns shows the ratio between the $h=1$ distribution and the isotropic
one ($\zeta=0$). The isotropic distribution is flat ($\phi$ independent), so
the ratio is merely a change in normalization of all curves. In both massless
and massive cases we see that a larger anisotropy results in a larger ratio
near $\pi/2$, where the distinction with respect to the isotropic case seems
the largest. Additionally, there is a small suppression near $0$ and $\pi$
which only slightly increases with the value of the anisotropy. For
$\hat{q}_{x}>\hat{q}_{y}$, i.e. $\zeta<0$, we would obtain the reverse
behaviour – peaks around $0$ and $\pi$. We further observe that the
distribution is approximately symmetric with respect to $\phi=\pi/2$, which
indicates that the helicity independent terms (symmetric for
$\phi\rightarrow\pi-\phi$) dominate over the dependent ones (anti-symmetric
for the same transformation). This results in an elliptic-type deformation of
the particle distribution in the transverse plane. On the right columns we
show the ratio between the $h=\pm 1$ distributions, which quantifies the
importance of the helicity dependent term. One can distinguish the two
helicity states better the greater the anisotropy and the softest the
splitting (owing to the $(1-2z)$ factor) and this is again more pronounced
around $\pi/2$. The symmetric behaviour is purely a consequence of two
distributions corresponding to $h=\pm 1$ being mirrored with respect to
$\phi=\pi/2$. The two peaks being in distinctly different regions indicates
that positive and negative helicity states will populate different regions
inside the jet. Regarding the massive case in specific, it seems that
including mass reduces the sensitivity to anisotropy values (Fig. 3 on the
left) and that the distinction between the two spin states becomes much less
pronounced (Fig. 3 on the right).

Figure 2: Particle number distribution (Eq. (30)) as a function of $\phi$ for
different values of the anisotropy parameter $\zeta$ and for fixed $r=5$.
Left: ratio of the $h=1$ case with respect to the isotropic ($\zeta=0$) case.
Right: ratio of the $h=1$ case with respect to the $h=-1$ case. The top row
corresponds to $z=0.1$ and bottom to $z=0.4$. The vacuum piece is subtracted.

Figure 3: Particle number distribution (Eq. (31)) as a function of $\phi$ for
different values of the anisotropy parameter $\zeta$ and for fixed $r=5$ and
$z=0.1$. Note that the spin index $r$ is renamed to $h$. Left: ratio of the
$h=1$ case with respect to the isotropic ($\zeta=0$) case. Right: ratio of the
$h=1$ case with respect to the $h=-1$ case. The top row corresponds to
$\mu=0.1$ and bottom to $\mu=0.4$. The vacuum piece is subtracted.

###  3.2 Fourier decomposition

The different symmetry properties of the two contributions, spin/helicity
dependent and independent, to the final state distribution suggests that one
can isolate the helicity terms through a harmonic Fourier decomposition. We
can further motivate this decomposition by performing an expansion of the
distribution in powers of $\zeta$. We find at leading order

$\begin{split}&\frac{dN^{h}}{dzd\phi}=\frac{\alpha_{s}}{2\pi}\text{Re}\,\int_{X_{v}^{+}}e^{-i\frac{m^{2}}{2\tilde{q}_{0}^{+}}\Delta
t}\Bigg{[}2c_{1}^{2}P_{qg}^{vac}(z)+\frac{2c_{1}^{2}}{c_{3}}\,\zeta\Bigg{(}ih\frac{1-2z}{2}\bar{\Omega}^{2}\Delta
t\sin(2\phi)+\\\ &P_{qg}^{vac}(z)\left(\frac{4c_{1}^{2}(c_{2}+c_{3})\Delta
t-c_{2}^{2}(1+c_{3}\Delta t)}{c_{3}}-(c_{2}+2c_{3}-4c_{1}^{2}\Delta
t)\right)\cos(2\phi)\Bigg{)}\Bigg{]}\end{split}$ (34)

where we used a new shorthand notation

$\begin{split}&c_{1}=\frac{\bar{\Omega}}{2i\sin{\bar{\Omega}\Delta
t}}\,,\quad\qquad c_{2}=\frac{\bar{\Omega}}{\tan{\bar{\Omega}\Delta
t}}\,,\quad
c_{3}=\frac{\Delta_{L}^{+}\bar{\Omega}^{2}}{2}-c_{2}\,,\end{split}$ (35)

where $\bar{\Omega}=\left(\frac{1-i}{\sqrt{2}}\right)\bar{\omega}$, and
$\bar{\omega}=(\omega_{y}+\omega_{x})/2$. Equation (34) neatly illustrates
that, at leading power in $\zeta$, the helicity independent terms are
proportional to $\cos(2\phi)$, i.e. the elliptic contribution, while the
helicity dependence goes with $\sin(2\phi)$, as argued above. In fact, this
observation can be promoted to all orders in $\zeta$. Performing the Fourier
decomposition for finite $\zeta$ then yields

$\displaystyle\frac{2\pi}{dN^{h}/dz}\frac{dN^{h}}{dzd\phi}=1+\sum_{n=1}^{\infty}v^{(h)}_{n}\cos\left(n\phi\right)+\sum_{n=1}^{\infty}w^{(h)}_{n}\sin(n\phi)\,,$
(36)

where $v_{n}\equiv(dv_{n}^{(h)}/dz)/(dv_{0}^{(h)}/dz)$,
$w_{n}=(dw_{n}^{(h)}/dz)/(dv_{0}^{(h)}/dz)$, and we have the relations

$\displaystyle v_{0}^{(h)}$
$\displaystyle=\frac{1}{2\pi}\int_{0}^{2\pi}d\phi\,\,2\pi\frac{dN^{h}}{dzd\phi}\,,$
$\displaystyle v_{n}^{(h)}$
$\displaystyle=\frac{1}{\pi}\int_{0}^{2\pi}d\phi\,\,2\pi\frac{dN^{h}}{dzd\phi}\,\cos\left(n\phi\right)\,,$
$\displaystyle w_{n}^{(h)}$
$\displaystyle=\frac{1}{\pi}\int_{0}^{2\pi}d\phi\,\,2\pi\frac{dN^{h}}{dzd\phi}\,\sin\left(n\phi\right)\,.$
(37)

At this point, it is important to note that several of these Fourier
coefficients can be immediately determined from the symmetries of the spectrum
in Eq. (30). First of all, in the case of an isotropic medium, all Fourier
coefficients except for $v_{0}$ vanish. For an anisotropic medium, both
$v^{(h)}_{n}$ and $w_{n}^{(h)}$ vanish for odd $n$ due to discrete symmetry
$\phi\rightarrow\phi-\pi$ of the particle distribution. This is expected since
the presence of the anisotropic $\hat{q}$ can not lead to any direct flow-like
contribution. The coefficients $v_{n}$ for even $n$ are non-vanishing but
independent of helicity, i.e., $v^{(h)}_{n}=v_{n}$, a consequence of
$\phi\rightarrow-\phi$ antisymmetry of the helicity dependent term. The
helicity dependence can only enter via the odd-parity terms $w_{n}^{(h)}$ in
the series. This is analogous to charge imbalance measurements in chiral
magnetic effect studies, see e.g. Kharzeev:2015znc .999Although we note that
in our context we do not need to look into fluctuations, as is the case in
studies of the chiral magnetic effect in heavy ions. It is also easy to
realize that $w_{n}^{(+1)}=-w_{n}^{(-1)}$, since $w_{n}^{(h)}\propto h$. Thus,
the presence of non-zero values of any of these coefficients gives information
about the underlying properties of the medium and its coupling to the jet in
different ways.

In Figs. 4 we show the evolution of the leading harmonics, $v_{2}$ and
$w_{2}^{(+1)}$, with $\zeta$ (top) and $z$ (bottom). An equivalent plot for
massive quarks is shown in Fig. 5. First, when comparing the $v_{2}$ with
$w_{2}^{(+1)}$ harmonics in absolute value, we get a possible explanation for
the near symmetry around $\pi/2$ on the left column plots of Fig. 2 and 3,
since this indicates that $\cos(2\phi)$ dominates. Then, $v_{2}$ is larger for
more imbalanced splittings and grows with the absolute value of the anisotropy
parameter, as one should expect. In contrast, the evolution of the $w_{2}$
harmonic is non-monotonic with the energy fraction; this is a consequence of
the fact that for $z=0,\,0.5$ there is no helicity/spin dependence and thus
the harmonic has to vanish as shown. As a result, the curves display a maximum
for intermediates values of $z$. The evolution in $\zeta$ is also surprising
since at large values, the curves tends to flatten and some even display a
maximum. These observations are even more striking in the massive case, where
the flow coefficient can become negative and it exhibits a more complex
evolution in $z$. We note that this behavior can be tied to the larger value
of $\mu=0.4$ considered, where some of the assumptions used can break down.
Additionally, we see quite smaller values of $w_{2}$ for the massive case,
explaining the less pronounced distinction between spin states on the right
column plots of Fig. 3. Overall, mass seems to confound the interesting
features of the azimuthal particle distribution.

To summarize, extracting the flow harmonics inside jets would provide a new
window into the jet substructure, with a clear sign for non-trivial medium
properties being isolated from more traditional medium effects, e.g. momentum
broadening, energy loss, which are mainly isotropic and thus can generate any
harmonic but $v_{0}$.

Figure 4: Left: Evolution of $v_{2}$ in Eq. (36) as a function of $\zeta$ (top
line) or of $z$ (bottom line) and for different values of $z$ or of $\zeta$,
respectively, with fixed $r=5$. Right: equivalent plots for $w_{2}^{(+1)}$ in
Eq. (36). The vacuum piece is subtracted.

Figure 5: Left: Evolution of $v_{2}$ in Eq. (36) as a function of $\zeta$ (top
line) or of $z$ (bottom line) and for different values of $z$ or of $\zeta$,
respectively, with fixed $r=5$ and $\mu=0.4$. Right: equivalent plots for
$w_{2}^{(+1)}$ in Eq. (36). The vacuum piece is subtracted.

###  3.3 Transverse polarization for a massive antenna

So far we have defined the quark and anti-quark spin states as the projection
along the $z$ axis. However, because we are trying to tap the interplay
between spin and the anisotropy the $q\bar{q}$ antenna experiences in the
transverse plane, it might be more interesting to consider spin projected
along some axis on that plane. One can do so by changing the expression for
the vertex following the steps detailed in Appendix A and inserting it in Eq.
(18). Carrying out the remaining of the calculation exactly as before, the
particle distribution reads

$\displaystyle\frac{dN^{rs}}{dzd\phi}=\frac{\alpha_{s}}{4(2\pi)^{2}}\text{Re}\,\int_{X_{v}^{+}}e^{-i\frac{m^{2}}{2\tilde{q}_{0}^{+}}\Delta
t}\frac{\sqrt{c_{1x}}\sqrt{c_{1y}}}{\sqrt{c_{3x}}\sqrt{c_{3y}}}\delta_{rs}$
$\displaystyle\left(\frac{(c_{1y}c_{2y}c_{3x}-c_{1x}c_{2x}c_{3y})(c_{3y}\cos^{2}{\phi}-c_{3x}\sin^{2}{\phi})}{\left(c_{3y}\cos^{2}{\phi}+c_{3x}\sin^{2}{\phi}\right)^{2}}+\left(\frac{m^{2}}{2\tilde{q}_{0}^{+}}+c_{4}\right)\frac{c_{3x}c_{3y}}{c_{3y}\cos^{2}{\phi}+c_{3x}\sin^{2}{\phi}}+\right.$
$\displaystyle\left.-\frac{(1-i)mr\sqrt{\pi}}{2\sqrt{\tilde{q}_{0}^{+}}}\sqrt{\frac{c_{3x}c_{3y}}{c_{3y}\cos^{2}\phi+c_{3x}\sin^{2}\phi}}\frac{(2ic_{1x}-c_{2x})c_{3y}\cos\phi\sin\alpha-(2ic_{1y}-c_{2y})c_{3x}\sin\phi\cos\alpha}{\left(c_{3y}\cos^{2}{\phi}+c_{3x}\sin^{2}{\phi}\right)}\right)\,,$
(38)

where $\alpha$ is the angle that the spin quantization axis makes with the
$x$-axis, i.e., $\alpha=0,\pi/2$ corresponds to measuring spin along the $x$
and $y$ axis, respectively. As discussed in Appendix A and as is clear from
the previous expression, when there is a very large anisotropy,
$\hat{q}_{y}\gg\hat{q}_{x}$ ($\zeta\rightarrow 1$), the outgoing $q\bar{q}$
pair is unpolarized along the $y$ direction ($\alpha=\pi/2$). However, when
measuring spin along the $x$-axis ($\alpha=0$), the resulting distribution is
polarized. Notice that this effect is solely driven by the mass dependent
term. This observation can be confirmed in Fig. 6 and 7, where we plot Eq.
(3.3) for several values of $\zeta$, $z$ and $\mu$ as a function of $\phi$.
Indeed, for very large anisotropies like $\zeta=0.95$ and $0.7$, one can
clearly observe that when $\alpha=\pi/2$ (right plots) the solid and dashed
curves are quite similar to each other, reflecting a result that is nearly
independent of $h$ and thus a $q\bar{q}$ that is very weakly polarized along
the $y$-axis. For $\alpha=0$, the two curves are significantly different
around $\phi=\pi/2$, giving rise to a significant polarization along the
$x$-axis. Even for intermediate values of anisotropy, for instance $\zeta=0.4$
(the green curve), the $h$ dependence, i.e., degree of polarization, is quite
clearly larger along the $x$-axis than along the $y$-axis. Regarding the
dependence in $z$, we see that the softer splittings ($z=0.1$ in Fig. 6 ) seem
to give rise to a more pronounced separation between distributions of
different spin states, enhancing the polarization effect. The mass $\mu$,
despite being the parameter enabling this interesting effect, seems to
suppress the degree of polarization when taking larger values (bottom plots of
Fig. 6 and 7). Additionally, one should note that in the isotropic case
(violet curve in all plots) the $q\bar{q}$ pair also obtains polarization due
to interactions with the medium. Nevertheless, the stark difference with
respect to the anisotropic case is that the degree of polarization is the
exact same in all directions $\alpha$. This is clear if one plots the
distribution as a function of $\phi-\alpha$ since, in the isotropic limit of
Eq. (3.3), the spin-dependent term is simply proportional to
$\cos(\phi-\alpha)$. Thus, for anisotropic matter, this direction dependent
polarization effect gives yet another feature which can be explored to extract
information from the medium, since it informs about the direction of maximal
$\hat{q}$.

Figure 6: Plots of the distribution in Eq. (3.3) as a function of $\phi$ for
different values of the anisotropy parameter $\zeta$, for $\alpha=0$ (left)
and $\alpha=\pi/2$ (right), for fixed $r=5$, $z=0.1$ and for two different
values of $\mu$ – $0.1$ (upper) and $0.4$ (lower). Note that the spin index
$r$ is renamed to $h$. The vacuum piece is subtracted.

Figure 7: Plots of the distribution in Eq. (3.3) as a function of $\phi$ for
different values of the anisotropy parameter $\zeta$, for $\alpha=0$ (left)
and $\alpha=\pi/2$ (right), for fixed $r=5$, $z=0.4$ and for two different
values of $\mu$ – $0.1$ (upper) and $0.4$ (lower). Note that the spin index
$r$ is renamed to $h$. The vacuum piece is subtracted.

##  4 Conclusions

In this work we have studied $q\bar{q}$ pair production from a gluon in the
presence of an anisotropic background QCD medium, as expected to exist in
several of the stages characterising the aftermath of URHICs, in particular in
the very initial stages, possibly dominated by strong longitudinal color
fields. We consider a simple medium model, where the background is static and
confined inside a box of fixed longitudinal size. We introduce anisotropic
effects by allowing the jet quenching parameter to differ in two orthogonal
directions in the transverse plane with respect to the jet axis. This leads to
the final particle distribution evolving differently in each of these
directions. Furthermore, we keep track of the helicity/spin dependence of the
final states, which we show directly couples to the medium anisotropy,
offering a new window to study directional effects inside jets in heavy ions
collisions. This possibility is very appealing as the direction of the
anisotropy is easy to identify experimentally. The fact that spin flip is
energy suppressed during the propagation of the parton indicates that these
spin-effects from the very early times could survive the subsequent evolution
of the medium.

Our final main results, given in Eqs. (2) and (31), can be straightforwardly
extended beyond the approximations we consider here, although we do not expect
the main physical picture portrayed to change qualitatively. Further, we
discuss their potential impact in several jet observables, indicating when
clear signatures of the medium anisotropy can be found in the jet
distribution. As examples of these, we propose the study of Fourier azimuthal
harmonics inside jets as a way to probe non-trivial correlations. Lastly, we
argue that spin measurements in different directions allow to probe the
polarization of partons inside the jet, further testing directional effects
induced by the medium. We leave a more detailed phenomenological study of such
observables for future work.

Finally, we note that the presented theoretical results can, in principle, be
extended to other hard probes, such as quarkonia or heavy flavour. We believe
that this would allow to expand the current toolkit of tomographic tools of
the matter produced in URHICs, allowing for a more complete energy scan, which
is not possible by just using high energy jets, see e.g. Kumar:2022ylt ;
Du:2023izb .

###### Acknowledgements.

We wish to thank A. V. Sadofyev and R. Szafron for helpful discussions. We are
especially thankful to A. V. Sadofyev for comments concerning the Fourier
decomposition of the particle distribution. This work is supported by the
European Research Council project ERC-2018-ADG-835105 YoctoLHC; by Maria de
Maeztu excellence unit grant CEX2023-001318-M and project PID2020-119632GB-I00
funded by MICIU/AEI/10.13039/501100011033; by ERDF/EU; and by Fundação para a
Ciência e a Tecnologia (FCT, Portugal) under project CERN/FIS-PAR/0032/2021.
It has received funding from Xunta de Galicia (CIGUS Network of Research
Centres); and from the European Union’s Horizon 2020 research and innovation
programme under grant agreement No. 824093. The Feynman diagrams in this work
were produced using the JaxoDraw software Binosi:2003yf . JB is supported by
the United States Department of Energy under Grant Contract DESC0012704. JMS
is supported by Fundação para a Ciência e Tecnologia (Portugal, FCT) under
contract PRT/BD/152262/2021.

##  Appendix A Effective in-medium Feynman rules in coordinate space

In this appendix we briefly review some elements entering the main text’s
calculation. Further details can be found in e.g. Mehtar-Tani:2013pia ;
Blaizot:2012fh .

###  A.1 Parton scalar propagators

The quark propagator in the presence of a classical background gauge field is
given at leading order in powers of $1/p^{+}$ by

$\mathcal{G}^{ij}(y^{+},\boldsymbol{y};x^{+},\boldsymbol{x}|p^{+})=\int^{\boldsymbol{r}(y^{+})=\boldsymbol{y}}_{\boldsymbol{r}(x^{+})=\boldsymbol{x}}\mathcal{D}\boldsymbol{r}(\xi)\exp{\left\\{i\frac{p^{+}}{2}\int_{x^{+}}^{y^{+}}ds^{+}\,\dot{\boldsymbol{r}}^{2}\right\\}}U^{ij}(x^{+},y^{+},\boldsymbol{r}(\xi))\,,$
(39)

where the fundamental Wilson line is given by

$U^{ij}(x^{+},y^{+},\boldsymbol{r})=\mathcal{P}\exp{\left\\{ig\int_{x^{+}}^{y^{+}}ds^{+}\,\mathcal{A}_{a}^{-}(s^{+},\boldsymbol{r}(s^{+}))t^{a}_{ij}\right\\}}\,.$
(40)

The anti-quark propagator $\overline{\mathcal{G}}$ is the same as the quark
propagator with the replacement $U\rightarrow U^{\dagger}$. As for the gluon,
one only needs to replace the fundamental Wilson line by an adjoint one
$U_{A}$,

$U_{A}^{ab}(x^{+},y^{+},\boldsymbol{r})=\mathcal{P}\exp{\left\\{g\int_{x^{+}}^{y^{+}}ds^{+}\,\mathcal{A}_{c}^{-}(s^{+},\boldsymbol{r}(s^{+}))f^{cab}\right\\}}\,.$
(41)

In practice, one can combine these effective propagators with standard Feynman
diagrammatic techniques with the following modifications compared to the
canonical rules for QCD:

$\begin{split}\text{{internal
line}}\qquad\rightarrow\qquad&\mathcal{G}^{ij}(y^{+},\boldsymbol{y};x^{+},\boldsymbol{x}|p^{+})/2p^{+}\\\\[8.61108pt]
\text{{incoming external
line}}\qquad\rightarrow\qquad&\mathcal{G}^{ij}(y^{+},\boldsymbol{y};x^{+},\boldsymbol{x}|p^{+})e^{-i(\frac{\boldsymbol{p}^{2}}{2p^{+}}x^{+}-\boldsymbol{p}\cdot\boldsymbol{x})}\\\\[8.61108pt]
\text{{outgoing external
line}}\qquad\rightarrow\qquad&\mathcal{G}^{ij}(y^{+},\boldsymbol{y};x^{+},\boldsymbol{x}|p^{+})e^{i(\frac{\boldsymbol{p}^{2}}{2p^{+}}y^{+}-\boldsymbol{p}\cdot\boldsymbol{y})}\\\
\end{split}$ (42)

For massive quarks, one only needs to account for the additional phases coming
from the intermediate integrations over $p_{n}^{-}$ in between each medium
scattering, all of which cancel except those at the light-cone time limits,
i.e., one makes the replacement:

$\mathcal{G}^{ij}(y^{+},\boldsymbol{y};x^{+},\boldsymbol{x}|p^{+})_{m\neq
0}\rightarrow\mathcal{G}^{ij}(y^{+},\boldsymbol{y};x^{+},\boldsymbol{x}|p^{+})_{m=0}\exp{\left\\{-i\frac{m^{2}}{2p^{+}}(y^{+}-x^{+})\right\\}}\,,$
(43)

where $m$ is the quark mass.

###  A.2 Spin/helicity-dependent quark-gluon vertex

We assume the quark-gluon vertex is unchanged by the presence of the medium
and we write it as

$(A^{\lambda rs})_{kl}^{b}(z,\boldsymbol{P})=igt_{kl}^{b}V^{\lambda
rs}(z,\boldsymbol{P})\,,$ (44)

where the relative transverse momentum reads
$\boldsymbol{P}=\boldsymbol{p}_{q}-z\boldsymbol{p}_{g}$, with $z$ the fraction
of $+$ momentum from the gluon carried by the quark, $\boldsymbol{p}_{g}$ the
former’s transverse momentum and $\boldsymbol{p}_{q}$ the latter’s. At leading
order in $1/p^{+}$, the only transverse momentum dependence of the vertex
comes from a linear term in ${\boldsymbol{P}}$, as we will see below. The
labels $\lambda,r,s$ represent, respectively, the gluon polarization state in
circular basis and the $q$ and $\bar{q}$ spin projection along a general spin
quantization axis or the helicity states.
To insert a momentum-space vertex via coordinate space rules, one replaces
transverse momenta by derivatives of $\delta$-functions. This can be explained
by going from the mixed space (transverse momentum + light-cone time)
formulation of the amplitude to the coordinate space one. The amplitude in
Fig. 1 would have the following form in mixed space:

$i\mathcal{M}=\int_{\boldsymbol{p}_{q},...}\boldsymbol{p}_{q}\text{
}\mathcal{G}^{ik}(L^{+},\boldsymbol{p_{1}};x_{v}^{+},\boldsymbol{p_{q}}|p_{1}^{+})R^{kj}(\boldsymbol{p}_{q},...)$
(45)

where $R^{kj}(\boldsymbol{p}_{q},...)$ represents the remaining part of the
amplitude. Let us focus only on the vertex term proportional to
${\boldsymbol{p}}_{q}$ (the other term in ${\boldsymbol{P}}$ is completely
analogous). We would like to write the vertex in such a way that one could
directly extract the coordinate space Feynman rule for it. To do so, we just
need to use the convolution theorem for Fourier transforms:

$\begin{split}i\mathcal{M}&=\int_{\boldsymbol{p}_{q},...}\mathcal{F}\left(\mathcal{F}^{-1}(\boldsymbol{p}_{q})\ast\mathcal{F}^{-1}(\mathcal{G}^{ik}(L^{+},\boldsymbol{p_{1}};x_{v}^{+},\boldsymbol{p_{q}}|p_{1}^{+}))\right)(\boldsymbol{p}_{q})R^{kj}(\boldsymbol{p}_{q},...)\\\
&=\int_{\boldsymbol{p}_{q},...}\int_{\boldsymbol{x_{v}}}e^{i\boldsymbol{x_{v}}\cdot\boldsymbol{p}_{q}}\int_{\boldsymbol{y}}i\delta^{{}^{\prime}}(\boldsymbol{y})\mathcal{G}^{ik}(L^{+},\boldsymbol{\boldsymbol{p_{1}}};x_{v}^{+},\boldsymbol{\boldsymbol{x_{v}-y}}|p_{1}^{+})R^{kj}(\boldsymbol{p}_{q},...)\\\
&=\int_{\boldsymbol{x_{v}},...}\int_{\boldsymbol{y}}i\delta^{{}^{\prime}}(\boldsymbol{y})\mathcal{G}^{ik}(L^{+},\boldsymbol{\boldsymbol{p_{1}}};x_{v}^{+},\boldsymbol{\boldsymbol{x_{v}-y}}|p_{1}^{+})R^{kj}(\boldsymbol{x_{v}},...)\,.\end{split}$
(46)

Hence, to write the vertex directly in coordinate space one should replace the
momentum coming out of the vertex as $\boldsymbol{p}\rightarrow
i\delta^{{}^{\prime}}(\boldsymbol{y})$ and displace the vertex coordinate
associated to the corresponding scalar propagator by $-\boldsymbol{y}$,
integrating over this displacement. For momenta coming into the vertex the
same reasoning holds but the correct momentum replacement is
$\boldsymbol{q}\rightarrow-i\delta^{{}^{\prime}}(\boldsymbol{y})$.
Accordingly, to write the vertex in Eq. (44) in coordinate space, one attaches
coordinate displacement ${\boldsymbol{x}}$ to the incoming gluon and
${\boldsymbol{y}}$ to the outgoing quark and makes the substitution:

${\boldsymbol{P}}\rightarrow
i(\delta({\boldsymbol{x}})\delta^{\prime}({\boldsymbol{y}})+z\delta^{\prime}({\boldsymbol{x}})\delta({\boldsymbol{y}}))\,.$
(47)

Further below, when presenting the squared expression for the vertex, we will
use $\overline{{\boldsymbol{P}}}$ for the relative momentum in the complex
conjugate amplitude. One should keep in mind the correct substitution for this
momentum is the complex conjugate of Eq. (47) with different (dummy)
coordinate displacements.

The helicity-dependent vertex for circularly polarized gluons and for massless
quarks is given by

$\begin{split}&V^{\lambda hh^{\prime}}(z,\boldsymbol{P})=2\gamma^{\lambda
h}(z)\delta_{h,-h^{\prime}}(\boldsymbol{\epsilon}^{\lambda}\cdot\boldsymbol{P})\\\
&\gamma^{\lambda h}(z)=\frac{1}{\sqrt{z(1-z)}}\left(z\delta^{\lambda
h}-(1-z)\delta^{\lambda-h}\right)\,.\end{split}$ (48)

As stated in the main text, the quark’s helicity is given by $h=\pm 1$, the
anti-quark’s by $-h$ and the gluon polarization vector by
$\boldsymbol{\epsilon}^{\lambda}=\frac{1}{\sqrt{2}}(1,\lambda i)$, i.e, the
circular basis with $\lambda=\pm 1$. When squaring and averaging over
polarizations one gets:

$\begin{split}&\frac{1}{2}\sum_{\lambda}(V^{\lambda
hh^{\prime}})(z,\boldsymbol{P})\left((V^{\lambda
hh^{\prime}})(z,\boldsymbol{\overline{P}})\right)^{*}=\\\
&\delta_{h,-h^{\prime}}\frac{2}{z(1-z)}\left(P_{qg}^{vac}(z)(\boldsymbol{P}\cdot\boldsymbol{\overline{P}})+\frac{ih}{2}(1-2z)({\boldsymbol{P}}\times\overline{{\boldsymbol{P}}})_{z}\right)\,,\end{split}$
(49)

where the $g\rightarrow q\bar{q}$ vacuum splitting function reads
$P_{qg}^{vac}(z)=T_{R}(z^{2}+(1-z)^{2})$ and
$({\boldsymbol{P}}\times\overline{{\boldsymbol{P}}})_{z}={\boldsymbol{P}}_{x}\overline{{\boldsymbol{P}}}_{y}-{\boldsymbol{P}}_{y}\overline{{\boldsymbol{P}}}_{x}$.

For massive quarks, a (non-suppressed) term where the quark and anti-quark
spins can be anti-parallel arises , i.e.,

$V^{\lambda rs}(z,\boldsymbol{P})=2r\,\gamma^{\lambda
r}(z)\delta_{rs}(\boldsymbol{\epsilon}^{\lambda}\cdot\boldsymbol{P})-\delta^{\lambda
r}\delta_{r,-s}\frac{2m}{\sqrt{2z(1-z)}}\,,$ (50)

where now we work with spin states $r,s$ instead of helicity states since the
latter ones are not Lorentz invariant for massive quarks. The parallel spins
term also contains mass dependent terms, but these are suppressed by an energy
denominator with respect to the remaining. When squaring and averaging over
polarizations one gets:

$\begin{split}&\frac{1}{2}\sum_{\lambda}(V^{\lambda
rs})(z,\boldsymbol{P})\left((V^{\lambda
rs})(z,\boldsymbol{\overline{P}})\right)^{*}=\\\
&\frac{2}{z(1-z)}\left\\{\delta_{rs}\left(P_{qg}^{vac}(z)(\boldsymbol{P}\cdot\boldsymbol{\overline{P}})+\frac{ir}{2}(1-2z)({\boldsymbol{P}}\times\overline{{\boldsymbol{P}}})_{z}\right)+\delta_{r,-s}\frac{m^{2}}{2}\right\\}\,.\end{split}$
(51)

It can also be interesting to measure the spins of the $q\bar{q}$ pair along
some particular axis, in which case one can generalize the previously obtained
results to a general spin-quantization axis
$\vec{n}=\left(\sin\theta\cos\alpha,\sin\theta\sin\alpha,\cos\theta\right)$
($\theta=0,\alpha=\pi$ gives us spin quantized along the $z$ axis, as in Eq.
(50)):

$\begin{split}&V^{\lambda rs}(z,\boldsymbol{P})=\frac{2}{\sqrt{z(1-z)}}\\\
&\left\\{\left.r\delta_{rs}\left(\left[\cos\theta\left(z\delta^{\lambda
r}-(1-z)\delta^{\lambda,-r}\right)+\sin^{2}(\theta/2)\left(\delta^{\lambda
r}-\delta^{\lambda,-r}\right)\right](\boldsymbol{\epsilon}^{\lambda}\cdot\boldsymbol{P})-\frac{m}{2\sqrt{2}}\sin\theta
e^{i\lambda\alpha}\right)+\right.\right.\\\
&\left.+\delta_{r,-s}\left(\sin\theta(z-1/2)(\boldsymbol{\epsilon}^{\lambda}\cdot\boldsymbol{P})+\frac{m}{\sqrt{2}}e^{i\lambda\alpha}\left[\delta^{\lambda
r}\cos^{2}(\theta/2)-\delta^{\lambda,-r}\sin^{2}(\theta/2)\right]\right)\right\\}\,.\end{split}$
(52)

Furthermore, if we consider $\theta\neq 0$, i.e, a spin-quantization axis
different from $z$, cross terms which are linear in the momenta (derivatives
in the coordinate space amplitude) and mass appear when squaring the vertex.
For instance, if we set $\theta=\pi/2$, i.e., measuring the spin of the quark
and anti-quark along some axis in the transverse plane, square Eq. (52) and
average over polarization states we get:

$\begin{split}&\frac{1}{2}\sum_{\lambda}(V^{\lambda
rs})(z,\boldsymbol{P})\left((V^{\lambda
rs})(z,\boldsymbol{\overline{P}})\right)^{*}=\\\
&\frac{1}{2z(1-z)}\left\\{(\boldsymbol{P}\cdot\boldsymbol{\overline{P}})\left(\delta_{rs}+(2z-1)^{2}\delta_{r,-s}\right)+m^{2}\right.\\\
&\left.+imr\left[\left(\boldsymbol{P}^{x}-\boldsymbol{\overline{P}}^{x}\right)\sin\alpha-\left(\boldsymbol{P}^{y}-\boldsymbol{\overline{P}}^{y}\right)\cos\alpha\right]\left(\delta_{rs}-(2z-1)\delta_{r,-s}\right)\right\\}\,.\end{split}$
(53)

The linear terms vanish in the vacuum limit since the momenta in the direct
and complex conjugate amplitudes are the same. This means that the difference
${\boldsymbol{P}}^{i}-\bar{\boldsymbol{P}}^{i}$ is controlled by $\hat{q}_{i}$
in the corresponding direction. If $\hat{q}^{y}\gg\hat{q}^{x}$, for instance,
then the spin dependent term is dominated by the $\cos\alpha$ term. This
implies that when measuring spin states along the $x$ axis ($\alpha=0$) the
difference between the spectrum for up and down ($r=\pm 1$) is much greater
than the difference obtained when measuring spin along the $y$ axis
($\alpha=\pi/2$). Put simply, in the limit where $\hat{q}^{x}\rightarrow 0$
and $\hat{q}^{y}$ is finite, then the $q\overline{q}$ pair is completely
unpolarized in the $y$ direction, since there is no spin dependence, and
polarized in the $x$ direction. The degree of polarization depends on the
quark mass and naturally grows with $\hat{q}^{y}$.

##  Appendix B Calculation details

In this appendix we provide further details on the treatment of the color
structure and medium averaging underlying the main calculation.

###  B.1 Color structure

Let us explicitly write the Wilson line correlators appearing in Eq. (16)
without prefactors and path integrals (to each $r_{i}$ corresponds a different
path integral and the numbering is irrelevant):

$\begin{split}&T^{\text{color}}=\frac{1}{N_{c}^{2}-1}\left\langle
U_{A}^{ba}({\boldsymbol{r}}_{5})(U_{A}^{\dagger})^{a\overline{b}_{1}}({\boldsymbol{r}}_{9})\right\rangle\\\
&\times
t_{kl}^{b}t_{\overline{l}\overline{k}}^{\overline{b}}\left\langle(U_{A}^{\dagger})^{\overline{b}_{1}\overline{b}}({\boldsymbol{r}}_{8})(U^{\dagger)})^{lj_{1}}({\boldsymbol{r}}_{4})U^{i_{1}k}({\boldsymbol{r}}_{3})\right\rangle\\\
&\times\left\langle
U^{ii_{1}}({\boldsymbol{r}}_{1})(U^{\dagger})^{\overline{k}i}({\boldsymbol{r}}_{6})(U^{\dagger})^{j_{1}j}({\boldsymbol{r}}_{2})U^{j\overline{l}}({\boldsymbol{r}}_{7})\right\rangle\,.\\\
\end{split}$ (54)

First note that this object is in a color singlet state. This is because we
are considering the final-state quark and anti-quark to have definite color
states which are then summed over when squaring the amplitude. Since color is
conserved throughout the process (the system as a whole does not exchange
color with the medium due to the diagonal color matrix in Eq. (4)), this means
the system is in a singlet state at all times. In particular, after expanding
the color structure inside each time region as a linear combination of the
state projectors in each of its irreducible representations, one can drop the
color non-singlet part (one can include the non-singlet part in the expansion,
but the contraction with the final region will always only pick out the
singlet parts). Finding out the color singlet states amounts to finding out
states for which the second-order Casimir operator in the appropriate
representation vanishes Fukushima:2007dy ; Fukushima:2017mko  . For the two-
point function the only option is the normalized projector
$\frac{\delta^{b\bar{b}_{1}}}{\sqrt{N_{c}^{2}-1}}$. Hence, we have

$\left\langle
U_{A}^{ba}(\boldsymbol{r}_{5})(U_{A}^{\dagger})^{a\bar{b}_{1}}(\boldsymbol{r}_{9})\right\rangle=\frac{\delta^{b\bar{b}_{1}}}{N_{c}^{2}-1}\left\langle\text{Tr}\left[U_{A}(\boldsymbol{r}_{5})U_{A}^{\dagger}(\boldsymbol{r}_{9})\right]\right\rangle\,.$
(55)

Propagating the $\delta^{b\bar{b}_{1}}$ to the rest of the squared amplitude,
writing the adjoint Wilson line in terms of fundamental ones and repeatedly
making use of the Fierz identity, i.e, using

$\displaystyle
U_{A}^{b\overline{b}}=2\text{Tr}\left[t^{b}U^{\dagger}t^{\overline{b}}U\right]\,,\quad
t^{b}_{kl}t^{b}{ij}=\frac{1}{2}(\delta^{kj}\delta^{li}-\frac{1}{N_{c}}\delta^{kl}\delta^{ij})\,,$
(56)

the other two regions read:

$\begin{split}&\frac{1}{2}\Bigg{(}\left\langle\left(U({\boldsymbol{r}}_{3})U^{\dagger}({\boldsymbol{r}}_{8})\right)^{i_{1}\bar{k}}\left(U({\boldsymbol{r}}_{8})U^{\dagger}({\boldsymbol{r}}_{4}))\right)^{\bar{l}j_{1}}\right\rangle_{(x_{v}^{+},\bar{x}_{v}^{+})}\left\langle\left(U^{\dagger}({\boldsymbol{r}}_{6})U({\boldsymbol{r}}_{1})\right)^{\bar{k}i_{1}}\left(U^{\dagger}({\boldsymbol{r}}_{2})U({\boldsymbol{r}}_{7}))\right)^{j_{1}\bar{l}}\right\rangle_{(\bar{x}_{v}^{+},L^{+})}\\\
&-\frac{1}{N_{c}^{2}}\left\langle\text{Tr}\left[U({\boldsymbol{r}}_{3})U^{\dagger}({\boldsymbol{r}}_{4})\right]\right\rangle_{(x_{v}^{+},\bar{x}_{v}^{+})}\left\langle\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{6})U({\boldsymbol{r}}_{1})U^{\dagger}({\boldsymbol{r}}_{2})U({\boldsymbol{r}}_{7})\right]\right\rangle_{(\bar{x}_{v}^{+},L+)}\Bigg{)}\,.\end{split}$
(57)

To fully contract the remaining color structure in the first line, one needs
to expand each time region’s contribution in a basis of color singlets, as we
have discussed. For both cases the number of such states is two and they can
be obtained through the contractions
$\delta^{i_{1}\bar{k}}\delta^{\bar{l}j_{1}}$ and
$\delta^{i_{1}j_{1}}\delta^{\bar{l}\bar{k}}$. An orthonormal basis for these
singlet projectors is straightforward to obtain: Fukushima:2007dy ;
Fukushima:2017mko :

$\begin{split}&s_{1}^{i_{1}\bar{k}\bar{l}j_{1}}=\frac{1}{N_{c}}\delta^{i_{1}\bar{k}}\delta^{\bar{l}j_{1}}\\\
&s_{2}^{i_{1}\bar{k}\bar{l}j_{1}}=\frac{1}{\sqrt{N_{c}^{2}-1}}\left(\delta^{i_{1}j_{1}}\delta^{\bar{l}\bar{k}}-\frac{1}{N_{c}}\delta^{i_{1}\bar{k}}\delta^{\bar{l}j_{1}}\right)\,.\end{split}$
(58)

After expanding each correlator in the first line of Eq. (57) as a linear
combination of the projectors $s_{1}$ and $s_{2}$, the coefficients are easily
obtained by taking advantage of orthonormality. For instance, we have

$\begin{split}&\left\langle\left(U^{\dagger}({\boldsymbol{r}}_{6})U({\boldsymbol{r}}_{1})\right)^{\bar{k}i_{1}}\left(U^{\dagger}({\boldsymbol{r}}_{2})U({\boldsymbol{r}}_{7}))\right)^{j_{1}\bar{l}}\right\rangle=\frac{1}{N_{c}}\left\langle\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{6})U({\boldsymbol{r}}_{1})\right]\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{2})U({\boldsymbol{r}}_{7})\right]\right\rangle
s_{1}^{i_{1}\bar{k}\bar{l}j_{1}}\\\
&+\frac{1}{\sqrt{N_{c}^{2}-1}}\left(\left\langle\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{6})U({\boldsymbol{r}}_{1})U^{\dagger}({\boldsymbol{r}}_{2})U({\boldsymbol{r}}_{7})\right]\right\rangle-\frac{1}{N_{c}}\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{6})U({\boldsymbol{r}}_{1})\right]\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{2})U({\boldsymbol{r}}_{7})\right]\right)s_{2}^{i_{1}\bar{k}\bar{l}j_{1}}\end{split}$
(59)

and for the remaining correlator it is completely analogous. Using these
expansions, the fully contracted result of the color part of the squared
amplitude in Eq. (54) finally simplifies to

$T^{\text{color}}=\mathcal{C}^{(2)}\mathcal{C}^{(3)}\mathcal{C}^{(4)}\,,$ (60)

where

$\begin{split}&\mathcal{C}^{(2)}=\frac{1}{N_{c}^{2}-1}\left\langle\text{Tr}\left[U_{A}(\boldsymbol{r}_{5})U_{A}^{\dagger}(\boldsymbol{r}_{9})\right]\right\rangle\\\
&\mathcal{C}^{(3)}=\frac{1}{2(N_{c}^{2}-1)}\left\langle\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{8})U({\boldsymbol{r}}_{3})\right]\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{4})U({\boldsymbol{r}}_{8})\right]-\frac{1}{N_{c}}\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{4})U({\boldsymbol{r}}_{3})\right]\right\rangle\\\
&\mathcal{C}^{(4)}=\frac{1}{N_{c}^{2}-1}\left\langle\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{6})U({\boldsymbol{r}}_{1})\right]\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{2})U({\boldsymbol{r}}_{7})\right]-\frac{1}{N_{c}}\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{6})U({\boldsymbol{r}}_{1})U^{\dagger}({\boldsymbol{r}}_{2})U({\boldsymbol{r}}_{7})\right]\right\rangle\,,\end{split}$
(61)

which resembles results in, e.g., Apolinario:2014csa ; Isaksen:2020npj ;
Dominguez:2012ad . Both $\mathcal{C}^{(2)}$ and $\mathcal{C}^{(3)}$ have exact
analytical solutions, which are presented in Eq. (64) and (77), while
$\mathcal{C}^{(4)}$ can only be obtained numerically Isaksen:2020npj . In the
large-$N_{c}$ limit, it has been argued Apolinario:2014csa ; Blaizot:2012fh
that the quadrupole (the second term in $\mathcal{C}^{(4)}$) can be
approximated by the product of two independent averages of two Wilson line
traces, an approximation which is kinematically favoured when the formation
time of the antenna is small compared to the medium’s length. Nevertheless, in
our case, this term is color suppressed by a power $N_{c}^{2}$, so we may drop
it. This simplification would not be possible in, for instance, the in-medium
$\gamma\rightarrow q\bar{q}$ splitting, where the quadrupole must be dealt
with it Dominguez:2019ges . As for the first term in $\mathcal{C}^{(4)}$, it
was shown in Isaksen:2020npj  that in the large-$N_{c}$ limit the differential
equation for the average of two traces (first term in $\mathcal{C}^{(4)}$)
decouples from the equation for the average of a single trace (second term in
$\mathcal{C}^{(4)}$). This allows the former to be calculated exactly in this
limit:

$\displaystyle\left\langle\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{6})U({\boldsymbol{r}}_{1})\right]\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{2})U({\boldsymbol{r}}_{7})\right]\right\rangle$
$\displaystyle=\left\langle\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{6})U({\boldsymbol{r}}_{1})\right]\right\rangle\left\langle\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{2})U({\boldsymbol{r}}_{7})\right]\right\rangle$
$\displaystyle+\mathcal{O}\left(\frac{1}{N_{c}}\right)$ (62)

Physically, we are approximating the four body average as the independent
broadening of the quark and anti-quark.

###  B.2 2-point functions

The first region $(x_{0}^{+},x_{v}^{+})$ contains the object
$\mathcal{S}^{(2)}$, which describes the momentum broadening of the gluon and
has a well known result in the literature. Let us review it:

$\mathcal{S}^{(2)}=\int_{\boldsymbol{r}_{1}(x_{0}^{+})=\boldsymbol{x_{g}}}^{\boldsymbol{r}_{1}(x_{v}^{+})=\boldsymbol{x_{v}}-\boldsymbol{x}}\mathcal{D}\boldsymbol{r}_{1}\int_{\boldsymbol{r}_{2}(x_{0}^{+})=\boldsymbol{\overline{x}_{g}}}^{\boldsymbol{r}_{2}(x_{v}^{+})=\boldsymbol{z_{1}}}\mathcal{D}\boldsymbol{r}_{2}\exp{\left\\{i\frac{q_{0}^{+}}{2}\int_{x_{0}^{+}}^{x_{v}^{+}}ds^{+}\left(\dot{\boldsymbol{r}}_{1}^{2}-\dot{\boldsymbol{r}}_{2}^{2}\right)\right\\}}\mathcal{C}^{(2)}\,.$
(63)

Following Eq. (61) and, e.g., Apolinario:2014csa , $\mathcal{C}^{(2)}$ reads

$\displaystyle\mathcal{C}^{(2)}$
$\displaystyle=\frac{1}{N_{c}^{2}-1}\left\langle\text{Tr}\left[U_{A}(\boldsymbol{r}_{1})U_{A}^{\dagger}({\boldsymbol{r}}_{2})\right]\right\rangle=\exp{\left\\{-\frac{C_{A}}{2}\int
ds^{+}\,n(s^{+})\sigma(\boldsymbol{r}_{1}-\boldsymbol{r}_{2})\right\\}}\,.$
(64)

By changing variables to
$\boldsymbol{u}=\boldsymbol{r}_{1}-\boldsymbol{r}_{2}$,
$\boldsymbol{v}=\boldsymbol{r}_{1}+\boldsymbol{r}_{2}$, we end up with

$\begin{split}\mathcal{S}^{(2)}&=\int\mathcal{D}\boldsymbol{u}\int\mathcal{D}\boldsymbol{v}\exp{\left\\{\int
ds^{+}\left(\frac{iq_{0}^{+}}{2}\dot{\boldsymbol{u}}\cdot\dot{\boldsymbol{v}}-\frac{C_{A}}{2}n(s^{+})\sigma(\boldsymbol{u})\right)\right\\}}\\\
&=\left(\frac{q_{0}^{+}}{2\pi\Delta_{0}^{+}}\right)^{2}\exp{\left\\{\frac{iq_{0}^{+}}{2\Delta_{0}^{+}}\left((\boldsymbol{x_{v}}-\boldsymbol{x}-\boldsymbol{x_{g}})^{2}-(\boldsymbol{z_{1}}-\boldsymbol{\overline{x}_{g}})^{2}\right)-\frac{C_{A}}{2}\int
ds^{+}\,n(s^{+})\sigma(\boldsymbol{u}^{0})\right\\}}\,,\end{split}$ (65)

where $\boldsymbol{u}^{0}$ is restricted to the path of a straight line due to
the integral in $\boldsymbol{v}$

$\boldsymbol{u}^{0}(s^{+})=\left((s^{+}-x_{0}^{+})(\boldsymbol{x_{v}}-\boldsymbol{x}-\boldsymbol{z_{1}})-(s^{+}-x_{v}^{+})(\boldsymbol{x_{g}}-\boldsymbol{\overline{x}_{g}})\right)/\Delta_{0}^{+}\,,$
(66)

and $\Delta_{0}^{+}=x_{v}^{+}-x_{0}^{+}$. Because the remaining of the squared
amplitude does not depend on the external transverse coordinates
$\boldsymbol{x_{g}}$ and $\boldsymbol{\overline{x}_{g}}$, one can explicitly
carry out these integrals. We start by changing variables to

$\begin{split}&\boldsymbol{\Delta_{g}}=\boldsymbol{x_{g}}-\boldsymbol{\overline{x}_{g}}\,,\quad\boldsymbol{\Sigma_{g}}=\boldsymbol{x_{g}}+\boldsymbol{\overline{x}_{g}}\,,\end{split}$
(67)

which is a transformation of jacobian $1/4$. Picking up the phases that depend
on $\boldsymbol{\Delta_{g}}$ and $\boldsymbol{\Sigma_{g}}$ in Eq. (16) and
realizing $\sigma(\boldsymbol{u}^{0})$ only depends on
$\boldsymbol{\Delta_{g}}$:

$\begin{split}&\int_{\boldsymbol{x_{g}},\boldsymbol{\overline{x}_{g}}}e^{i\boldsymbol{p}\cdot\boldsymbol{x_{g}}}e^{-i\bar{\boldsymbol{p}}\cdot\overline{x}_{g}}\mathcal{S}^{(2)}=\frac{1}{4}\left(\frac{q_{0}^{+}}{2\pi\Delta_{0}^{+}}\right)^{2}\int_{\boldsymbol{\Delta_{g}},\boldsymbol{\Sigma_{g}}}e^{i\boldsymbol{\Delta_{g}}\cdot\boldsymbol{\Sigma_{p}}/2}e^{i\boldsymbol{\Sigma_{g}}\cdot\boldsymbol{\Delta_{p}}/2}\\\
&\times\exp{\left\\{\frac{iq_{0}^{+}}{2\Delta_{0}^{+}}\Big{(}(\boldsymbol{x_{v}}-\boldsymbol{x})^{2}-\boldsymbol{z_{1}}^{2}-\boldsymbol{\Delta_{g}}\cdot(\boldsymbol{x_{v}}-\boldsymbol{x}+\boldsymbol{z_{1}})-\boldsymbol{\Sigma_{g}}\cdot(\boldsymbol{x_{v}}-\boldsymbol{x}-\boldsymbol{z_{1}}-\boldsymbol{\Delta_{g}})\Big{)}\right\\}}\\\
&\times\exp{\left\\{-\frac{C_{A}}{2}\int
ds^{+}\,n(s^{+})\sigma(\boldsymbol{u}^{0})\right\\}}\\\
&=e^{i\boldsymbol{\Sigma_{p}}\cdot({\boldsymbol{x}}_{v}-{\boldsymbol{x}}-{\boldsymbol{z}}_{1}-\Delta_{0}^{+}(\boldsymbol{\Delta_{p}}/q_{0}^{+}))/2}\,e^{i\boldsymbol{\Delta_{p}}\cdot(\boldsymbol{x_{v}}-{\boldsymbol{x}}+\boldsymbol{z_{1}})/2}\\\
&\times\exp{\left\\{-\frac{C_{A}}{2}\int
ds^{+}\,n(s^{+})\sigma(\boldsymbol{x_{v}}-\boldsymbol{x}-\boldsymbol{z_{1}}+(s^{+}-x_{v}^{+})\frac{\boldsymbol{\Delta_{p}}}{q_{0}^{+}})\right\\}}\,,\end{split}$
(68)

where we defined the placeholders
$\boldsymbol{\Sigma_{p}}={\boldsymbol{p}}+\bar{\boldsymbol{p}}$ and
$\boldsymbol{\Delta_{p}}={\boldsymbol{p}}-\bar{\boldsymbol{p}}$.

The last region $(\bar{x}_{v}^{+},L^{+})$ contains the object
$\mathcal{S}^{(4)}$, which in the large-$N_{c}$ limit simplifies to the
independent broadening of the quark and anti-quark (see Appendix B.1) and thus
to the product of two 2-point functions. It then reads

$\mathcal{S}^{(4)}=\mathcal{S}^{(2)}_{q}\mathcal{S}^{(2)}_{\bar{q}}+\mathcal{O}\left(\frac{1}{N_{c}}\right)\,,$
(69)

where

$\displaystyle\mathcal{S}^{(2)}_{q}=\int_{\boldsymbol{r}_{1}(\bar{x}_{v}^{+})=\boldsymbol{z_{3}}}^{\boldsymbol{r}_{1}(L^{+})=\boldsymbol{x_{q}}}\mathcal{D}\boldsymbol{r}_{1}\int_{\boldsymbol{r}_{2}(\bar{x}_{v}^{+})=\boldsymbol{\bar{x}_{v}-\bar{y}}}^{\boldsymbol{r}_{2}(L^{+})=\boldsymbol{\bar{x}_{q}}}\mathcal{D}\boldsymbol{r}_{2}\exp{\left\\{i\frac{q_{0}^{+}}{2}\int_{x_{0}^{+}}^{x_{v}^{+}}ds^{+}\left(\dot{\boldsymbol{r}}_{1}^{2}-\dot{\boldsymbol{r}}_{2}^{2}\right)\right\\}}\mathcal{C}^{(2)}_{q}$
$\displaystyle\mathcal{S}^{(2)}_{\bar{q}}=\int_{\boldsymbol{r}_{1}(\bar{x}_{v}^{+})=\boldsymbol{z_{2}}}^{\boldsymbol{r}_{1}(L^{+})=\boldsymbol{x_{\bar{q}}}}\mathcal{D}\boldsymbol{r}_{1}\int_{\boldsymbol{r}_{2}(\bar{x}_{v}^{+})=\boldsymbol{\bar{x}_{v}}}^{\boldsymbol{r}_{2}(L^{+})=\boldsymbol{\bar{x}_{\bar{q}}}}\mathcal{D}\boldsymbol{r}_{2}\exp{\left\\{i\frac{q_{0}^{+}}{2}\int_{x_{0}^{+}}^{x_{v}^{+}}ds^{+}\left(\dot{\boldsymbol{r}}_{1}^{2}-\dot{\boldsymbol{r}}_{2}^{2}\right)\right\\}}\mathcal{C}^{(2)}_{\bar{q}}\,,$
(70)

and

$\displaystyle\mathcal{C}^{(2)}_{q}=\frac{1}{N_{c}}\left\langle\text{Tr}\left[U({\boldsymbol{r}}_{1})U^{\dagger}({\boldsymbol{r}}_{2})\right]\right\rangle=\exp{\left\\{-\frac{C_{F}}{2}\int
ds^{+}\,n(s^{+})\sigma(\boldsymbol{r}_{1}-\boldsymbol{r}_{2})\right\\}}$
$\displaystyle\mathcal{C}^{(2)}_{\bar{q}}=\frac{1}{N_{c}}\left\langle\text{Tr}\left[U({\boldsymbol{r}}_{2})U^{\dagger}({\boldsymbol{r}}_{1})\right]\right\rangle=\exp{\left\\{-\frac{C_{F}}{2}\int
ds^{+}\,n(s^{+})\sigma(\boldsymbol{r}_{2}-\boldsymbol{r}_{1})\right\\}}\,.$
(71)

Picking up the phases in Eq. (16) attached to ${\boldsymbol{p}}_{1}$ and
${\boldsymbol{p}}_{2}$ and carrying out the calculation similarly to what was
done for the gluon broadening $\mathcal{S}^{2}$ but with the same momenta in
the direct and complex conjugate amplitudes, we get

$\begin{split}&\int_{\boldsymbol{x_{\bar{q}}},\boldsymbol{\overline{x}_{\bar{q}}},\boldsymbol{x_{q}},\boldsymbol{\overline{x}_{q}}}e^{-i\boldsymbol{p_{1}}\cdot(\boldsymbol{x_{q}}-\boldsymbol{\overline{x}_{q}})}e^{-i\boldsymbol{p_{2}}\cdot(\boldsymbol{x_{\bar{q}}}-\boldsymbol{\overline{x}_{\bar{q}}})}\mathcal{S}^{(4)}=e^{-i\boldsymbol{p_{1}}\cdot(\boldsymbol{z_{3}}-\boldsymbol{\overline{x}_{v}}+\boldsymbol{\overline{y}})}e^{-i\boldsymbol{p_{2}}\cdot(\boldsymbol{z_{2}}-\boldsymbol{\overline{x}_{v}})}\\\
&\times\exp{\left\\{-\frac{C_{F}}{2}\int_{\overline{x}_{v}^{+}}^{L^{+}}ds^{+}n(s^{+})(\sigma(\boldsymbol{z_{3}}-\boldsymbol{\overline{x}_{v}}+\boldsymbol{\overline{y}})+\sigma(\boldsymbol{\overline{x}_{v}}-\boldsymbol{z_{2}}))\right\\}}+\mathcal{O}\left(\frac{1}{N_{c}}\right)\,.\end{split}$
(72)

###  B.3 3-point function

The middle region $(x_{v}^{+},\bar{x}_{v}^{+})$ contains the object
$\mathcal{S}^{(3)}$ describing the production of the $q\bar{q}$ antenna inside
the medium. After changing the intermediate transverse variables from
$\\{{\boldsymbol{z}}_{1},{\boldsymbol{z}}_{2},{\boldsymbol{z}}_{3}\\}$ to

$\begin{split}&\boldsymbol{z_{1}}\rightarrow\boldsymbol{v_{1}}=\boldsymbol{x_{v}-{\boldsymbol{x}}-z_{1}}\\\
&\boldsymbol{z_{2}}\rightarrow\boldsymbol{v_{2}}=\boldsymbol{z_{2}-\overline{x}_{v}}\\\
&\boldsymbol{z_{3}}\rightarrow\boldsymbol{v_{3}}=\boldsymbol{z_{3}-\overline{x}_{v}+\overline{{\boldsymbol{y}}}}\,,\end{split}$
(73)

we write the three body average as

$\displaystyle\mathcal{S}^{(3)}$
$\displaystyle=\int_{\boldsymbol{r_{0}^{i}}}^{\boldsymbol{r_{0}^{f}}}\mathcal{D}\boldsymbol{r}_{0}\int_{\boldsymbol{r_{1}^{i}}}^{\boldsymbol{r_{1}^{f}}}\mathcal{D}\boldsymbol{r}_{1}\int_{\boldsymbol{r_{2}^{i}}}^{\boldsymbol{r_{2}^{f}}}\mathcal{D}\boldsymbol{r}_{2}\exp{\left\\{\frac{i}{2}\int_{x_{v}^{+}}^{\overline{x}_{v}^{+}}ds^{+}\left(-q_{0}^{+}\dot{\boldsymbol{r}}_{0}^{2}+p_{1}^{+}\dot{\boldsymbol{r}}_{1}^{2}+p_{2}^{+}\dot{\boldsymbol{r}}_{2}^{2}\right)\right\\}}\mathcal{C}^{(3)}\,,$
(74)

where

$\displaystyle\mathcal{C}^{(3)}=\frac{1}{2(N_{c}^{2}-1)}\left\langle\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{8})U({\boldsymbol{r}}_{3})\right]\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{4})U({\boldsymbol{r}}_{8})\right]-\frac{1}{N_{c}}\text{Tr}\left[U^{\dagger}({\boldsymbol{r}}_{4})U({\boldsymbol{r}}_{3})\right]\right\rangle$
(75)

and

$\begin{split}&\boldsymbol{r_{0}^{i}}=\boldsymbol{x_{v}-x-v_{1}}\\\
&\boldsymbol{r_{1}^{i}}=\boldsymbol{x_{v}-y}\\\
&\boldsymbol{r_{2}^{i}}=\boldsymbol{x_{v}}\end{split}\qquad\qquad\begin{split}&\boldsymbol{r_{0}^{f}}=\boldsymbol{\overline{x}_{v}-\overline{x}}\\\
&\boldsymbol{r_{1}^{f}}=\boldsymbol{v_{3}+\overline{x}_{v}-\overline{y}}\\\
&\boldsymbol{r_{2}^{f}}=\boldsymbol{v_{2}+\overline{x}_{v}}\,.\end{split}$
(76)

The explicit form of $\mathcal{C}^{(3)}$ is known in the literature (cf.
Isaksen:2020npj ; Apolinario:2014csa ). Hence, we have:

$\begin{split}\mathcal{S}^{(3)}&=\frac{1}{2}\int\mathcal{D}\boldsymbol{r}_{0}\int\mathcal{D}\boldsymbol{r}_{1}\int\mathcal{D}\boldsymbol{r}_{2}\exp{\left\\{\frac{i}{2}\int
ds^{+}\left(-q_{0}^{+}\dot{\boldsymbol{r}}_{0}^{2}+p_{1}^{+}\dot{\boldsymbol{r}}_{1}^{2}+p_{2}^{+}\dot{\boldsymbol{r}}_{2}^{2}\right)\right\\}}\\\
&\times\exp{\left\\{-\frac{1}{2}\int
ds^{+}\,n(s^{+})(N_{c}(\sigma(\boldsymbol{r}_{2}-\boldsymbol{r}_{0})+\sigma(\boldsymbol{r}_{1}-\boldsymbol{r}_{0}))-\frac{1}{N_{c}}\sigma(\boldsymbol{r}_{1}-\boldsymbol{r}_{2}))\right\\}}\,.\end{split}$
(77)

To further simplify these path integrals, we can make the following change of
variables:

$\boldsymbol{r}_{1}\rightarrow\boldsymbol{u}_{1}=\boldsymbol{r}_{2}-\boldsymbol{r}_{1}\,,\quad\boldsymbol{r}_{2}\rightarrow\boldsymbol{u}_{2}=-\boldsymbol{r}_{0}+z\boldsymbol{r}_{1}+(1-z)\boldsymbol{r}_{2}\\\
$ (78)

where $z=p_{1}^{+}/q_{0}^{+}$ as stated already. This leaves us with no
kinetic term and no $\sigma$ term on $\boldsymbol{r}_{0}$, allowing to perform
the integrals in $\boldsymbol{r}_{0}$ and $\boldsymbol{u}_{2}$, restricting
the latter to a straight line path (much like the integration in Eq. (67)):

$\begin{split}\mathcal{S}^{(3)}&=\frac{1}{2}\int\mathcal{D}\boldsymbol{r}_{0}\int\mathcal{D}\boldsymbol{u}_{1}\int\mathcal{D}\boldsymbol{u}_{2}\exp{\left\\{\frac{iq_{0}^{+}}{2}\int
ds^{+}\left(2\dot{\boldsymbol{r}}_{0}\cdot\dot{\boldsymbol{u}}_{2}+\dot{\boldsymbol{u}}_{2}^{2}+z(1-z)\dot{\boldsymbol{u}}_{1}^{2}\right)\right\\}}\\\
&\times\exp{\left\\{-\frac{1}{2}\int
ds^{+}\,n(s^{+})(N_{c}(\sigma(\boldsymbol{u}_{2}+z\boldsymbol{u}_{1})+\sigma(\boldsymbol{u}_{2}-(1-z)\boldsymbol{u}_{1}))-\frac{1}{N_{c}}\sigma(-\boldsymbol{u}_{1}))\right\\}}\\\
&=\frac{1}{2}\left(\frac{q_{0}^{+}}{2\pi\Delta
t}\right)^{2}\exp{\left\\{\frac{iq_{0}^{+}}{2}\frac{\Delta\boldsymbol{u}_{2}}{\Delta
t}\left(2\Delta\boldsymbol{r}_{0}+\Delta\boldsymbol{u}_{2}\right)\right\\}}\mathcal{K}(\boldsymbol{u}_{1}^{i},\boldsymbol{u}_{1}^{f},\boldsymbol{u}_{2}^{i},\boldsymbol{u}_{2}^{f})\,,\end{split}$
(79)

with $\tilde{q}_{0}^{+}=z(1-z)q_{0}^{+}$, $\Delta
t=\overline{x}_{v}^{+}-x_{v}^{+}$ and

$\begin{split}\boldsymbol{u}_{2}^{0}(s^{+})=((s^{+}-x_{v}^{+})\boldsymbol{u}_{2}^{f}-(s^{+}-\overline{x}_{v}^{+})\boldsymbol{u}_{2}^{i})/\Delta
t\,,\\\ \end{split}$
$\begin{split}&\boldsymbol{r_{0}^{i}}=\boldsymbol{x_{v}-x-v_{1}}\\\
&\boldsymbol{u_{1}^{i}}=\boldsymbol{y}\\\
&\boldsymbol{u_{2}^{i}}=\boldsymbol{v_{1}}+{\boldsymbol{x}}-z\boldsymbol{y}\end{split}\qquad\qquad\begin{split}&\boldsymbol{r_{0}^{f}}=\boldsymbol{\overline{x}_{v}-\overline{x}}\\\
&\boldsymbol{u_{1}^{f}}=\boldsymbol{v_{2}}-\boldsymbol{v_{3}}+\boldsymbol{\overline{y}}\\\
&\boldsymbol{u_{2}^{f}}=\boldsymbol{v_{2}}+\overline{{\boldsymbol{x}}}-z(\boldsymbol{v_{2}}-\boldsymbol{v_{3}}+\boldsymbol{\overline{y}})\,.\\\
\end{split}$ (80)

We have also defined the emission kernel $\mathcal{K}$ as

$\displaystyle\mathcal{K}(\boldsymbol{u}_{1}^{i},$
$\displaystyle\boldsymbol{u}_{1}^{f},\boldsymbol{u}_{2}^{i},\boldsymbol{u}_{2}^{f})\equiv\int\mathcal{D}\boldsymbol{u}_{1}\exp{\left\\{\frac{i\tilde{q}_{0}^{+}}{2}\int
ds^{+}\dot{\boldsymbol{u}}_{1}^{2}\right\\}}$
$\displaystyle\times\exp{\left\\{-\frac{1}{2}\int
ds^{+}\,n(s^{+})(N_{c}(\sigma(\boldsymbol{u}_{2}^{0}+z\boldsymbol{u}_{1})+\sigma(\boldsymbol{u}_{2}^{0}-(1-z)\boldsymbol{u}_{1}))-\frac{1}{N_{c}}\sigma(-\boldsymbol{u}_{1}))\right\\}}\,.$
(81)

It is interesting to note that $\mathcal{S}^{(3)}$ only depends on
${\boldsymbol{x}}_{v}$ and $\bar{\boldsymbol{x}}_{v}$ through the difference
$\Delta_{v}=\bar{\boldsymbol{x}}_{v}-{\boldsymbol{x}}_{v}$ in
$\boldsymbol{\Delta r}_{0}$. Aside from this dependence, after applying the
change of variables in Eq. (73) to the results in Eq. (68) and (72), we see
that the only dependence on the vertex transverse positions that is left is
via a complex phase
$e^{i({\boldsymbol{p}}-\bar{\boldsymbol{p}})\cdot{\boldsymbol{x}}_{v}}$ in Eq.
(68). This allows one to integrate over both ${\boldsymbol{x}}_{v}$ and
$\bar{\boldsymbol{p}}$ to set ${\boldsymbol{p}}=\bar{\boldsymbol{p}}$, i.e.,
the same transverse momentum flows out of the initial current in amplitude and
its complex conjugate. This leads to the result in Eq. (18).

## References

  * (1) PHENIX Collaboration, K. Adcox et al., Suppression of hadrons with large transverse momentum in central Au+Au collisions at $\sqrt{s_{NN}}$ = 130-GeV, Phys. Rev. Lett. 88 (2002) 022301, [nucl-ex/0109003].
  * (2) STAR Collaboration, J. Adams et al., Experimental and theoretical challenges in the search for the quark gluon plasma: The STAR Collaboration’s critical assessment of the evidence from RHIC collisions, Nucl. Phys. A 757 (2005) 102–183, [nucl-ex/0501009].
  * (3) STAR Collaboration, J. Adam et al., Measurement of inclusive charged-particle jet production in Au + Au collisions at $\sqrt{s_{NN}}=$200 GeV, Phys. Rev. C 102 (2020), no. 5 054913, [arXiv:2006.00582].
  * (4) B. Muller, J. Schukraft, and B. Wyslouch, First Results from Pb+Pb collisions at the LHC, Ann. Rev. Nucl. Part. Sci. 62 (2012) 361–386, [arXiv:1202.3233].
  * (5) ATLAS Collaboration, M. Aaboud et al., Measurement of the nuclear modification factor for inclusive jets in Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}}=5.02$ TeV with the ATLAS detector, Phys. Lett. B 790 (2019) 108–128, [arXiv:1805.05635].
  * (6) CMS Collaboration, A. M. Sirunyan et al., First measurement of large area jet transverse momentum spectra in heavy-ion collisions, JHEP 05 (2021) 284, [arXiv:2102.13080].
  * (7) W. Busza, K. Rajagopal, and W. van der Schee, Heavy Ion Collisions: The Big Picture, and the Big Questions, Ann. Rev. Nucl. Part. Sci. 68 (2018) 339–376, [arXiv:1802.04801].
  * (8) S. Schlichting and D. Teaney, The First fm/c of Heavy-Ion Collisions, Ann. Rev. Nucl. Part. Sci. 69 (2019) 447–476, [arXiv:1908.02113].
  * (9) J. Berges, M. P. Heller, A. Mazeliauskas, and R. Venugopalan, QCD thermalization: Ab initio approaches and interdisciplinary connections, Rev. Mod. Phys. 93 (2021), no. 3 035003, [arXiv:2005.12299].
  * (10) A. Jaiswal and V. Roy, Relativistic hydrodynamics in heavy-ion collisions: general aspects and recent developments, Adv. High Energy Phys. 2016 (2016) 9623034, [arXiv:1605.08694].
  * (11) G.-Y. Qin and X.-N. Wang, Jet quenching in high-energy heavy-ion collisions, Int. J. Mod. Phys. E 24 (2015), no. 11 1530014, [arXiv:1511.00790].
  * (12) L. Cunqueiro and A. M. Sickles, Studying the QGP with Jets at the LHC and RHIC, Prog. Part. Nucl. Phys. 124 (2022) 103940, [arXiv:2110.14490].
  * (13) L. Apolinário, Y.-J. Lee, and M. Winn, Heavy quarks and jets as probes of the QGP, Prog. Part. Nucl. Phys. 127 (2022) 103990, [arXiv:2203.16352].
  * (14) E. Wang and X.-N. Wang, Jet tomography of dense and nuclear matter, Phys. Rev. Lett. 89 (2002) 162301, [hep-ph/0202105].
  * (15) I. Vitev, Jet tomography, J. Phys. G 30 (2004) S791–S800, [hep-ph/0403089].
  * (16) N. Armesto, C. A. Salgado, and U. A. Wiedemann, Measuring the collective flow with jets, Phys. Rev. Lett. 93 (2004) 242301, [hep-ph/0405301].
  * (17) S. Hauksson and E. Iancu, Jet polarisation in an anisotropic medium, JHEP 08 (2023) 027, [arXiv:2303.03914].
  * (18) A. Baty, P. Gardner, and W. Li, Novel observables for exploring QCD collective evolution and quantum entanglement within individual jets, Phys. Rev. C 107 (2023), no. 6 064908, [arXiv:2104.11735].
  * (19) W. Zhao, Z.-W. Lin, and X.-N. Wang, Collectivity inside high-multiplicity jets in high-energy proton-proton collisions, arXiv:2401.13137.
  * (20) K. Boguslavski, A. Kurkela, T. Lappi, F. Lindenbauer, and J. Peuron, Jet quenching parameter in QCD kinetic theory, arXiv:2312.00447.
  * (21) K. Boguslavski, A. Kurkela, T. Lappi, F. Lindenbauer, and J. Peuron, Jet momentum broadening during initial stages in heavy-ion collisions, Phys. Lett. B 850 (2024) 138525, [arXiv:2303.12595].
  * (22) A. Ipp, D. I. Müller, and D. Schuh, Anisotropic momentum broadening in the 2+1D Glasma: analytic weak field approximation and lattice simulations, Phys. Rev. D 102 (2020), no. 7 074001, [arXiv:2001.10001].
  * (23) A. Ipp, D. I. Müller, and D. Schuh, Jet momentum broadening in the pre-equilibrium Glasma, Phys. Lett. B 810 (2020) 135810, [arXiv:2009.14206].
  * (24) D. Avramescu, V. Băran, V. Greco, A. Ipp, D. I. Müller, and M. Ruggieri, Simulating jets and heavy quarks in the glasma using the colored particle-in-cell method, Phys. Rev. D 107 (2023), no. 11 114021, [arXiv:2303.05599].
  * (25) M. E. Carrington, A. Czajka, and S. Mrowczynski, Jet quenching in glasma, Phys. Lett. B 834 (2022) 137464, [arXiv:2112.06812].
  * (26) M. E. Carrington, A. Czajka, and S. Mrowczynski, Transport of hard probes through glasma, Phys. Rev. C 105 (2022), no. 6 064910, [arXiv:2202.00357].
  * (27) J. Barata, S. Hauksson, X. Mayo López, and A. V. Sadofyev, Jet quenching in the glasma phase: medium-induced radiation, arXiv:2406.07615.
  * (28) M. V. Kuzmin, X. Mayo López, J. Reiten, and A. V. Sadofyev, Jet quenching in anisotropic flowing matter, Phys. Rev. D 109 (2024), no. 1 014036, [arXiv:2309.00683].
  * (29) J. Barata, J. G. Milhano, and A. V. Sadofyev, Picturing QCD jets in anisotropic matter: from jet shapes to energy energy correlators, Eur. Phys. J. C 84 (2024), no. 2 174, [arXiv:2308.01294].
  * (30) J. Barata, X. Mayo López, A. V. Sadofyev, and C. A. Salgado, Medium induced gluon spectrum in dense inhomogeneous matter, Phys. Rev. D 108 (2023), no. 3 034018, [arXiv:2304.03712].
  * (31) J. Barata, A. V. Sadofyev, and X.-N. Wang, Quantum partonic transport in QCD matter, Phys. Rev. D 107 (2023), no. 5 L051503, [arXiv:2210.06519].
  * (32) C. Andres, F. Dominguez, A. V. Sadofyev, and C. A. Salgado, Jet broadening in flowing matter: Resummation, Phys. Rev. D 106 (2022), no. 7 074023, [arXiv:2207.07141].
  * (33) J. Barata, A. V. Sadofyev, and C. A. Salgado, Jet broadening in dense inhomogeneous matter, Phys. Rev. D 105 (2022), no. 11 114010, [arXiv:2202.08847].
  * (34) A. V. Sadofyev, M. D. Sievert, and I. Vitev, Ab initio coupling of jets to collective flow in the opacity expansion approach, Phys. Rev. D 104 (2021), no. 9 094044, [arXiv:2104.09513].
  * (35) Y. He, L.-G. Pang, and X.-N. Wang, Gradient Tomography of Jet Quenching in Heavy-Ion Collisions, Phys. Rev. Lett. 125 (2020), no. 12 122301, [arXiv:2001.08273].
  * (36) Y.-X. Xiao, Y. He, L.-G. Pang, H. Zhang, and X.-N. Wang, Asymmetric jet shapes with 2D jet tomography, arXiv:2402.00264.
  * (37) Y. Fu, J. Casalderrey-Solana, and X.-N. Wang, Asymmetric transverse momentum broadening in an inhomogeneous medium, Phys. Rev. D 107 (2023), no. 5 054038, [arXiv:2204.05323].
  * (38) Y. He, W. Chen, T. Luo, S. Cao, L.-G. Pang, and X.-N. Wang, Event-by-event jet anisotropy and hard-soft tomography of the quark-gluon plasma, Phys. Rev. C 106 (2022), no. 4 044904, [arXiv:2201.08408].
  * (39) M. V. Kuzmin and X. Mayo López, Gluon radiation inside a flowing medium, arXiv:2406.14628.
  * (40) S. Hauksson, S. Jeon, and C. Gale, Probes of the quark-gluon plasma and plasma instabilities, Phys. Rev. C 103 (2021) 064904, [arXiv:2012.03640].
  * (41) S. Hauksson, S. Jeon, and C. Gale, Momentum broadening of energetic partons in an anisotropic plasma, Phys. Rev. C 105 (2022), no. 1 014914, [arXiv:2109.04575].
  * (42) B. G. Zakharov, Fully quantum treatment of the Landau-Pomeranchuk-Migdal effect in QED and QCD, JETP Lett. 63 (1996) 952–957, [hep-ph/9607440].
  * (43) R. Baier, Y. L. Dokshitzer, A. H. Mueller, S. Peigne, and D. Schiff, Radiative energy loss of high-energy quarks and gluons in a finite volume quark - gluon plasma, Nucl. Phys. B 483 (1997) 291–320, [hep-ph/9607355].
  * (44) M. Attems, J. Brewer, G. M. Innocenti, A. Mazeliauskas, S. Park, W. van der Schee, and U. A. Wiedemann, The medium-modified $g\to c\overline{c}$ splitting function in the BDMPS-Z formalism, JHEP 01 (2023) 080, [arXiv:2203.11241].
  * (45) Y. Mehtar-Tani, J. G. Milhano, and K. Tywoniuk, Jet physics in heavy-ion collisions, Int. J. Mod. Phys. A 28 (2013) 1340013, [arXiv:1302.2579].
  * (46) J.-P. Blaizot and Y. Mehtar-Tani, Jet Structure in Heavy Ion Collisions, Int. J. Mod. Phys. E 24 (2015), no. 11 1530012, [arXiv:1503.05958].
  * (47) J.-P. Blaizot, F. Dominguez, E. Iancu, and Y. Mehtar-Tani, Medium-induced gluon branching, JHEP 01 (2013) 143, [arXiv:1209.4585].
  * (48) J. Barata, J.-P. Blaizot, and Y. Mehtar-Tani, Quantum to classical parton dynamics in QCD media, Phys. Rev. D 108 (2023), no. 1 014039, [arXiv:2305.10476].
  * (49) J.-P. Blaizot and M. A. Escobedo, Quantum and classical dynamics of heavy quarks in a quark-gluon plasma, JHEP 06 (2018) 034, [arXiv:1711.10812].
  * (50) J.-P. Blaizot and E. Iancu, Lifetimes of quasiparticles and collective excitations in hot QED plasmas, Phys. Rev. D 55 (1997) 973–996, [hep-ph/9607303].
  * (51) E. Braaten and R. D. Pisarski, Calculation of the gluon damping rate in hot QCD, Phys. Rev. D 42 (1990) 2156–2160.
  * (52) Y. Mehtar-Tani, Gluon bremsstrahlung in finite media beyond multiple soft scattering approximation, JHEP 07 (2019) 057, [arXiv:1903.00506].
  * (53) J. Barata and Y. Mehtar-Tani, Improved opacity expansion at NNLO for medium induced gluon radiation, JHEP 10 (2020) 176, [arXiv:2004.02323].
  * (54) L. Apolinário, N. Armesto, J. G. Milhano, and C. A. Salgado, Medium-induced gluon radiation and colour decoherence beyond the soft approximation, JHEP 02 (2015) 119, [arXiv:1407.0599].
  * (55) H. Kleinert, Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets, .
  * (56) R. Snellings, Elliptic Flow: A Brief Review, New J. Phys. 13 (2011) 055008, [arXiv:1102.3010].
  * (57) D. E. Kharzeev, J. Liao, S. A. Voloshin, and G. Wang, Chiral magnetic and vortical effects in high-energy nuclear collisions—A status report, Prog. Part. Nucl. Phys. 88 (2016) 1–28, [arXiv:1511.04050].
  * (58) A. Kumar, B. Müller, and D.-L. Yang, Spin polarization and correlation of quarks from the glasma, Phys. Rev. D 107 (2023), no. 7 076025, [arXiv:2212.13354].
  * (59) X. Du, Heavy quark drag and diffusion coefficients in the prehydrodynamic QCD plasma, Phys. Rev. C 109 (2024), no. 1 014901, [arXiv:2306.02530].
  * (60) D. Binosi and L. Theussl, JaxoDraw: A Graphical user interface for drawing Feynman diagrams, Comput. Phys. Commun. 161 (2004) 76–86, [hep-ph/0309015].
  * (61) K. Fukushima and Y. Hidaka, Light projectile scattering off the color glass condensate, JHEP 06 (2007) 040, [arXiv:0704.2806].
  * (62) K. Fukushima and Y. Hidaka, General formulae for dipole Wilson line correlators with the Color Glass Condensate, JHEP 11 (2017) 114, [arXiv:1708.03051].
  * (63) J. H. Isaksen and K. Tywoniuk, Wilson line correlators beyond the large-Nc, JHEP 21 (2020) 125, [arXiv:2107.02542].
  * (64) F. Dominguez, C. Marquet, A. M. Stasto, and B.-W. Xiao, Universality of multiparticle production in QCD at high energies, Phys. Rev. D 87 (2013) 034007, [arXiv:1210.1141].
  * (65) F. Domínguez, J. G. Milhano, C. A. Salgado, K. Tywoniuk, and V. Vila, Mapping collinear in-medium parton splittings, Eur. Phys. J. C 80 (2020), no. 1 11, [arXiv:1907.03653].