aainstitutetext: Department of Physics and Technology, University of Bergen,
5007 Bergen, Norway

# Wilson line correlators beyond the large-$N_{c}$

Johannes Hamre Isaksen a and Konrad Tywoniuk<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We study hard $1\to 2$ final-state parton splittings in the medium, and put
special emphasis on calculating the Wilson line correlators that appear in
these calculations. As partons go through the medium their color continuously
rotates, an effect that is encapsulated in a Wilson line along their
trajectory. When calculating observables, one typically has to calculate
traces of two or more medium-averaged Wilson lines. These are usually dealt
with in the literature by invoking the large-$N_{c}$ limit, but exact
calculations have been lacking in many cases. In our work, we show how
correlators of multiple Wilson lines appear, and develop a method to calculate
them numerically to all orders in $N_{c}$. Initially, we focus on the trace of
four Wilson lines, which we develop a differential equation for. We will then
generalize this calculation to a product of an arbitrary number of Wilson
lines, and show how to do the exact calculation numerically, and even
analytically in the large-$N_{c}$ limit. Color sub-leading corrections, that
are suppressed with a factor $N_{c}^{-2}$ relative to the leading scaling, are
calculated explicitly for the four-point correlator and we discuss how to
extend this method to the general case. These results are relevant for
high-$p_{T}$ jet processes and initial stage physics at the LHC.

##  1 Introduction

One of the primary reasons for colliding heavy ions with ultra-relativistic
energies is to probe QCD matter in an extremely hot and dense phase, called
the quark-gluon plasma (QGP). There are multiple ways to probe and learn about
the properties of QGP, utilizing the properties of bulk particle production
and rare probes. In one example of the latter category, the heavy-ion
collision involves a hard partonic sub-collision that produces hard partons
that propagate through the medium and escape to the detectors as jets. The
study of how the properties of these jets change as they go through the
medium, colloquially referred to as “jet quenching,” is a versatile tool to
study hot QCD matter dEnterria:2009xfs ; Majumder:2010qh ; Mehtar-Tani:2013pia
.

Experiments at RHIC Adams:2005dq ; Adcox:2004mh  and the LHC Aamodt:2010jd ;
Khachatryan:2016odn ; Chatrchyan:2011sx ; Aad:2010bu ; Abelev:2013kqa
colliders have found strong suppression of high-$p_{T}$ particles in heavy-ion
collisions compared to proton-proton collisions, which is interpreted as a
clear sign of the energy loss of jets that suffer final-state interactions
with the surrounding QGP. On the theoretical side, this is interpreted in
terms of radiative energy loss, where particles in the jet lose energy through
medium-induced emission of gluons that end up outside of the reconstructed jet
cone, and elastic drag. For large media, as typically encountered in central
to semi-central lead-lead collisions, it is the former process that dominates
the total lost energy.

The energy loss process for single partons is well understood since many
years, see, e.g., Baier:1994bd ; Baier:1996sk ; Baier:1996kr ; Baier:1998yf ;
Zakharov:1996fv ; Zakharov:1997uu ; Wiedemann:2000za . However, a jet is a
more complicated composite object consisting of several hard partons. A hard
parton propagating through the medium will typically undergo several
splittings, resulting in a multi-parton state that will interact differently
with the medium compared to how the individual partons would. Such splittings
can occur as long as the scale of the splittings, for instance the generated
relative transverse momentum in the splitting, is bigger than what the medium
can supply through multiple scattering. In particular, the modifications of
effects of color coherence play an important role in determining which
emissions will be resolved by the medium and contribute toward the total
energy loss Mehtar_Tani_2018 ; Mehtar-Tani:2019tvy ; Caucal:2018dla ;
Caucal:2020xad . Instead of focusing on single partons, we will study a hard
parton splitting into two, and their subsequent propagation through the
medium. This is certainly a better approximation of a real jet than a single
parton, and has the additional advantage that one can build up jets from
several partons by consecutive $1\to 2$ splittings.

Previous studies of such processes focused mostly on a hard photon splitting
into a quark-antiquark pair Mehtar_Tani_2018 ; Dom_nguez_2020  and invoked the
large-$N_{c}$ approximation to obtain analytical formulas. In this work, we
consider three generic QCD splitting processes that involve up to eight
correlated Wilson lines in the fundamental representation, in the case of
gluon splitting into two daughter gluons. Our specific improvement concerns a
more precise way to calculate correlators of Wilson lines that often appear in
these calculations, and it can, in principle, be extended for an arbitrary
number of propagating particles through the medium.

To give a general flavor of how our procedure works, recall that a matrix
element generally involves several propagators that resum multiple scattering
through Wilson lines $V$, which extend along the trajectories in the medium.
Ignoring some factors irrelevant for the present discussion, the matrix
element squared will take the following simplified form

$\displaystyle\left\langle|{\cal M}|^{2}\right\rangle$
$\displaystyle\sim\langle\operatorname{tr}[V^{\dagger}V\ldots
V^{\dagger}V]\ldots\operatorname{tr}[V^{\dagger}V\ldots
V^{\dagger}V]\rangle\,,$ (1)

where the angular brackets denote an average over medium configurations. For a
generic $1\to 2$ process, the amplitude squared can be reduced to a product of
two-, three- and four-point correlators Blaizot_2013 ; Apolinario:2014csa . To
calculate these processes it is imperative to know the form of the Wilson line
correlator appearing on the right hand side of (1), which we will denote by
the letter $C^{K}$ for correlator, where the subscript $K$ refers to the
number of traces. If you assume that the number of colors $N_{c}$ is large the
calculation of these correlators usually simplifies sufficiently to be
possible to calculate. Namely, the leading $N_{c}$ scaling emerges from
simplifying the medium averages to
$\langle\operatorname{tr}[V^{\dagger}V\ldots
V^{\dagger}V]\rangle\ldots\langle\operatorname{tr}[V^{\dagger}V\dots
V^{\dagger}V]\rangle$, which scales like $N_{c}^{K}$. However, since $N_{c}=3$
is not a very large number it is sensible to ask whether this approximation is
sound or not. As we will see, the terms that are discarded by performing the
large-$N_{c}$ approximation will be smaller than the other terms by a factor
$\sim 1/N_{c}^{2}\simeq 10\%$ for typical situations. However, evaluating the
correlators at large times, could lead to big discrepancies between the finite
and large-$N_{c}$ calculations.

In this paper we will develop a method for calculating correlators of an
arbitrary number of Wilson lines at finite $N_{c}$, which casts their
evolution and mixing in terms of a coupled evolution equation in time
(referring to their trajectories through the medium). This reduces the
complexity of the formulation compared to previous calculations of multi-
Wilson line correlators, see Kovner_2001  for a technique based on
diagonalization of the evolution matrix and Dominguez:2011wm ;
Apolinario:2014csa  for an iterative procedure. The derivation of the
evolution matrix culminates in Eq. (4.2). This allows us to evaluate these
correlators at an arbitrary time, and can be addressed using numerical
techniques. We also consider in detail the large-$N_{c}$ approximation, which
leads to a striking simplification of the dynamics since all higher-order
correlators can be calculated using two-point correlators (dipoles) and their
convolutions. Furthermore, we have computed the sub-leading correction in
color. Considering again the generic correlator $C^{K}$ from Eq. (1) above,
the generic expansion in $N_{c}$ takes the following form,

$C^{K}=N_{c}^{K}\hat{C}^{K}_{\text{leading
$N_{c}$}}+N_{c}^{K-2}\hat{C}^{K}_{\text{sub-leading
$N_{c}$}}+\mathcal{O}(N_{c}^{K-4})\,,$ (2)

where the two first terms can be found analytically (the hat over the
correlators imply that we have explicitly extracted their leading $N_{c}$
behavior). It turns out that, in many cases, $C^{K}_{\text{sub-leading
$N_{c}$}}$ is essential to recover the correct long-time behavior of the
correlators.

We will explore how big the error is by comparing the exact results to the
large-$N_{c}$ approximation in realistic settings in high-energy jet
splittings. We mainly consider hard emissions early in the medium, i.e. at
scales much larger than those provided by the medium, and therefore we neglect
any broadening of the particles. The daughters are traversing the medium at a
fixed angle, or “tilt”, given by the kinematics of the hard splitting (we fix
our coordinate system so that the parent particle has zero angle). For
splittings where at least one of the daughters becomes very soft or is being
emitted at a large angle, one should also allow for additional transverse
momentum broadening, as done in Blaizot_2013 ; Apolinario:2014csa , albeit
only in the large-$N_{c}$ approximation. We have left this additional
complication for future work.

Our calculation is also very pertinent for improving our understanding of
color dynamics in the medium, for instance in the context of multi-gluon
emissions with overlapping formation times Arnold:2019qqc  and to understand
hadronization after exiting the QGP Zakharov:2018hfz . In the process of
evaluation of the multi-Wilson line correlators, the only assumption made is
the exact form of the medium average, see Eq. (7), which is also employed in
other contexts than for a thermal medium, see, e.g., Hatta:2020wre  for
calculating such correlators on the lattice. Therefore, although we have
derived our method of calculating Wilson line correlators in the context of
jet quenching, it is a general result that can be applied in more branches of
QCD. One concrete example refer to initial state physics, where multi-particle
production is considered an important channel to verify saturation effects in
the nuclei Kovner_2001 ; Jalilian-Marian:2004vhw ; Iancu:2011ns . Furthermore,
sub-leading corrections in color have also been considered in the context of
high-energy QCD evolution at next-to-leading order Lappi:2020srm . Finally,
the generic color structure of high-energy QCD events is actively studied
Dominguez:2011wm ; Dominguez:2012ad . It is also interesting to note that sub-
leading color corrections have been considered in the context of improving
parton showers in the vacuum, see, e.g., Nagy:2012bt ; Hamilton:2020rcu .

Let us briefly outline the structure of the paper. Section 2 introduces the
notation and formalism we will make use of throughout the paper. In Sec. 3 we
will consider three examples of splitting processes that lead to Wilson line
correlators: a photon producing a quark-antiquark pair, a quark emitting a
gluon and a gluon splitting into two gluons. Those processes will provide the
motivation for the rest of the calculation in the paper. In Sec. 4.1, we will
go into detail about calculating the simplest of the Wilson lines structures
from Sec. 3, which is a trace of four lines. Here, we also develop a method to
compute the color sub-leading corrections, corresponding to the second term on
the right hand side in Eq. (2). Thereafter, in Sec. 4 we will generalize the
method used in Sec. 4.1 to correlators of an arbitrary number of Wilson lines,
and show how one can always make a system of differential equations to
describe these structures. This section contains the main theoretical results
of the paper. The formulas developed in Sec. 4 are used to calculate the more
complicated Wilson line structures appearing in Sec. 3. We will show how the
calculations simplify in the large-$N_{c}$ approximation, and use numerical
evaluation to compare the approximate results to the exact ones.

##  2 Basic elements and notation

We will assume that the partons propagating through the medium are highly
energetic and travelling on the light-cone almost strictly in the positive $z$
direction. In light-cone (LC) coordinates it will have momentum
$(p^{+},p^{-},{\bm{p}})$, where $p^{+}=(p^{0}+p^{3})/2$ is identified with the
LC energy $E\equiv p^{+}$, $p^{-}=p^{0}-p^{3}$ is negligible and ${\bm{p}}$ is
the transverse momentum. The parton interacts with the medium, which is
modelled by a classical background gauge field $A^{\mu,a}(t,{\bm{r}})$. The
interaction of the parton with the classical field leads to transverse
momentum broadening and energy loss. The interactions can be resummed using a
framework developed by Baier-Dokshitzer-Mueller-Peigné-Schiff Baier:1994bd ;
Baier:1996sk ; Baier:1996kr ; Baier:1998yf  and Zakharov Zakharov:1996fv ;
Zakharov:1997uu , and is known as the BDMPS-Z formalism. For small media,
where interactions are rare, this is equivalent with considering only one
interaction, known as the Gyulassy-Levai-Vitev (GLV) Gyulassy:2000er
approximation.

It is possible to construct Feynman rules from the BDMPS-Z approach, with
special in-medium propagators and vertices Mehtar_Tani_2018 . In this
formulation a highly energetic parton travelling through the medium can be
described by the propagator

$({\bm{x}}|{\cal
G}_{R}(t,t_{0})|{\bm{x}}_{0})=\Theta(t-t_{0})\,\int^{{\bm{x}}}_{{\bm{x}}_{0}}{\cal
D}{\bm{r}}\exp\left[i\frac{E}{2}\int_{t_{0}}^{t}{\rm
d}s\,\dot{\bm{r}}^{2}(s)\right]\,V_{R}(t,t_{0};{\bm{r}}(t))\,.$ (3)

In this expression $V_{R}$ is a Wilson line in the representation $R$, which
is given by

$V_{R}\left(t,t_{0};\mathbf{r}(t)\right)=\mathcal{P}\exp\left[ig\int_{t_{0}}^{t}\mathrm{d}s\,A^{a}(s,\mathbf{r}(s))T^{a}_{R}\right]\,,$
(4)

where the symbol $\mathcal{P}$ enforces path ordering. A quark transforms in
the fundamental representation, so the group generator is $T_{F}^{a}\equiv
t^{a}_{ij}$. Similarly, a gluon transforms in the adjoint representation, and
its group generator is $T_{A}^{a}\equiv(T^{a})^{bc}=-if^{abc}$. The final
results in this paper will mainly concern fundamental lines, which we will
denote by $V\equiv V_{F}$. Similarly, we will write the adjoint lines as
$U\equiv V_{A}$. Focusing on fundamental lines is sufficient, since one can
always transform adjoint Wilson lines to fundamental ones through the identity

$U^{ab}=2\operatorname{tr}\left[t^{a}Vt^{b}V^{\dagger}\right]=U^{\dagger
ba}\,.$ (5)

In the absence of interactions, i.e. when the Wilson line is evaluated at
$g=0$, we simply get

$({\bm{x}}|{\cal G}_{0}(t,t_{0})|{\bm{x}}_{0})\equiv{\cal
G}_{0}({\bm{x}}-{\bm{x}}_{0},t-t_{0})=\Theta(t-t_{0})\frac{E}{2\pi
i(t-t_{0})}{\rm
e}^{i\frac{E}{2}\frac{({\bm{x}}-{\bm{x}}_{0})^{2}}{(t-t_{0})}}\,,$ (6)

which is a representation of the retarded part of the Feynman propagator
($E>0$).

As mentioned in the introduction, the matrix element describing final-state
interactions in the QGP will involve one or more propagators of the form in
Eq. (3). Hence, on the level of the matrix element squared, we have to compute
correlators of such lines averaged over all possible medium configurations.
The medium average is indicated by $\langle\dots\rangle$, and we assume that
the correlator of the medium fields takes the form

$\langle
A^{a}(t,{\bm{r}})A^{b}(t^{\prime},{\bm{r}}^{\prime})\rangle=\delta^{ab}n(t)\delta(t-t^{\prime})\gamma\left(\bm{r}-\bm{r}^{\prime}\right)\,,$
(7)

which corresponds to the Gaussian noise approximation. Here, $n(t)$ is the
(time-dependent) density of scattering centers in the medium and

$\gamma({\bm{r}})=\int\frac{{\rm d}^{2}{\bm{q}}}{(2\pi)^{2}}\,{\rm
e}^{i{\bm{q}}\cdot{\bm{r}}}\frac{{\rm d}^{2}\sigma_{\rm el}}{{\rm
d}^{2}{\bm{q}}}\sim g^{2}\int\frac{{\rm
d}^{2}{\bm{q}}}{(2\pi)^{2}}\,\frac{{\rm
e}^{i{\bm{q}}\cdot{\bm{r}}}}{{\bm{q}}^{4}}\,,$ (8)

is the Fourier transform of the in-medium elastic scattering potential, where
the infrared behavior of the potential is regulated by an in-medium screening
mass. The delta function in time indicates that we have assumed the medium
interactions to be instantaneous. In many cases it will be convenient to
define

$\sigma({\bm{r}})=g^{2}\big{[}\gamma(\bm{0})-\gamma(\bm{r})\big{]}\,.$ (9)

The form of the function $\sigma$ depends on how the medium is modelled. The
two main ways of calculating this is through the Gyulassy-Wang model
Wang:1991xy  or through Hard Thermal Loop theory Aurenche:2002pd . These
models differ mainly in how infrared screening is implemented when
$q_{\perp}\to 0$. In this paper, we will however work in the harmonic
oscillator approximation, which accounts for multiple soft interactions. In
this case, the potential $\sigma({\bm{r}})$ can be cast as

$\displaystyle
C_{R}n\sigma(\bm{r})\simeq\frac{1}{4}{\bm{r}}^{2}\hat{q}_{R}(t)\,,$ (10)

where

$\hat{q}_{R}=C_{R}ng^{2}\int^{q^{\rm max}_{\perp}}\frac{{\rm
d}^{2}{\bm{q}}}{(2\pi)^{2}}\,{\bm{q}}^{2}\frac{{\rm d}^{2}\sigma_{\rm
el}}{{\rm d}^{2}{\bm{q}}}\,,$ (11)

is the jet quenching coefficient where $R$ denotes the color representation of
the Wilson lines. For the fundamental and adjoint representations we have
$C_{F}=\frac{N_{c}^{2}-1}{2N_{c}}$ and $C_{A}=N_{c}$, respectively. In this
paper we will use ${\hat{q}}={\hat{q}}_{F}$ unless otherwise stated. In Eq.
(11) we have explicitly introduced a UV cut-off to regularize the integral. A
more systematic approach to the regularization of the integral, and the
extension beyond the soft scattering approximation, has been pursued in Refs.
Mehtar-Tani:2019tvy ; Mehtar_Tani_2020 ; Barata:2020rdn .

We stress that the approximation in (10) is not necessary in order to solve
numerically the system of equations for arbitrary $n$-point correlators, but
it is very useful to employ to compare these exact results to analytical
calculations of the leading and sub-leading color correlators.

In the current work, we will focus on hard $1\to 2$ splitting processes in the
medium, where the initial particle has energy $E$ and the two splitting
products carry, respectively, $\omega_{1}=(1-z)E$ and $\omega_{2}=zE$. This is
formally equivalent to setting the energy of the mother particle,
$E\to\infty$, and considering a finite momentum sharing fraction $0\ll z\ll
1$. These conditions enforce that both the mother and daughter particles
travel on classical paths. Concretely, the trajectory of a particle in the
medium between time $t_{0}$ and $t$, given by the the propagator
$({\bm{x}}|{{\cal G}}(t,t_{0})|{\bm{x}}_{0})$, in configuration space, for
$E\gg(t-t_{0})^{-1}$ gets strongly constrained to the classical path
connecting the initial and final transverse positions, see Eq. (3), and leads
to

$({\bm{x}}|{\cal G}_{R}(t,t_{0})|{\bm{x}}_{0})\simeq{\cal
G}_{0}({\bm{x}}-{\bm{x}}_{0},t-t_{0})\,V_{R}(t,t_{0};[{\bm{x}}_{cl}(s)])\,,$
(12)

where the classical trajectory is given by ${\bm{x}}_{\rm
cl}(s)={\bm{x}}_{0}+\frac{s-t_{0}}{t-t_{0}}({\bm{x}}-{\bm{x}}_{0})$. This
corresponds to the product of a Wilson line, trailing the direction of the
particle, times a vacuum propagator, see Eq. (6). Corrections to this limit
can also be systematically be calculated Altinoluk_2014 . In the mixed
representation, this leads to,

$({\bm{p}}|{\cal
G}_{R}(t,t_{0})|{\bm{p}}_{0})\simeq(2\pi)^{2}\delta(\bm{p}-\bm{p}_{0})V_{R}\left(t,t_{0};\left[\bm{x}_{\mathrm{cl}}(s)=\bm{n}s\right]\right)\mathrm{e}^{-i\frac{{\bm{p}}^{2}}{2E}\left(t-t_{0}\right)}\,.$
(13)

where ${\bm{n}}={\bm{p}}/E$, see Mehtar_Tani_2018 ; Dom_nguez_2020 . The last
term in this product is simply the Fourier transform of the vacuum propagator.

In detail, the $1\to 2$ partonic processes we consider are: 1) $\gamma\to
q+\bar{q}$, 2) $q\to q+g$, 3) $g\to g+g$. These will, at most, involve
correlators of 4, 6 and 8 Wilson lines (in the fundamental representation). We
also write out the relevant correlators for $g\to q+\bar{q}$, but we do not
explicitly evaluate the spectrum in this case. All three processes consist of
one (off-shell) particle111We will however only consider physical
polarizations/spin states for the initial particle, since other contributions
do not propagate. traversing the medium splitting into two particles. While we
derive formulas for a generic medium profile, our numerical calculations apply
to a medium with constant density (aka the “brick”), where the splitting can
occur either inside the medium or outside.

As mentioned above, the first particle, with LC energy $E$, is produced at
initial time $t_{0}=0$ and is propagating along the light-cone in the positive
$z$ direction. It splits at times $t_{1}$ in the amplitude and $t_{2}$ in the
complex conjugate amplitude, see Fig. 1. The two daughter particles, which now
carry LC energies $(1-z)E$ and $zE$, respectively, then propagate on the
classical paths ${\bm{r}}_{1}(t)$ (${\bm{r}}_{\bar{1}}(t)$) and
${\bm{r}}_{2}(t)$ (${\bm{r}}_{\bar{2}}(t)$) in the amplitude (complex
conjugate amplitude) to the end of the medium at $L$ . In the high-energy,
eikonal approximation these paths are classical and are given by

$\begin{split}{\bm{r}}_{1}(t)&=\bm{n}_{1}(t-t_{1})\,,\\\
{\bm{r}}_{2}(t)&=\bm{n}_{2}(t-t_{1})\,,\\\
{\bm{r}}_{\bar{1}}(t)&=\bm{n}_{1}(t-t_{2})\,,\\\
{\bm{r}}_{\bar{2}}(t)&=\bm{n}_{2}(t-t_{2})\,,\end{split}$ (14)

where $\bm{n}_{1}\equiv\frac{\bm{p}_{1}}{zE}$ and
$\bm{n}_{2}\equiv\frac{\bm{p}_{2}}{(1-z)E}$ are the transverse velocity
vectors. To slightly compress the notation we will usually refer to the
coordinates as numbers, meaning that we will write $V({\bm{r}}_{1})\equiv
V_{1}$ and $\gamma({\bm{r}}_{1}-{\bm{r}}_{\bar{2}})\equiv\gamma_{1\bar{2}}$,
etc.

Finally, in the harmonic approximation, we need the square of the differences
of the transverse coordinates. Using the eikonal approximation this is

$\begin{split}({\bm{r}}_{1}-{\bm{r}}_{2})^{2}&=(t-t_{1})^{2}\theta^{2}\,,\\\
({\bm{r}}_{\bar{1}}-{\bm{r}}_{\bar{2}})^{2}&=(t-t_{2})^{2}\theta^{2}\,,\\\
({\bm{r}}_{1}-{\bm{r}}_{\bar{1}})^{2}&=(1-z)^{2}(t_{2}-t_{1})^{2}\theta^{2}\,,\\\
({\bm{r}}_{2}-{\bm{r}}_{\bar{2}})^{2}&=z^{2}(t_{2}-t_{1})^{2}\theta^{2}\,,\\\
({\bm{r}}_{1}-{\bm{r}}_{\bar{2}})^{2}&=(t-(1-z)t_{1}-zt_{2})^{2}\theta^{2}\,,\\\
({\bm{r}}_{\bar{1}}-{\bm{r}}_{2})^{2}&=(t-zt_{1}-(1-z)t_{2})^{2}\theta^{2}\,,\end{split}$
(15)

where we have assumed that the angle $\theta$ is small.

##  3 Emission spectra

In this section we will present the results for the in-medium emission spectra
$\frac{{\rm d}I}{{\rm d}z{\rm d}\theta}$ for the in-medium splitting
processes. We refer to appendix A for the details of the calculations. All of
the Wilson line correlators in this section were calculated using the methods
developed in Sec. 4. For more details about the calculation of correlators of
six and eight Wilson lines we refer to appendix B.

One can define the vacuum spectrum as

$\frac{\mathrm{d}I^{\textrm{vac}}}{\mathrm{d}z\,\mathrm{d}\theta}=\frac{\alpha}{\pi}\frac{P(z)}{\theta}\,,$
(16)

where $\alpha$ can be $\alpha_{\textrm{em}}$ or $\alpha_{s}$ depending on the
process, and $P(z)$ is the relevant Altarelli-Parisi splitting function. Then
one can write the full spectrum on the form Dom_nguez_2020

$\displaystyle\frac{\mathrm{d}I^{\textrm{full}}}{\mathrm{d}z\,\mathrm{d}\theta}$
$\displaystyle=\frac{\mathrm{d}I^{\textrm{vac}}}{\mathrm{d}z\,\mathrm{d}\theta}+\frac{\mathrm{d}I^{\textrm{med}}}{\mathrm{d}z\,\mathrm{d}\theta}$
$\displaystyle=\frac{\mathrm{d}I^{\textrm{vac}}}{\mathrm{d}z\,\mathrm{d}\theta}\left(1+F_{\textrm{med}}(z,\theta)\right)\,.$
(17)

The term $F_{\textrm{med}}(z,\theta)$ contains the medium modification to the
processes. For a generic medium profile the medium radiation reads

$\frac{{\rm d}I^{\rm med}}{{\rm d}z\,{\rm d}\theta}=\frac{{\rm d}I^{\rm
vac}}{{\rm d}z\,{\rm d}\theta}\,2{\rm Re}\int_{0}^{L}\frac{{\rm
d}t_{1}}{t_{\rm f}}\int_{t_{1}}^{L}\frac{{\rm d}t_{2}}{t_{\rm f}}\,{\rm
e}^{-i\frac{t_{2}-t_{1}}{t_{\rm
f}}}\mathcal{C}^{(4)}(L,t_{2})\mathcal{C}^{(3)}(t_{2},t_{1})\,,$ (18)

where, in the high-energy limit employed in this paper, the medium-induced
spectrum is proportional to the vacuum spectrum. This proportionality does not
_a priori_ hold in all the phase space, in particular whenever the transverse
momentum in the splitting $k_{\perp}=z(1-z)E\theta$ is comparable to the
transverse momentum accumulated in the medium $Q_{s}\sim\hat{q}L$ Blaizot_2013
. Finally, the factors $\mathcal{C}^{(n)}(t_{b},t_{a})$ appearing in (18) are
$n$-particle correlators that have support during time $t_{a}<t<t_{b}$, and
$t_{\rm f}=\frac{2}{z(1-z)E\theta^{2}}$ is the formation time of the process.
The splitting process is illustrated in Fig. 1.

Figure 1: The process of a photon splitting to a quark-antiquark pair. The
amplitude is on the top and the complex conjugate amplitude on the bottom. The
splitting happens at time $t_{1}$ in the amplitude and at a later time $t_{2}$
in the complex conjugate amplitude.

For a medium with fixed density and extension $L$, we have
$\hat{q}(t)=\hat{q}\,\Theta(L-t)$. In this case, the integrals over the
emission times $t_{1}$ and $t_{2}$ in (18) can be split, so that

$\displaystyle\frac{{\rm d}I^{\rm in-in}}{{\rm d}z\,{\rm d}\theta}$
$\displaystyle=\frac{{\rm d}I^{\rm vac}}{{\rm d}z\,{\rm d}\theta}\,2{\rm
Re}\int_{0}^{L}\frac{{\rm d}t_{1}}{t_{\rm f}}\int_{t_{1}}^{L}\frac{{\rm
d}t_{2}}{t_{\rm f}}\,{\rm e}^{-i\frac{t_{2}-t_{1}}{t_{\rm
f}}}\mathcal{C}^{(4)}(L,t_{2})\mathcal{C}^{(3)}(t_{2},t_{1})\,,$ (19)
$\displaystyle\frac{{\rm d}I^{\rm in-out}}{{\rm d}z\,{\rm d}\theta}$
$\displaystyle=\frac{{\rm d}I^{\rm vac}}{{\rm d}z\,{\rm d}\theta}\,2{\rm
Im}\int_{0}^{L}\frac{{\rm d}t_{1}}{t_{\rm f}}\,{\rm
e}^{-i\frac{L-t_{1}}{t_{\rm f}}}\mathcal{C}^{(3)}(L,t_{1})\,,$ (20)

where ${\rm d}N^{\rm med}/({\rm d}z{\rm d}\theta)={\rm d}N^{\rm in-in}/({\rm
d}z{\rm d}\theta)+{\rm d}N^{\rm in-out}/({\rm d}z{\rm d}\theta)$. Taking into
account that the Wilson line correlators are real the medium modification term
can be written Dom_nguez_2020

$F_{\textrm{med}}=2\int_{0}^{L}\frac{{\rm d}t_{1}}{t_{\rm
f}}\left[\int_{t_{1}}^{L}\frac{{\rm d}t_{2}}{t_{\rm
f}}\cos\left(\frac{t_{2}-t_{1}}{t_{\rm
f}}\right)\mathcal{C}^{(4)}(L,t_{2})\mathcal{C}^{(3)}(t_{2},t_{1})-\sin\left(\frac{L-t_{1}}{t_{\rm
f}}\right)\mathcal{C}^{(3)}(L,t_{1})\right]\,.$ (21)

We now will discuss three concrete cases that are relevant for jet quenching
phenomenology. We will compute the double-differential spectrum for a wide
range of LC energy sharing fraction $z$ and angles $\theta$ to map the regions
where medium-induced corrections appear, as quantified by the factor $F_{\rm
med}(z,\theta)$. Our focus here is to provide a test bed for evaluating
precisely the multi-Wilson line correlators appearing in (21), and we will
therefore not worry about the validity of the eikonal approximation (12) of
the splitting products. Including non-eikonal corrections on the particle
trajectories will be postponed to future work.

###  3.1 Derivation of the splitting functions

##### Photon splitting

We will start with the case of a photon splitting into a quark-antiquark pair,
i.e. $\gamma\to q+\bar{q}$. Due to the least number of fundamental Wilson
lines, this is the simplest process to analyze. We will therefore treat it in
more detail, taking the advantage to discuss the relevant medium and jet
scales appearing in the calculation.

In this case, the vacuum emission spectrum is given by (16) with the QED
coupling constant $\alpha_{\rm em}$ and the Altarelli-Parisi splitting
function being $P_{q\gamma}(z)=n_{f}N_{c}[z^{2}+(1-z)^{2}]$, where $n_{f}$ is
the number of active flavors. Furthermore, the correlator $\mathcal{C}^{(3)}$
reduces to an effective two-point function because the photon does not carry
color charge. We have

$\displaystyle\mathcal{C}_{q\gamma}^{(4)}(L,t_{2})$
$\displaystyle=\frac{1}{N_{c}}\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle\,,$
(22) $\displaystyle\mathcal{C}_{q\gamma}^{(3)}(t_{2},t_{1})$
$\displaystyle=\frac{1}{N_{c}}\operatorname{tr}\langle
V_{1}V_{2}^{\dagger}\rangle\,,$ (23)

where the time extension of each of the medium-averaged color correlators on
the right hand side is implied by the time argument on left hand side of the
equation.

The correlator of two Wilson lines, which in this case corresponds to
$\mathcal{C}_{q\gamma}^{(3)}(t,t_{1})=\mathcal{S}_{12}(t,t_{1})$, is generally
referred to as a dipole correlator, and is known to be

$\mathcal{S}_{12}(t,t_{1})\equiv\frac{1}{N_{c}}\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\rangle={\rm
e}^{-C_{F}\int_{t_{1}}^{t}{\rm d}s\,n(s)\sigma({\bm{r}})}\,,$ (24)

where ${\bm{r}}={\bm{r}}_{1}-{\bm{r}}_{2}$ is the difference of transverse
positions of the two Wilson lines. For a fixed separation, i.e. ${\bm{r}}=$
const., in the HO approximation and in a medium with constant density, it
simply reads $\mathcal{S}_{12}(t,t_{1})={\rm
e}^{-\frac{1}{4}\hat{q}(t-t_{1})\,{\bm{r}}^{2}}$, where, as a reminder, we
have denoted $\hat{q}\equiv\hat{q}_{F}$. However, for the kinematics we
consider, see Eq. (15), this becomes

$\mathcal{S}_{12}(t,t_{1})={\rm
e}^{-\frac{1}{12}\hat{q}(t-t_{1})^{3}\theta^{2}}\,.$ (25)

Then, assuming that $t=t_{2}$ and $t_{2}-t_{1}\sim t_{\rm f}$, we find

$\mathcal{S}_{12}\approx{\rm
e}^{-\frac{2}{3}\frac{\hat{q}}{\omega^{3}\theta^{4}}}\,,$ (26)

with $\omega=z(1-z)E$. This implies that medium modifications appear, in this
term, whenever $\omega^{3}\theta^{4}\lesssim\hat{q}$.

The correlator of four Wilson lines $\mathcal{C}^{(4)}(t,t_{2})$, referred to
as the quadrupole (in the fundamental representation), can only be calculated
numerically at finite-$N_{c}$. We will later show how this can be achieved
through the differential equation Eq. (50). In the large-$N_{c}$ limit,
however, it can be calculated analytically through the simplified differential
equation Eq. (52). There are only two ways of connecting the Wilson lines at
the final time (their connection at initial time is given by the vacuum
splitting process). We can therefore define
$C_{1\bar{2}}(t,t_{2})\equiv\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle$
and
$C_{\bar{2}1}(t,t_{2})\equiv\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle$,
note the absence of explicit normalization factors at this stage. In the
large-$N_{c}$ approximation, the first of these correlators reads simply,

$\frac{1}{N_{c}^{2}}C_{1\bar{2}}(t,t_{2})\simeq{\rm
e}^{-\frac{1}{12}\hat{q}\theta^{2}\left[(t-t_{2})^{3}+(t-t_{1})^{3}-\tau^{3}\right]}\,,$
(27)

where $\tau\equiv t_{2}-t_{1}$, and the dependence on $t_{1}$ appears as a
consequence of the fixed trajectories. Here we have used the eikonal (14) and
harmonic oscillator approximations (10). Similarly, at large-$N_{c}$, the
second correlator is

$\displaystyle\frac{1}{N_{c}}C_{\bar{2}1}(t,t_{2})$ $\displaystyle\simeq{\rm
e}^{-\frac{1}{4}\hat{q}\theta^{2}\xi(t-t_{2})\tau^{2}}$
$\displaystyle-\frac{1}{2}\hat{q}\theta^{2}z(1-z)\tau^{2}\int^{t}_{t_{2}}\mathrm{d}s\,e^{-\frac{1}{4}\hat{q}\theta^{2}\xi(t-s)\tau^{2}}e^{-\frac{1}{12}\hat{q}\theta^{2}\left[(s-t_{2})^{3}+(s-t_{1})^{3}-\tau^{3}\right]}\,,$
(28)

where we defined $\xi\equiv z^{2}+(1-z)^{2}$. Only the latter of these
correlators appears in the spectrum, cf. Eq. (22), but we include both for
completeness.222In Dom_nguez_2020  $C_{\bar{2}1}$ is also calculated in the
large-$N_{c}$ limit. In their Eq. (29) they get the same as (3.1), except they
lack the factor of $1/2$ in front of the second term. Assuming the dominance
of the first term in (3.1), setting $t=L$ and assuming that $L\gg t_{\rm f}$,
we find that

$\frac{1}{N_{c}}C_{\bar{2}1}\approx{\rm
e}^{-\frac{1}{6}\frac{\hat{q}L}{(\omega\theta)^{2}}}\,,$ (29)

where we put $\xi\approx 2/3$. The factor in the exponential becomes large
whenever $\omega\theta<\sqrt{\hat{q}L}$. This factor is related to momentum
broadening of the quark and anti-quark after they have been produced.

Let us compare the two conditions when exponential suppression arise either in
the dipole $\mathcal{S}_{12}(t_{2},t_{1})$ or quadrupole
$C_{\bar{2}1}(L,t_{2})$. For a fixed energy $\omega$, the two conditions are
equal at the critical angle

$\theta_{c}\sim\left(\hat{q}L^{3}\right)^{-1/2}\,.$ (30)

Let us also define the characteristic energies $\omega_{\rm
d}=(\hat{q}/\theta^{4})^{1/3}$ and $\omega_{\rm
broad}=\sqrt{\hat{q}L}/\theta$. At large angles $\theta>\theta_{c}$, the
condition from the dipole starts affecting soft gluon emissions, i.e.
$\omega_{\rm d}<\omega_{\rm broad}$. This reflects the length-dependence color
coherence. On the one hand, the dipole, which has support only during the
formation time $t_{\rm f}\lesssim L$, needs a large angle to resolve the two
particles within that time scale. On the other hand, the quadrupole, which
extends up to $L$, will ultimately resolve even narrower configurations.

We also plot the dependence on the latest time of both
$\mathcal{C}_{q\gamma}^{(3)}(t,t_{1})$ and
$\mathcal{C}_{q\gamma}^{(4)}(t,t_{2})$ in Fig. 2(a), keeping $t_{1}=0.3$ fm
fixed, in the case of the dipole, and both $t_{2}=1$ fm and $t_{1}=0.3$ fm
fixed, in the case of the quadrupole. The other parameters are chosen as
$\hat{q}=1.5$ GeV${}^{2}/$fm, $\theta=0.5$ and $z=0.5$. We notice the fast
decay of the dipole, that goes like $\sim{\rm e}^{-t^{3}}$ according to (25),
compared to the exponential decay of the quadrupole, i.e. $\sim{\rm e}^{-t}$,
at large times. Finally, we notice that the large-$N_{c}$ approximation to the
full quadrupole, given in Eq. (3.1), is very good up very late times.

##### Quark-gluon splitting

Next we consider the slightly more complicated problem of a quark-gluon
splitting. This was also outlined in Dom_nguez_2020 , but not calculated
explicitly. For this process, the vacuum emission spectrum is given by (16),
with the QCD coupling constant $\alpha_{s}$ and the Altarelli-Parisi splitting
function $P_{gq}(z)=C_{F}\frac{1+(1-z)^{2}}{z}$. The four- and three-point
functions read

$\displaystyle\mathcal{C}_{gq}^{(4)}(L,t_{2})$
$\displaystyle=\frac{1}{N_{c}^{2}-1}\left\langle\operatorname{tr}[V_{\bar{1}}^{\dagger}V_{1}V_{2}^{\dagger}V_{\bar{2}}]\operatorname{tr}[V_{\bar{2}}^{\dagger}V_{2}]-\frac{1}{N_{c}}\operatorname{tr}[V_{\bar{1}}^{\dagger}V_{1}]\right\rangle\,,$
(31) $\displaystyle\mathcal{C}_{gq}^{(3)}(t_{2},t_{1})$
$\displaystyle=\frac{1}{N_{c}^{2}-1}\left\langle\operatorname{tr}[V_{2}^{\dagger}V_{1}]\operatorname{tr}[V_{0}^{\dagger}V_{2}]-\frac{1}{N_{c}}\operatorname{tr}[V_{0}^{\dagger}V_{1}]\right\rangle\,.$
(32)

The emission spectrum is composed of correlators of two, four and six Wilson
lines. The three-point function can be solved exactly, see (B.1), resulting in

$\displaystyle\mathcal{C}_{gq}^{(3)}(t_{2},t_{1})$ $\displaystyle={\rm
e}^{-\frac{1}{2}\int_{t_{1}}^{t_{2}}{\rm
d}s\,n(s)[N_{c}(\sigma_{02}+\sigma_{12})-\frac{1}{N_{c}}\sigma_{01}]}$
$\displaystyle={\rm
e}^{-\frac{1}{12}{\hat{q}}(t_{2}-t_{1})^{3}\theta^{2}\left(1+z^{2}+\frac{2z}{N_{c}^{2}-1}\right)}\,.$
(33)

This expression is very similar to the dipole term in Eq. (25) and the same
scale analysis applies.

The four-point correlator involving six and two Wilson lines can only be
calculated numerically at finite $N_{c}$. In the large-$N_{c}$ limit, the
former can be calculated analytically, and reads

$\displaystyle\frac{1}{N_{c}^{2}}\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\operatorname{tr}[V_{2}V_{\bar{2}}^{\dagger}]\rangle\simeq{\rm
e}^{-\frac{1}{4}{\hat{q}}\theta^{2}(t-t_{2})(t_{2}-t_{1})^{2}(1-2z+3z^{2})}$
$\displaystyle\times\left(1-\frac{1}{2}{\hat{q}}\theta^{2}z(1-z)(t_{2}-t_{1})^{2}\int_{t_{2}}^{t}{\rm
d}s\,{\rm
e}^{-\frac{1}{12}{\hat{q}}\theta^{2}\left[(s-t_{2})^{2}(2s-3t_{1}+t_{2})+6z(1-z)(s-t_{2})(t_{2}-t_{1})^{2}\right]}\right)\,.$
(34)

Once again, the first term in the correlator above has a form very similar to
the four-point function relevant for photon splitting, see Eq. (3.1).

##### Gluon-gluon splitting

The last process of interest is the case of a gluon splitting into two other
gluons. This process was discussed quite extensively in Blaizot_2013 . For
this process, the vacuum emission spectrum is given by (16) with the QCD
coupling constant $\alpha_{s}$ and the Altarelli-Parisi splitting function
$P_{gg}(z)=2N_{c}\big{[}z(1-z)+\frac{1-z}{z}+\frac{z}{1-z}\big{]}$. In this
case the 4- and 3-point functions read

$\displaystyle\mathcal{C}_{gg}^{(4)}(L,t_{2})$
$\displaystyle=\frac{1}{N_{c}(N_{c}^{2}-1)}\left\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}]\operatorname{tr}[V_{2}V_{\bar{2}}^{\dagger}V_{\bar{1}}V_{1}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{2}^{\dagger}]-\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}V_{2}V_{\bar{2}}^{\dagger}V_{\bar{1}}V_{1}^{\dagger}V_{\bar{2}}V_{2}^{\dagger}]\right\rangle\,,$
(35) $\displaystyle\mathcal{C}_{gg}^{(3)}(t_{2},t_{1})$
$\displaystyle=\frac{1}{N_{c}(N_{c}^{2}-1)}\left\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{0}V_{1}^{\dagger}]\operatorname{tr}[V_{2}V_{0}^{\dagger}]-\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{0}V_{1}^{\dagger}V_{2}V_{0}^{\dagger}]\right\rangle\,.$
(36)

When cast as correlators of Wilson lines in the fundamental representation,
the $\mathcal{C}^{(4)}_{gg}$ involves 8-point correlators, which is the
largest number we will calculate in detail.

The 3-point function can be solved exactly, either by the differential
equation (4.2) or by writing it in terms of adjoint Wilson lines (93). In the
end, the result reads

$\displaystyle\mathcal{C}_{gg}^{(3)}(t_{2},t_{1})$ $\displaystyle={\rm
e}^{-\frac{N_{c}}{2}\int_{t_{1}}^{t_{2}}{\rm
d}t\,n(t)[\sigma_{01}+\sigma_{02}+\sigma_{12}]}$ $\displaystyle={\rm
e}^{-\frac{1}{12}{\hat{q}}(t_{2}-t_{1})^{3}\theta^{2}\frac{N_{c}}{C_{F}}(1-z+z^{2})}\,.$
(37)

Note the similarity to the previous results, see Eqs. (25) and (3.1).

The 4-point function consists of two different correlators of eight Wilson
lines. They can be calculated through the differential equation in Eq. (4.2).
Interestingly, the four-point function $\mathcal{C}^{(4)}_{gg}$ involves a
eight-point correlator, see the second term in (35), which cannot be reduced
further in the large-$N_{c}$ approximation. This can nevertheless still be
exactly solved in the large-$N_{c}$ approximation, which we present in the
figures below, but the expression is too long to extract any meaningful
approximation. Anticipating the numerical results, we can mention that it is
for this correlator that the large-$N_{c}$ approximation gives the biggest
deviations with respect to the exact result.

##### Gluon-quark splitting

We now consider a gluon that splits into a quark-antiquark pair. The
Altarelli-Parisi splitting function is $P_{qg}(z)=n_{f}T_{R}[z^{2}+(1-z)^{2}]$
and the correlators read

$\displaystyle\mathcal{C}_{qg}^{(4)}(L,t_{2})$
$\displaystyle=\frac{1}{N_{c}}\left\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]-\frac{1}{N_{c}}\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{2}^{\dagger}]\right\rangle\,,$
(38) $\displaystyle\mathcal{C}_{qg}^{(3)}(t_{2},t_{1})$
$\displaystyle=\frac{1}{N_{c}^{2}-1}\left\langle\operatorname{tr}[V_{1}V_{0}^{\dagger}]\operatorname{tr}[V_{0}V_{2}^{\dagger}]-\frac{1}{N_{c}}\operatorname{tr}[V_{1}V_{2}^{\dagger}]\right\rangle\,.$
(39)

Since these expressions involve only quadrupoles and dipoles, that were
previously encountered and analyzed in detail above, we will not present
further results for this splitting process.

###  3.2 Numerical results

Here we present the numerical calculations of the results from the previous
section. We focus first on the details of the three- and four-point functions
for each of the three splitting processes, and proceed with calculating the
double-differential spectrum in the momentum sharing fraction $z$ and angle
$\theta$.

(a) Photon splitting.

(b) Quark-gluon splitting.

(c) Gluon-gluon splitting.

Figure 2: The time evolution of $\mathcal{C}^{(3)}(L,t_{1})$ and
$\mathcal{C}^{(4)}(L,t_{2})$ for the three processes. For
$\mathcal{C}^{(4)}(L,t_{2})$ both the exact and large-$N_{c}$ versions are
plotted.

In Fig. 2, we show how $\mathcal{C}^{(3)}_{ij}(t,t_{1})$, with blue, solid
curves, and $\mathcal{C}^{(4)}_{ij}(t,t_{2})$, with orange, solid curves, for
the three processes evolve with time. For the four-point functions, we also
plot the large-$N_{c}$ approximation with orange, dashed curves. We fix both
$t_{1}=0.3$ fm and $t_{2}=1$ fm and plot for the latest time $t=L$. For the
other parameters we choose $\hat{q}=1.5$ GeV${}^{2}/$fm, $\theta=0.5$ and
$z=0.5$.

While this approximation turns out work extremely well for the photon
splitting, see Fig. 2(a), we note that it has a more limited range of
applicability for both the quark-gluon, see Fig. 2(b), and gluon-gluon, see
Fig. 2(c), splitting processes, respectively. In all of the cases the exact
value is slightly higher than the approximate one. As we derived analytically,
the $\mathcal{C}^{(3)}_{ij}$ terms all decay as $\sim{\rm
e}^{-{\hat{q}}(t-t_{1})^{3}\tau^{2}\theta^{2}f(z)}$, where $f(z)$ is a process
dependent regular function. The $\mathcal{C}^{(4)}_{ij}$ terms are more
complicated, especially at early times where all terms contribute, but at late
times the dominant contribution comes from $\sim{\rm
e}^{-{\hat{q}}(t-t_{2})\theta^{2}}$.

(a) Photon splitting.

(b) Quark-gluon splitting.

(c) Gluon-gluon splitting.

Figure 3: The medium modification factor $F_{\textrm{med}}(z,\theta)$ for
three splitting processes as a function of $\theta$ and $z$ with $L=2$ fm and
$E=100$ GeV at finite $N_{c}$.

The ratio of double-differential in-medium to vacuum spectrum reveals the
medium modification factor $F_{\rm med}(z,\theta)={\rm d}I^{\rm med}/({\rm
d}z{\rm d}\theta)\big{/}{\rm d}I^{\rm vac}/({\rm d}z{\rm d}\theta)$. We plot
this factor, calculated at finite $N_{c}$, for the three processes in Fig. 3.
These results have been obtained for the medium parameters $\hat{q}=1.5$
GeV2/fm and $L=2$ fm and an energy of the initial particle, before splitting,
of $E=100$ GeV.

As one can see from Fig. 3, the medium modification factor
$F_{\textrm{med}}(z,\theta)$ has roughly the same characteristic shape for all
three processes. The medium modifications appear at large angles
$\theta>\theta_{c}$, in between the characteristic lines
$\omega^{3}\theta^{3}<\hat{q}$ and $\omega^{2}\theta^{2}<\hat{q}L$ which we
have identified for the three- and four-point functions in Sec. 3.1. This
corresponds to formation times smaller than the medium length, $t_{\rm f}<L$.
In fact, we can recast these conditions in terms of the formation time of the
process, namely $t_{\rm f}<t_{\rm d}$ and $t_{\rm f}<t_{\rm broad}$, where

$t_{\rm d}\sim\left(\frac{1}{\hat{q}\theta^{2}}\right)^{1/3}\,,\qquad\text{and
}\qquad t_{\rm broad}\sim\left(\frac{1}{\hat{q}\theta^{2}L}\right)^{1/2}\,.$
(40)

The modifications appear for the range of formation times $t_{\rm
broad}<t_{\rm f}<t_{\rm d}$ and $\theta>\theta_{c}$ Dom_nguez_2020 . There
also seems to be a trend that both the magnitude and the region of the
modifications grow with the number of Wilson lines. This can be traced back to
the finite terms, $f(z)$, in the exponents that modify the scaling behavior.
Naively, we would expect the relevant jet quenching parameter to be roughly a
factor $N_{c}/C_{F}\approx 2$ larger for gluon splitting than for the photon.

(a) Photon splitting.

(b) Quark-gluon splitting.

(c) Gluon-gluon splitting.

Figure 4: The ratio $F_{\rm med}(z,\theta)|_{{\rm large-}N_{c}}/F_{\rm
med}(z,\theta)$ for three splitting processes as a function of $\theta$ and
$z$ with $L=2$ fm and $E=100$ GeV.

Our main focus in this work is to highlight the differences between the
finite-$N_{c}$ results versus their large-$N_{c}$ approximated counterparts.
To illustrate this we have plotted the ratio of the exact and large-$N_{c}$
medium modification factors, i.e. $F_{\rm med}(z,\theta)|_{{\rm
large-}N_{c}}/F_{\rm med}(z,\theta)$ in Fig. 4. The difference between the
exact and approximate result is small in the whole phase space in the photon
splitting case, where there is a correlator of four Wilson lines. However, in
the cases of quark-gluon and especially gluon-gluon splitting, which contain
correlators of six and eight Wilson lines, the error can be relatively big,
maximally of the order of 16% in case of the latter process. This is the
reflection of the behavior observed previously in Fig. 2. From these
calculations it seems like the more complicated color structure, the bigger
the error is by using the large-$N_{c}$ approximation. Once again, the error
becomes most sizable at relatively large in-medium formation times, i.e.
$t_{\rm f}\sim t_{\rm broad}$ and $t_{\rm f}\sim t_{\rm d}$, but at the same
time $t_{\rm f}<L$. This is most clearly seen in the gluon-gluon splitting,
cf. Fig. 4(c). Finally, we note that the finite-$N_{c}$ corrections come as a
modulation along the previously established scaling lines which hints that
such corrections could perhaps be absorbed into an _effective_ jet quenching
parameter.

To summarize, we have calculated the the double-differential spectrum
$\frac{{\rm d}I}{{\rm d}z{\rm d}\theta}$ for three different splitting
processes, and shown that the resulting expressions factorize into three- and
four-point functions that contain medium-averaged products of 2, 4, 6 and 8
fundamental Wilson lines. In the coming Sec. 4 we will detail how these are
calculated. Strikingly, the three- and four-point functions all take a very
similar scaling form as was derived analytically exactly, for the former, and
in the large-$N_{c}$ approximation, for the latter. This corresponds to the
identification of two characteristic time-scales in the medium, related to
broadening along the length of the medium, $t_{\rm broad}$, and decoherence
during the formation of the splitting, $t_{\rm d}$. These were identified
first in Dom_nguez_2020  for the photon splitting process, and we have here
extended their validity to all other splitting QCD processes. Finally, we have
seen that finite-$N_{c}$ corrections play an increasingly important role the
bigger the total color charge involved in the splitting process.

##  4 Calculating Wilson line correlators

In this section we will present our method for calculating Wilson line
correlators. As an illustration we will first show how it is done in the
simple case of four Wilson lines in the fundamental representation. Thereafter
this process will be generalized to an arbitrary number of Wilson lines.

###  4.1 Four Wilson lines

The simplest Wilson line correlator comes from the pair production process
(22), where there is a trace of four Wilson lines
$\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle$.
In this section, we will show how to derive a system of differential equations
to calculate this. Let the Wilson lines have support from $t_{0}$ to some
arbitrary time $t+\epsilon$. Then, following Kovner_2001 , we expand them
between $t$ and $t+\epsilon$ to get

$\displaystyle
V(t+\epsilon,t_{0};\mathbf{r})=V(t+\epsilon,t;\mathbf{r})V(t,t_{0};\mathbf{r})$
$\displaystyle=\left(1+ig\int^{t+\epsilon}_{t}{\rm
d}s\,A^{a}(s,\mathbf{r})t^{a}-\frac{g^{2}}{2!}\int^{t+\epsilon}_{t}{\rm
d}s\int^{t+\epsilon}_{t}{\rm
d}s^{\prime}\,A^{a}(s,\mathbf{r})A^{b}(s^{\prime},\mathbf{r}^{\prime})t^{a}t^{b}+\mathcal{O}(\epsilon^{2})\right)$
$\displaystyle\times V(t,t_{0};\mathbf{r}),$ (41)

where we have kept some of the color indices implicit. All four Wilson lines
are expanded in this manner. We end up with having to take the medium average
of the integrals over two medium fields, traced over the relevant color
indices, which is dealt in the following way

$\displaystyle\int^{t+\epsilon}_{t}{\rm d}s\int^{t+\epsilon}_{t}{\rm
d}s^{\prime}\left\langle\,A^{a}(s,\mathbf{r})A^{b}(s^{\prime},\mathbf{r}^{\prime})\operatorname{tr}[t^{a}t^{b}]\right\rangle$
$\displaystyle=\int^{t+\epsilon}_{t}{\rm
d}s\,n(s)\gamma({\bm{r}}-{\bm{r}}^{\prime})t_{ij}^{a}t_{ji}^{a}$
$\displaystyle\simeq\epsilon\,C_{F}n(t)\gamma({\bm{r}}-{\bm{r}}^{\prime})\,,$
(42)

where in the first step we applied the medium average (7). Then, keeping terms
up to the first order of $\epsilon$ this becomes

$\displaystyle\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle_{(t+\epsilon)}$
$\displaystyle=\left(1+\epsilon
g^{2}n(t)C_{F}[\gamma_{1\bar{1}}+\gamma_{2\bar{2}}-2\gamma_{0}]\right)\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle_{(t)}$
$\displaystyle-\epsilon
g^{2}n(t)[\gamma_{1\bar{2}}-\gamma_{12}-\gamma_{\bar{1}\bar{2}}+\gamma_{2\bar{1}}]\langle\operatorname{tr}[t^{a}V_{1}V_{2}^{\dagger}t^{a}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle_{(t)}\,.$
(43)

Using the Fierz identity

$t_{ij}^{a}t_{kl}^{a}=\frac{1}{2}\left(\delta_{il}\delta_{jk}-\frac{1}{N_{c}}\delta_{ij}\delta_{kl}\right)\,,$
(44)

this results in the differential equation

$\displaystyle\frac{{\rm d}}{{\rm
d}t}\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle_{(t)}=\lim_{\epsilon\to
0}\frac{\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle_{(t+\epsilon)}-\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle_{(t)}}{\epsilon}$
$\displaystyle=g^{2}n(t)\left[C_{F}(\gamma_{1\bar{1}}+\gamma_{2\bar{2}}-2\gamma_{0})+\frac{1}{2N_{c}}(\gamma_{1\bar{2}}-\gamma_{12}-\gamma_{\bar{1}\bar{2}}+\gamma_{2\bar{1}})\right]\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle_{(t)}$
$\displaystyle-\frac{1}{2}g^{2}n(t)(\gamma_{1\bar{2}}-\gamma_{12}-\gamma_{\bar{1}\bar{2}}+\gamma_{2\bar{1}})\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle_{(t)}\,.$
(45)

It is evident that the original term mixes with another four-point correlator,
given in the term on the last line. To understand this, let us look closer at
the term
$\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle$.
Since the process is happening in the medium the quarks and antiquarks can at
any time exchange gluons, so their color is continuously rotating. In the case
of four Wilson lines there are two possible ways of connecting the color at
time $t$ to ensure color conservation, namely as shown in Fig. 5. The second
way is exactly the term
$\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle$
that appeared in equation (4.1). The inclusion of this term in the
differential equation (4.1) just represents the possibility for color rotation
to happen at each time.

Figure 5: The two possible ways of color connecting the four Wilson lines. On
the left is
$C_{\bar{2}1}\equiv\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle$,
while on the right is
$C_{1\bar{2}}\equiv\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle$,
both going from times $t_{2}$ to an arbitrary time $t$. The grey lines at the
beginning and end indicate the colour connections.

To continue one can find a complementary differential equation for
$\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle$
and see if we can find a solution for the set. Going through the same
procedure as above gives

$\displaystyle\frac{{\rm d}}{{\rm
d}t}\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle_{(t)}$
$\displaystyle=g^{2}n(t)\left[C_{F}(\gamma_{12}+\gamma_{\bar{1}\bar{2}}-2\gamma_{0})+\frac{1}{2N_{c}}(\gamma_{1\bar{2}}-\gamma_{1\bar{1}}-\gamma_{2\bar{2}}+\gamma_{\bar{1}2})\right]\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle_{(t)}$
$\displaystyle-\frac{1}{2}g^{2}n(t)(\gamma_{1\bar{2}}-\gamma_{1\bar{1}}-\gamma_{2\bar{2}}+\gamma_{\bar{1}2})\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle_{(t)}\,.$
(46)

We now have a set of two coupled differential equations. To save space the
following notation will be used
$C_{1\bar{2}}(t)\equiv\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle_{(t)}$
and
$C_{\bar{2}1}(t)\equiv\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle_{(t)}$.
This notation warrants some more explanation. Both of these expressions are
composed of the two pairs of Wilson lines, namely $V_{1}V_{2}^{\dagger}$ and
$V_{\bar{2}}V_{\bar{1}}^{\dagger}$. The only difference is how to connect
them. The two subscripts in the $C$’s tell which Wilson line that comes
immediately after the two pairs. So $C_{\bar{2}1}$ means that
$V_{1}V_{2}^{\dagger}$ is connected to $V_{\bar{2}}$ and
$V_{\bar{2}}V_{\bar{1}}^{\dagger}$ connects to $V_{1}$. The result is
$\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle$.
This notation might seem overly complicated, but it will prove to be useful
when considering more than four Wilson lines.

The two differential equations (4.1) and (4.1) can be gathered into the
following system,

$\frac{{\rm d}}{{\rm d}t}\begin{bmatrix}C_{1\bar{2}}(t)\\\
C_{\bar{2}1}(t)\end{bmatrix}=-\frac{n(t)}{2}\mathbb{M}\begin{bmatrix}C_{1\bar{2}}(t)\\\
C_{\bar{2}1}(t)\end{bmatrix}\,,$ (47)

where the evolution matrix takes the following form,

$\mathbb{M}=\begin{bmatrix}2C_{F}(\sigma_{12}+\sigma_{\bar{2}\,\bar{1}})+\frac{1}{N_{c}}\Sigma_{1}&-\Sigma_{1}\\\
-\Sigma_{2}&2C_{F}(\sigma_{1\bar{1}}+\sigma_{\bar{2}2})+\frac{1}{N_{c}}\Sigma_{2}\end{bmatrix}\,.$
(48)

Here we have used Eq. (9) to define
$\sigma_{12}=\sigma({\bm{r}}_{1}-{\bm{r}}_{2})$, and introduced

$\displaystyle\Sigma_{1}$
$\displaystyle\equiv\sigma_{1\bar{2}}+\sigma_{2\bar{1}}-\sigma_{1\bar{1}}-\sigma_{2\bar{2}}$
$\displaystyle\Sigma_{2}$
$\displaystyle\equiv\sigma_{1\bar{2}}+\sigma_{\bar{1}2}-\sigma_{12}-\sigma_{\bar{1}\bar{2}}\,.$
(49)

To proceed, we employ the the harmonic approximation (10). For the eikonal,
straight-line trajectories, given in Eqs. (14), the evolution matrix becomes

$-\frac{n(t)}{2}\mathbb{M}=-\frac{\hat{q}\theta^{2}}{4C_{F}}\begin{bmatrix}C_{F}[(t-t_{1})^{2}+(t-t_{2})^{2}]-\frac{1}{N_{c}}(t-t_{1})(t-t_{2})&-(t-t_{1})(t-t_{2})\\\
z(1-z)\tau^{2}&C_{F}\tau^{2}\xi-\frac{1}{N_{c}}z(1-z)\tau^{2}\end{bmatrix}\,,$
(50)

where we have defined $\tau\equiv t_{2}-t_{1}$, $\xi=z^{2}+(1-z)^{2}$ and
assumed that the angle between the two particles $\theta$ is small.
Unfortunately, since the matrix elements depend on time in our setup, we can
only solve this system of differential equations exactly by using numerical
methods.

The authors of Dom_nguez_2020  calculated the four-point function
$\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle$
in the large-$N_{c}$ limit, which is interesting to compare with our results.
This example is illustrative of the general structure of the hierarchy between
the different correlators, and we will therefore go through it in detail. To
take the large-$N_{c}$ limit you start the system of differential equations
(47) and count the powers of $N_{c}$ in each term in the evolution matrix and
the vector of correlators, taking into account that $C_{1\bar{2}}\sim
N_{c}^{2}$ and $C_{\bar{2}1}\sim N_{c}^{1}$. In this limit we also have
$C_{F}\sim N_{c}/2$. The terms on the right-hand side of (47) then have the
following powers of $N_{c}$,

$\displaystyle\begin{bmatrix}\mathcal{O}(N_{c}^{0})+\mathcal{O}(N_{c}^{-2})&\mathcal{O}(N_{c}^{-1})\\\
\mathcal{O}(N_{c}^{-1})&\mathcal{O}(N_{c}^{0})+\mathcal{O}(N_{c}^{-2})\end{bmatrix}\begin{bmatrix}\mathcal{O}(N_{c}^{2})\\\
\mathcal{O}(N_{c}^{1})\end{bmatrix}$
$\displaystyle\xrightarrow[]{\text{large}-N_{c}}\begin{bmatrix}\mathcal{O}(N_{c}^{0})&0\\\
\mathcal{O}(N_{c}^{-1})&\mathcal{O}(N_{c}^{0})\end{bmatrix}\begin{bmatrix}\mathcal{O}(N_{c}^{2})\\\
\mathcal{O}(N_{c}^{1})\end{bmatrix}\,.$ (51)

The large-$N_{c}$ approximation amounts to dropping all the terms in the
matrix that are not scaling with the same power of $N_{c}$ as the original
vector, given by the second term in (51). We see that the the next-to-leading
power of $N_{c}$ turns out to be a factor $N_{c}^{-2}$ smaller compared to the
leading terms. This scaling has also been corroborated generally for $n$-line
correlators in Sec. 4.2.

Hence, employing the large-$N_{c}$ approximation leads to the simplified
system of equations

$\frac{{\rm d}}{{\rm d}t}\begin{bmatrix}C_{1\bar{2}}(t)\\\
C_{\bar{2}1}(t)\end{bmatrix}\simeq-\frac{\hat{q}\theta^{2}}{4N_{c}}\begin{bmatrix}N_{c}[(t-t_{1})^{2}+(t-t_{2})^{2}]&0\\\
2z(1-z)\tau^{2}&N_{c}\tau^{2}\xi\end{bmatrix}\begin{bmatrix}C_{1\bar{2}}(t)\\\
C_{\bar{2}1}(t)\end{bmatrix}\,.$ (52)

Now it is evident that the differential equation for $C_{1\bar{2}}$ is
separable and can be solved easily, which means that $C_{\bar{2}1}$ also can
be solved. This leads to the equations (27) and (3.1). The physical picture of
this differential equation is quite transparent. The correlator of the two
particles (described by two lines in the amplitude and two lines in the
complex conjugate amplitude) can be in either of the states shown in Fig. 5,
and there is a possibility of exchanging a gluon and transferring from one
state to the other. This is encoded in the off-diagonal terms in the matrix
(48), and is associated with a factor of $\sim\sigma$, which scales as
$N_{c}^{-1}$. Say you start in the state $C_{1\bar{2}}$ shown on the right in
figure 5, scaling as $N_{c}^{2}$. If you exchange a gluon you pick up a factor
$N_{c}^{-1}$ from the $\sigma$, and go to the state $C_{\bar{2}1}$, which is a
single trace correlator that scales as $N_{c}^{1}$, so in total this
transition is associated with a factor $N_{c}^{0}$. This is a factor
$N_{c}^{-2}$ smaller compared to the starting point so it can safely be
dropped in the large-$N_{c}$ limit. However, starting with $C_{\bar{2}1}$ and
going to $C_{1\bar{2}}$ you go from a state that scales as $N_{c}^{1}$ to one
scaling as $N_{c}^{2}$, but you lose a power of $N_{c}$ from the $\sigma$, so
in this case the end result has the same $N_{c}$ scaling as the starting
point. Hence, in the large-$N_{c}$ limit you can drop the upper right term in
the matrix, but must keep the lower left one, see Sec. 4.2 for a general
argument for $n$-point correlators.

Figure 6: The exact and large-$N_{c}$ version of
$C_{1\bar{2}}(t)=\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle$
(blue, solid and blue, dashed lines, respectively) and
$C_{\bar{2}1}(t)=\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle$
(orange, solid and orange, dashed lines, respectively). We also plot only the
leading, diagonal term of the large-$N_{c}$ approximation of $C_{\bar{2}1}$
(orange, dotted line) which exhibits the correct large-time asymptotic
behavior.

The solutions to (50) and (52) were plotted in Fig. 6 (solid and dashed lines,
respectively). For this particular case, the agreement between the
large-$N_{c}$ approximation and the exact, finite-$N_{c}$ result is strikingly
good for the $C_{\bar{2}1}$ correlator. At late times, we observe an
exponential suppression, $\propto{\rm e}^{-t}$, with a slope that is in good
agreement with the first term of Eq. (3.1). At early times, there is an
interplay between $C_{\bar{2}1}$ and $C_{1\bar{2}}$ that leads to a more rapid
decrease initially. This is however well captured by the large-$N_{c}$
approximation, given by both terms in Eq. (3.1).

The $C_{1\bar{2}}$ correlator is described well within the large-$N_{c}$
approximation at early times. However, at late times it exhibits a long tail
that is not captured within this approximation. This can be remedied by
including sub-leading corrections in color.

Sub-leading corrections in color can be incorporated to improve on the
sometimes crude large-$N_{c}$ calculation above. To do this write the full
correlators as the sum of their large-$N_{c}$ versions calculated through (52)
and some smaller correction term,

${\bm{C}}={\bm{C}}^{(0)}+{\bm{C}}^{(1)}\,,$ (53)

where ${\bm{C}}=\big{(}C_{1\bar{2}},C_{\bar{2}1}\big{)}^{\intercal}$ is a
vector of the correlators in question, so that ${\bm{C}}^{(1)}$ is a factor
$\mathcal{O}(N_{c}^{-2})$ smaller than ${\bm{C}}^{(0)}$. We can also write the
matrix $\mathbb{M}$ in a form that isolates the large-$N_{c}$ terms from the
finite-$N_{c}$ corrections, i.e.

$\mathbb{M}=\mathbb{M}^{(0)}+\mathbb{M}^{\rm corr.}\,,$ (54)

where the first term strictly corresponds to the leading terms in the
large-$N_{c}$ limit. In our example above, we find that

$-\frac{n(t)}{2}\mathbb{M}^{(0)}=-\frac{\hat{q}\theta^{2}}{4}\begin{bmatrix}(t-t_{1})^{2}+(t-t_{2})^{2}&0\\\
\frac{2}{N_{c}}z(1-z)\tau^{2}&\tau^{2}\xi\end{bmatrix}\,,$ (55)

while

$-\frac{n(t)}{2}\mathbb{M}^{\rm
corr.}\simeq\frac{\hat{q}\theta^{2}}{2N_{c}^{2}}\begin{bmatrix}(t-t_{1})(t-t_{2})&N_{c}(t-t_{1})(t-t_{2})\\\
-\frac{1}{N_{c}}z(1-z)\tau^{2}&z(1-z)\tau^{2}\end{bmatrix}\,,$ (56)

where we expanded the correction matrix to find the leading terms in $N_{c}$.
It can be confirmed that the overall correction to both correlators is of the
order $N_{c}^{-2}$.

The correlators at leading color, i.e. ${\bm{C}}^{(0)}$, are known. They solve
the simplified set of equations ${\rm d}{\bm{C}}^{(0)}(t)/{\rm
d}t=-\frac{n(t)}{2}\mathbb{M}^{(0)}{\bm{C}}^{(0)}(t)$, and are given
explicitly in (27) and (3.1). This can now be used to calculate the color sub-
leading contributions ${\bm{C}}^{(1)}$. Simply plugging this into the full
differential equation (47) results in the following differential equation for
the first correction

$\displaystyle\frac{{\rm d}}{{\rm d}t}{\bm{C}}^{(1)}(t)\simeq$
$\displaystyle-\frac{n(t)}{2}\mathbb{M}^{(0)}{\bm{C}}^{(1)}(t)-\frac{n(t)}{2}\mathbb{M}^{\rm
corr}{\bm{C}}^{(0)}(t)\,,$ (57)

where we have neglected terms that are even more sub-leading, i.e. resulting
from $\mathbb{M}^{\rm corr.}{\bm{C}}^{(1)}$. This is an nonhomogeneous version
of the large-$N_{c}$ system of differential equations (52), and can also be
solved exactly. As an example the first correction to $C^{1}_{1\bar{2}}(t)$ is

$\displaystyle C^{(1)}_{1\bar{2}}(t)$
$\displaystyle=\frac{{\hat{q}}\theta^{2}}{2N_{c}^{2}}\int_{t_{2}}^{t}{\rm
d}s\,(s-t_{1})(s-t_{2})\left[C^{(0)}_{1\bar{2}}(s)+N_{c}C^{(0)}_{\bar{2}1}(s)\right]$
$\displaystyle\times{\rm
e}^{\frac{{\hat{q}}\theta^{2}}{12}[(t-t_{1})^{3}-(s-t_{1})^{3}+(t-t_{2})^{3}-(s-t_{2})^{3}]}\,.$
(58)

The first correction contains $C^{(0)}_{\bar{2}1}(t)$, given in (3.1) which as
can be seen in Fig. 6 has a linear tail at long times. One would therefore
expect that this correction will rectify the difference between the exact
calculation and the large-$N_{c}$ version of $C_{1\bar{2}}(t)$ at long times
which can be seen in the same plot. On Fig. 6, we have plotted this
correction, and it is indeed clear that it contains this linear tail. It is
also worth noticing that the $N_{c}$-scaling of the correction is
$C^{(1)}_{1\bar{2}}\sim N_{c}^{0}$, since there is an $N_{c}^{-1}$ in the pre-
factor and $C^{(0)}_{\bar{2}1}\sim N_{c}^{1}$. As expected the correction is
lower by a factor $N_{c}^{-2}$ compared to the large-$N_{c}$ result.

It is possible to calculate higher order corrections going as $N_{c}^{-4}$,
$N_{c}^{-6}$ etc. compared to the large-$N_{c}$ expression using the same
technique recursively.

###  4.2 General method for Wilson line correlators

In Sec. 3, we showed that doing similar calculations starting with a quark or
a gluon emitting a gluon leads to correlators of six and eight fundamental
Wilson lines, respectively. We will now generalize the procedure demonstrated
in the preceding section and develop a method of calculating correlators of an
arbitrary number of fundamental Wilson lines. To be more precise we get
systems of differential equations like in (47), and will show how to easily
calculate all the matrix elements in the $K!\times K!$ matrices. The system
can then be solved numerically or, as we will see, analytically in the
large-$N_{c}$ limit.

The correlators of six and eight Wilson lines that appeared in Sec. 3 are

  * •

$\langle\operatorname{tr}[V_{\bar{1}}^{\dagger}V_{1}V_{2}^{\dagger}V_{\bar{2}}]\operatorname{tr}[V_{\bar{2}}^{\dagger}V_{2}]\rangle$,

  * •

$\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}]\operatorname{tr}[V_{2}V_{\bar{2}}^{\dagger}V_{\bar{1}}V_{1}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{2}^{\dagger}]\rangle$,

  * •

and
$\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}V_{2}V_{\bar{2}}^{\dagger}V_{\bar{1}}V_{1}^{\dagger}V_{\bar{2}}V_{2}^{\dagger}]\rangle$.

Note that one can divide the correlators of these Wilson lines into pairs on
the form $[V_{n}V_{\bar{m}}^{\dagger}]_{i_{n}j_{m}}$ times some Kronecker
deltas that connect the indices.

To start, consider the special case of calculating a correlator involving $K$
pairs of a Wilson line in the amplitude times the same Wilson line in the
complex conjugate amplitude

$\langle[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{1}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{2}}\dots[V_{K}V_{\bar{K}}^{\dagger}]_{i_{K}j_{K}}\rangle=\langle\prod_{n=1}^{K}[V_{n}V_{\bar{n}}^{\dagger}]_{i_{n}j_{n}}\rangle\,.$
(59)

This is very useful to consider, even though none of the correlators mentioned
above are of this exact form. The reason is that in this form all of the
formulas derived in this section become much nicer. In addition, it is easy to
generalize this to include all cases simply by changing the labels of the
Wilson lines in (59) to whatever is needed in the specific problem at hand.
For example, choosing $K=3$ and changing labels
$(1,\bar{1},2,\bar{2},3,\bar{3})\to(1,2,\bar{2},\bar{1},2,\bar{2})$ gives the
structure needed in (31), while $K=4$ and changing labels
$(1,\bar{1},2,\bar{2},3,\bar{3},4,\bar{4})\to(1,\bar{1},2,\bar{2},\bar{1},1,\bar{2},2)$
reproduces the correlators in (35). So even though it seems we are calculating
a special case, simply changing the labels in the equations in this section
leads to all possible cases.

To compress the notation a bit we will write the k’th instance of a Wilson
line pair as

$W^{k}_{i_{k}j_{k}}\equiv[V_{k}V_{\bar{k}}^{\dagger}]_{i_{k}j_{k}}\,.$ (60)

It is possible to generalize the method of reaching a system of differential
equations showed in the previous section to an arbitrary number $K$ pairs of
Wilson lines. The steps are outlined in App. C. This procedure leads to the
differential equation,

$\displaystyle\frac{2N_{c}}{g^{2}}\frac{{\rm d}}{{\rm
d}t}\left\langle\prod_{n=1}^{K}W^{n}_{i_{n}j_{n}}\right\rangle$
$\displaystyle=n(t)\left[\sum_{k=1}^{K-1}\sum_{l>k}^{K}(\gamma_{kl}+\gamma_{\bar{k}\bar{l}}-\gamma_{k\bar{l}}-\gamma_{\bar{k}l})-\sum_{k=1}^{K}\gamma_{k\bar{k}}-K(N_{c}^{2}-1)\gamma_{0}\right]\left\langle\prod_{n=1}^{K}W^{n}_{i_{n}j_{n}}\right\rangle$
$\displaystyle+n(t)\sum_{k=1}^{K}\left[\gamma_{k\bar{k}}\left\langle\operatorname{tr}(W^{k})\delta_{i_{k}j_{k}}\left(\prod_{n\neq
k}^{K}W^{n}_{i_{n}j_{n}}\right)\right\rangle\right]$
$\displaystyle+n(t)\sum_{k=1}^{K-1}\sum_{l>k}^{K}\Bigg{\langle}\left(\gamma_{k\bar{l}}\delta_{i_{k}j_{l}}[W^{l}W^{k}]_{i_{l}j_{k}}+\gamma_{\bar{k}l}\delta_{i_{l}j_{k}}[W^{k}W^{l}]_{i_{k}j_{l}}-\gamma_{kl}W^{k}_{i_{l}j_{k}}W^{l}_{i_{k}j_{l}}-\gamma_{\bar{k}\bar{l}}W^{k}_{i_{k}j_{l}}W^{l}_{i_{l}j_{k}}\right)$
$\displaystyle\quad\times\prod_{n\neq k,n\neq
l}^{K}W^{n}_{i_{n}j_{n}}\Bigg{\rangle}\,.$ (61)

One can see that the term on the first line has the same index structure as
the original, while the subsequent lines contain mixing terms. Notice that, in
the mixing terms, only at most two $W$’s change place. The rest stay the same
as before.

The above equation is a step in the right direction. It makes it possible to
quite easily project out all the different differential equations by
contraction with the product of $K$ Kronecker deltas. For example starting
with (4.2) and projecting out with
$\delta_{j_{1}i_{1}}\delta_{j_{2}i_{2}}\dots\delta_{j_{K}i_{K}}$ turns it into
a differential equation for $\frac{{\rm d}}{{\rm d}t}\langle
tr[W^{1}]\operatorname{tr}[W^{2}]\dots\operatorname{tr}[W^{K}]\rangle$, while
$\delta_{j_{1}i_{2}}\delta_{j_{2}i_{3}}\dots\delta_{j_{K}i_{1}}$ leads to
$\frac{{\rm d}}{{\rm d}t}\langle\operatorname{tr}[W^{1}W^{2}\dots
W^{K}]\rangle$. We will denote these two possibilities by $C_{12\dots
K}\equiv\langle\operatorname{tr}[W^{1}]\operatorname{tr}[W^{2}]\dots\operatorname{tr}[W^{K}]\rangle$
and $C_{23\dots K1}\equiv\langle\operatorname{tr}[W^{1}W^{2}\dots
W^{K}]\rangle$. The general version of this is $C_{m_{1}m_{2}\dots m_{K}}$,
where $m_{1}m_{2}\dots m_{K}$ is one of the $K!$ permutations of the numbers
between $1$ and $K$. The idea behind this notation is that $W^{1}$ is
connected to $W^{m_{1}}$, $W^{2}$ is connected to $W^{m_{2}}$ etc.333One final
example to clarify the notation can, for instance, be the correlator
$C_{213\ldots
K}\equiv\langle\operatorname{tr}[W^{1}W^{2}]\operatorname{tr}[W^{3}]\ldots\operatorname{tr}[W^{K}]\rangle$.

Although it is possible to use (4.2) to project out all the necessary
differential equations, there are actually $K!$ such projections, which
quickly becomes a huge number. It would be much preferable to write this
system in matrix form, like in Eq. (47). Making use of the notation we
described above we want to write the system of differential equations for $K$
pairs of Wilson lines as

$\frac{{\rm d}}{{\rm d}t}C_{m_{1}m_{2}\dots
m_{K}}=-\frac{1}{2}n(t)\sum_{p_{1}p_{2}\dots p_{K}}\mathbb{M}_{m_{1}m_{2}\dots
m_{K}}^{p_{1}p_{2}\dots p_{K}}C_{p_{1}p_{2}\dots p_{K}}\,,$ (62)

where $p_{1}p_{2}\dots p_{K}$ also is one of the $K!$ permutations of $12\dots
K$.

Starting from (4.2), one can deduce the general form of the matrix elements
$\mathbb{M}_{m_{1}m_{2}\dots m_{K}}^{p_{1}p_{2}\dots p_{K}}$. For details on
how this is done, we refer to App. C. Fortunately, most of the matrix elements
are zero, and those that are not have quite simple expressions. The $K!$
diagonal entries are

$\displaystyle\mathbb{M}_{m_{1}m_{2}\dots m_{K}}^{m_{1}m_{2}\dots m_{K}}$
$\displaystyle=N_{c}\sum_{k=1}^{K}\sigma_{\bar{k}m_{k}}+\frac{1}{N_{c}}\underbrace{\sum_{k=1}^{K}\sum_{l>k}^{K}(\sigma_{kl}+\sigma_{\bar{k}\bar{l}}-\sigma_{k\bar{l}}-\sigma_{\bar{k}l})-\frac{1}{N_{c}}\sum_{k=1}^{K}\sigma_{k\bar{k}}}_{A_{K}}$
(63)

Note here that only the first sum depends on the exact permutation we use. The
two latter sums are independent of this, and are common to all the diagonal
terms, so we call it $A_{K}$. The only other non-zero matrix elements
$\mathbb{M}_{m_{1}m_{2}\dots m_{K}}^{p_{1}p_{2}\dots p_{K}}$ are those where
$p_{1}p_{2}\dots p_{K}$ is just $m_{1}m_{2}\dots m_{K}$, but with exactly two
entries swapped places. If our original sequence is $m_{1}m_{2}\dots
m_{i}\dots m_{j}\dots m_{K}$, and its entries in positions $i$ and $j$ have
changed places it becomes $m_{1}m_{2}\dots m_{j}\dots m_{i}\dots m_{K}$. Then
we get $K!\frac{K(K-1)}{2}$ entries of the form

$\mathbb{M}_{m_{1}m_{2}\dots m_{i}\dots m_{j}\dots m_{K}}^{m_{1}m_{2}\dots
m_{j}\dots m_{i}\dots
m_{K}}=\sigma_{\bar{i}m_{j}}+\sigma_{m_{i}\bar{j}}-\sigma_{m_{i}m_{j}}-\sigma_{\bar{i}\bar{j}}\,.$
(64)

Finally, we have

$\mathbb{M}_{m_{1}m_{2}\dots m_{K}}^{p_{1}p_{2}\dots p_{K}}=0\,,$ (65)

for $p_{1}p_{2}\dots p_{K}$ being any other permutation of $m_{1}m_{2}\dots
m_{K}$. This means that out of the $K!^{2}$ matrix elements, only
$\frac{1}{2}K!(K^{2}-K+2)$ are non-zero. These are given by the relatively
simple formulas (63) and (64). Putting it all together this becomes

$\displaystyle\frac{{\rm d}}{{\rm d}t}C_{m_{1}m_{2}\dots m_{i}\dots m_{j}\dots
m_{K}}=$
$\displaystyle-\frac{1}{2}n(t)\big{(}N_{c}\sum_{k=1}^{K}\sigma_{\bar{k}m_{k}}+\frac{1}{N_{c}}A_{K}\big{)}C_{m_{1}m_{2}\dots
m_{i}\dots m_{j}\dots m_{K}}$
$\displaystyle-\frac{1}{2}n(t)\sum_{i=1}^{K-1}\sum_{j>i}^{K}(\sigma_{\bar{i}m_{j}}+\sigma_{m_{i}\bar{j}}-\sigma_{m_{i}m_{j}}-\sigma_{\bar{i}\bar{j}})C_{m_{1}m_{2}\dots
m_{j}\dots m_{i}\dots m_{K}}\,.$ (66)

Of course, for all differential equations you need to specify some initial
conditions. It is clear from the definition of the Wilson line (4) that
$V_{ij}(t_{0},t_{0})=\delta_{ij}$. The trace of this is
$\operatorname{tr}V(t_{0},t_{0})$=$N_{c}$. This means that the initial
condition of a Wilson line correlator is $N_{c}$ to the power of traces it
contains. A few illustrative examples of this are $C_{12\dots
K}=\langle\operatorname{tr}[W^{1}]\operatorname{tr}[W^{2}]\dots\operatorname{tr}[W^{K}]\rangle\sim
N_{c}^{K}$, $C_{213\dots
K}=\langle\operatorname{tr}[W^{1}W^{2}]\operatorname{tr}[W^{3}]\dots\operatorname{tr}[W^{K}]\rangle\sim
N_{c}^{K-1}$ and $C_{23\dots K1}=\langle\operatorname{tr}[W^{1}W^{2}\dots
W^{K}]\rangle\sim N_{c}^{1}$.

To better understand what this system of differential equations looks like, it
is useful to view it in matrix form. Generally, there will be several
correlators that go as the same power of $N_{c}$. It is useful to gather these
in vectors $\bm{C^{M}}$, where the superscript $M$ is meant to indicate that
this scales as $N_{c}^{M}$. Then, Eq. (4.2) can be represented as

$\frac{{\rm d}}{{\rm d}t}\begin{bmatrix}C^{K}\\\ \bm{C^{K-1}}\\\
\bm{C^{K-2}}\\\ \vdots\\\ \bm{C^{2}}\\\
\bm{C^{1}}\end{bmatrix}\sim\left(\textrm{diag}\left(N_{c}\sigma+\frac{1}{N_{c}}\sigma,\dots,N_{c}\sigma+\frac{1}{N_{c}}\sigma\right)+\begin{bmatrix}0&\bm{\sigma}&\bm{0}&\dots&\dots&\bm{0}\\\
\bm{\sigma}&\bm{0}&\bm{\sigma}&\bm{0}&\dots&\bm{0}\\\
\bm{0}&\bm{\sigma}&\bm{0}&\bm{\sigma}&\bm{0}&\bm{0}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\bm{0}&\dots&\bm{0}&\bm{\sigma}&\bm{0}&\bm{\sigma}\\\
\bm{0}&\dots&\dots&\bm{0}&\bm{\sigma}&\bm{0}\\\
\end{bmatrix}\right)\begin{bmatrix}C^{K}\\\ \bm{C^{K-1}}\\\ \bm{C^{K-2}}\\\
\vdots\\\ \bm{C^{2}}\\\ \bm{C^{1}}\end{bmatrix}\,.$ (67)

The first matrix contains the diagonal elements, written in detail in (63).
The second matrix represents the non-diagonal elements, and $\bm{\sigma}$ is a
block containing non-zero elements, which we get from Eq. (64).

###  4.3 Wilson line correlators in the large-$N_{c}$ limit

The system of differential equations given by Eq. (4.2) is to our knowledge
not possible to solve analytically in the case where $\sigma$ is a function of
time, so we have to turn to numerical techniques. However, in the
large-$N_{c}$ limit, the system simplifies in a way that makes it possible to
solve it exactly. This can be seen from the matrix representation in Eq. (67).
Since $\sigma\sim N_{c}^{-1}$, the diagonal matrix elements go as $\sim
N_{c}^{0}+N_{c}^{-2}$, while the non-diagonal ones go as $\sim N_{c}^{-1}$.
Multiplying in the vector on the end and representing every term by its
$N_{c}$ scaling this becomes

$\frac{{\rm d}}{{\rm d}t}\begin{bmatrix}C^{K}\\\ \bm{C^{K-1}}\\\
\bm{C^{K-2}}\\\ \vdots\\\ \bm{C^{2}}\\\
\bm{C^{1}}\end{bmatrix}\sim\begin{bmatrix}(N_{c}^{0}+N_{c}^{-2})C^{K}\\\
(\bm{N_{c}^{0}}+\bm{N_{c}^{-2}})\bm{C^{K-1}}\\\
(\bm{N_{c}^{0}}+\bm{N_{c}^{-2}})\bm{C^{K-2}}\\\ \vdots\\\
(\bm{N_{c}^{0}}+\bm{N_{c}^{-2}})\bm{C^{2}}\\\
(\bm{N_{c}^{0}}+\bm{N_{c}^{-2}})\bm{C^{1}}\end{bmatrix}+\begin{bmatrix}\bm{N_{c}^{-1}}\bm{C^{K-1}}\\\
\bm{N_{c}^{-1}}(\bm{C^{K}}+\bm{C^{K-2}})\\\
\bm{N_{c}^{-1}}(\bm{C^{K-1}}+\bm{C^{K-3}})\\\ \vdots\\\
\bm{N_{c}^{-1}}(\bm{C^{3}}+\bm{C^{1}})\\\
\bm{N_{c}^{-1}}\bm{C^{2}}\end{bmatrix}\,.$ (68)

Taking the large-$N_{c}$ limit is equivalent to keeping only the leading order
of $N_{c}$ in each row, and dropping terms going as $N_{c}^{-2}$ compared to
the leading term. Translating this back to the form in Eq. (67), this becomes

$\frac{{\rm d}}{{\rm d}t}\begin{bmatrix}C^{K}\\\ \bm{C^{K-1}}\\\
\bm{C^{K-2}}\\\ \vdots\\\ \bm{C^{2}}\\\
\bm{C^{1}}\end{bmatrix}\sim\left(\textrm{diag}\left(N_{c}\sigma,\dots,N_{c}\sigma\right)+\begin{bmatrix}0&\bm{0}&\bm{0}&\dots&\dots&\bm{0}\\\
\bm{\sigma}&\bm{0}&\bm{0}&\bm{0}&\dots&\bm{0}\\\
\bm{0}&\bm{\sigma}&\bm{0}&\bm{0}&\bm{0}&\bm{0}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\bm{0}&\dots&\bm{0}&\bm{\sigma}&\bm{0}&\bm{0}\\\
\bm{0}&\dots&\dots&\bm{0}&\bm{\sigma}&\bm{0}\\\
\end{bmatrix}\right)\begin{bmatrix}C^{K}\\\ \bm{C^{K-1}}\\\ \bm{C^{K-2}}\\\
\vdots\\\ \bm{C^{2}}\\\ \bm{C^{1}}\end{bmatrix}\,.$ (69)

Hence, in the large-$N_{c}$ limit, all of the terms above the diagonal go to
zero, and the system simplifies drastically. To get some more intuition into
why this is true physically it is useful to imagine being in some color
configuration that scales as $\sim N_{c}^{M}$. At any point it is possible to
exchange one gluon, after which the possible resulting color configurations of
the system will change its $N_{c}$ power by exactly one, and go as $\sim
N_{c}^{M+1}$ or $\sim N_{c}^{M-1}$. The gluon exchange comes with a factor
$\sigma\sim N_{c}^{-1}$, so in total the overall $N_{c}$ power of going to
these systems are $N_{c}^{M}$ and $N_{c}^{M-2}$. In the large-$N_{c}$
approximation the latter possibility is discarded, which is equivalent to
dropping all the terms above the diagonal in the matrix (69).

It is clear from this discussion that the system of differential equations
(4.2) simplifies, in the large-$N_{c}$ limit, to

$\displaystyle\frac{{\rm d}}{{\rm d}t}C_{m_{1}m_{2}\dots m_{i}\dots m_{j}\dots
m_{K}}^{M}\simeq$
$\displaystyle-\frac{1}{2}n(t)N_{c}\sum_{k=1}^{K}\sigma_{\bar{k}m_{k}}C_{m_{1}m_{2}\dots
m_{i}\dots m_{j}\dots m_{K}}^{M}$
$\displaystyle-\frac{1}{2}n(t)\sum_{i=1}^{K-1}\sum_{j>i}^{K}(\sigma_{\bar{i}m_{j}}+\sigma_{m_{i}\bar{j}}-\sigma_{m_{i}m_{j}}-\sigma_{\bar{i}\bar{j}})C^{M+1}_{m_{1}m_{2}\dots
m_{j}\dots m_{i}\dots m_{K}}\,.$ (70)

Here we have included superscripts to show the $N_{c}$-scaling. In the second
line, we have indicated that only the correlators scaling as $N_{c}^{M+1}$
should be included in the sum. This means that in the large-$N_{c}$ limit the
correlators with $M$ traces only depend on the correlators with $M+1$ traces.
Similarly, the correlators with $M+1$ traces depend on the correlators with
$M+2$ traces and so on. This continues all the way up to the correlators with
$K-1$ traces, which depend on the correlators with $K$ traces. Using (4.3) the
differential equation for the correlator scaling as $N_{c}^{K}$ is

$\displaystyle\frac{{\rm d}}{{\rm d}t}C_{12\dots K}$
$\displaystyle\simeq-\frac{1}{2}n(t)N_{c}\sum_{k=1}^{K}\sigma_{\bar{k}k}C_{12\dots
K}\,.$ (71)

Since this is exactly solvable,

$C_{12\dots K}(t)=N_{c}^{K}\,{\rm e}^{-\frac{1}{2}N_{c}\int^{t}_{t_{2}}{\rm
d}s\,n(s)\sum_{k=1}^{K}\sigma_{\bar{k}k}(s)}\,,$ (72)

this provides a “bootstrap” for the whole system of equations. The above
argument shows that in principle all the correlators can be solved exactly in
the large-$N_{c}$ limit.

As a side note, we can also understand the large-$N_{c}$ approximation as a
simplification of the operation of performing medium averages on multiple
traced correlators. Given that a dipole in the large-$N_{c}$ is given by

$\mathcal{S}_{1\bar{1}}(t,t_{2})\equiv\frac{1}{N_{c}}\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}]\rangle={\rm
e}^{-\frac{1}{2}N_{c}\int_{t_{2}}^{t}{\rm d}s\,n(s)\sigma_{1\bar{1}}}\,,$ (73)

the answer for $C_{1\ldots K}(t)$ is just given by the product of $K$ dipoles,
i.e. $C_{1\dots K}\simeq
N_{c}^{K}\mathcal{S}_{1\bar{1}}\dots\mathcal{S}_{K\bar{K}}$. On the level of
the full correlator, this corresponds to the simplification

$\displaystyle\langle\operatorname{tr}[W^{1}]\operatorname{tr}[W^{2}]\dots\operatorname{tr}[W^{K}]\rangle\approx\langle\operatorname{tr}[W^{1}]\rangle\langle\operatorname{tr}[W^{2}]\rangle\dots\langle\operatorname{tr}[W^{K}]\rangle\,.$
(74)

This argument can also be extended to any of the other correlators discussed
above, e.g. $\langle\operatorname{tr}[W^{1}]\operatorname{tr}[W^{2}\ldots
W^{K}]\rangle\approx\langle\operatorname{tr}[W^{1}]\rangle\langle\operatorname{tr}[W^{2}\ldots
W^{K}]\rangle$.

The simplified differential equation, Eq. (4.3), can also be solved directly
to get the recursive formula

$\displaystyle C_{m_{1}m_{2}\dots m_{i}\dots m_{j}\dots
m_{K}}^{M}=N_{c}^{M}\,{\rm e}^{-\frac{1}{2}N_{c}\int_{t_{2}}^{t}{\rm
d}s\,n(s)\sum_{k=1}^{K}\sigma_{\bar{k}m_{k}}}$
$\displaystyle-\frac{1}{2}\int_{t_{2}}^{t}{\rm
d}s\,n(s)\sum_{i=1}^{K-1}\sum_{j>i}^{K}\left((\sigma_{\bar{i}m_{j}}+\sigma_{m_{i}\bar{j}}-\sigma_{m_{i}m_{j}}-\sigma_{\bar{i}\bar{j}})C^{M+1}_{m_{1}m_{2}\dots
m_{j}\dots m_{i}\dots m_{K}}\right)$ $\displaystyle\times{\rm
e}^{-\frac{1}{2}N_{c}\int_{s}^{t}{\rm
d}s^{\prime}\,n(s^{\prime})\sum_{k=1}^{K}\sigma_{\bar{k}m_{k}}}\,.$ (75)

This can also be written in terms of dipoles, namely

$\displaystyle C_{m_{1}m_{2}\dots m_{i}\dots m_{j}\dots
m_{K}}^{M}=N_{c}^{M}\prod_{k=1}^{K}\mathcal{S}_{m_{k}\bar{k}}(t,t_{2})$
$\displaystyle-\frac{1}{2}\int_{t_{2}}^{t}{\rm
d}s\,n(s)\sum_{i=1}^{K-1}\sum_{j>i}^{K}\left((\sigma_{\bar{i}m_{j}}+\sigma_{m_{i}\bar{j}}-\sigma_{m_{i}m_{j}}-\sigma_{\bar{i}\bar{j}})C^{M+1}_{m_{1}m_{2}\dots
m_{j}\dots m_{i}\dots
m_{K}}\right)\prod_{k=1}^{K}\mathcal{S}_{m_{k}\bar{k}}(t,s)\,.$ (76)

From this equation it is clear that all of the Wilson line correlators can be
written in terms of dipoles in the large-$N_{c}$ limit. That is because Eq.
(4.3) is a recursive relation (“bootstrap”) that stops when you reach the term
with $K$ traces, which is given in terms of dipoles in (74). Since the only
Wilson line structure that appears in both (74) and (4.3) is dipoles, it means
all the correlator can be written in terms of dipoles. In Ref.
Dominguez:2012ad  it was pointed out that that all higher-order correlators
can be reduced to dipoles and quadrupoles at large-$N_{c}$. The result in this
section directly confirms this, and show that really only dipoles are needed.

We could, in principle, also devise a scheme to compute sub-leading color
corrections, that scale like $N_{c}^{-2}$ relative to the leading terms,
following the steps in Eqs. (53) and (54), and below. We have nevertheless not
pursued this program further in this work.

##  5 Conclusion and outlook

In this paper we have developed a general method for calculating correlators
involving an arbitrary number of Wilson lines in the fundamental
representation. This culminated in the system of differential equations in Eq.
(4.2). This system can be solved numerically. We showed that in the
large-$N_{c}$ limit the resulting simplified system of differential equations,
Eq. (4.3), can be solved exactly. We also provided a general way to compute
color sub-leading corrections, suppressed by $N_{c}^{-2}$ relative to the
leading terms. This was done in detail for the four-point correlator, in Eqs.
(53) and (54), but can easily be extended to any higher-order correlator. All
the results can then be written in terms of dipoles and their convolutions.

This technique was applied on three different cases of $1\to 2$ parton
splittings in the medium, which were shown to involve correlators containing
up to eight (fundamental) Wilson lines. We used our method to calculate these
both at finite and large $N_{c}$. Comparisons of the results are shown in Fig.
4. From these plots it is clear that in this exact case the large-$N_{c}$
approximations works quite well for small $\theta$, but the differences become
bigger as $\theta$ grows. In certain areas of the phase space the error in
using the large-$N_{c}$ limit might be as high as $16\%$. This is expected
given that the corrections we find generically scale as $N_{c}^{-2}$.

Since our method deals with a generic set of correlated Wilson lines,
representing particles moving on eikonal trajectories through a background
field, it could easily be extended to many other physical situations. For
future work it would be interesting to apply our results in initial state
physics, where similar correlators of Wilson lines also appear, and for soft
contributions to event or jet observables in electron-positron or proton-
proton collisions. Finally, we plan on extending the formulation to account
for non-eikonal corrections to the particle trajectories.

###### Acknowledgements.

We would like to thank A. Takacs for useful discussions. This work is
supported by a Starting Grant from Trond Mohn Foundation (BFS2018REK01) and
the University of Bergen.

##  Appendix A Calculation of spectrums

Here we will show the calculations leading up to the for the emission spectra
$\frac{{\rm d}I}{{\rm d}z{\rm d}\theta}$. The Feynman rules from
Mehtar_Tani_2018  have been used to calculate the matrix elements.

###  A.1 Pair production

We start with the process of a photon producing a quark-antiquark pair. This
process has been calculated in Dom_nguez_2020  but we will restate some of the
results. The amplitude is

$\displaystyle{\cal M}_{s_{1},s_{2}}^{ij}$
$\displaystyle=\int_{{\bm{p}}_{0},{\bm{p}}_{1}^{\prime},{\bm{p}}_{2}^{\prime}}\int_{t_{0}}^{L}{\rm
d}t_{1}\,(2\pi)^{2}\delta({\bm{p}}_{0}-{\bm{p}}_{1}^{\prime}-{\bm{p}}_{2}^{\prime})[({\bm{p}}_{1}|{\cal
G}_{F}(L,t_{1})|{\bm{p}}_{1}^{\prime})({\bm{p}}_{2}^{\prime}|\bar{{\cal
G}}_{F}(L,t_{1})|{\bm{p}}_{2})]^{ij}$ $\displaystyle\times
A_{\lambda,s_{1},s_{2}}({\bm{p}}_{2}^{\prime}-z{\bm{p}}_{0},z)\frac{1}{2E}{\rm
e}^{-i\frac{{\bm{p}}_{0}^{2}}{2E}(t_{1}-t_{0})}{\cal
M}_{0\lambda}({\bm{p}}_{0})\,,$ (77)

where the photon-quark vertex is given by

$A_{\lambda,s_{1},s_{2}}({\bm{q}},z)=\frac{2ie}{\sqrt{z(1-z)}}\delta_{-s_{2}s_{1}}(z\delta_{\gamma
s_{1}}-(1-z)\delta_{\gamma s_{2}}){\bm{q}}\cdot\epsilon_{\gamma}\,.$ (78)

The initial hard process is represented by the amplitude ${\cal M}_{0}$. After
using the eikonal approximation (13) this becomes (up to some phase that
cancels when we take the square)

$\displaystyle{\cal M}_{s_{1},s_{2}}^{ij}$
$\displaystyle=\frac{1}{2E}\int_{t_{0}}^{L}{\rm
d}t_{1}\,e^{i\frac{1}{2z(1-z)E}((1-z){\bm{p}}_{2}-z{\bm{p}}_{1})^{2}t_{1}}[V_{1}(L,t_{1})V_{2}^{\dagger}(t_{1},L)]^{ij}$
$\displaystyle\times
A_{\lambda,s_{1},s_{2}}((1-z){\bm{p}}_{2}-z{\bm{p}}_{1},z){\cal
M}_{0\lambda}({\bm{p}}_{1}+{\bm{p}}_{2})\,.$ (79)

We have used the more compact notation to write $V_{F}({\bm{r}}_{1})\equiv
V_{1}$, $V_{F}({\bm{r}}_{2})\equiv V_{2}$. The goal is to calculate

$\frac{{\rm d}I}{{\rm d}z\,{\rm
d}\theta}=\frac{z(1-z)E^{2}\theta}{8\pi^{2}}\frac{\left\langle|{\cal
M}|^{2}\right\rangle}{\left\langle|{\cal M}_{0}|^{2}\right\rangle}\,.$ (80)

The Wilson lines can be split using $V(L,t_{1})=V(L,t_{2})V(t_{2},t_{1})$.
Then we only need to deal with the two time intervals $(L,t_{2})$ and
$(t_{2},t_{1})$. After squaring the amplitude, averaging over initial
polarization, summing the final spins, flavor and colors and taking the medium
average this becomes (18) with (22) and (23).

###  A.2 Quark-gluon splitting

The amplitude was calculated in Mehtar_Tani_2018  and is

$\displaystyle{\cal M}_{\lambda,s}^{ai}$
$\displaystyle=\int_{{\bm{p}}_{0},{\bm{p}}_{0}^{\prime},{\bm{k}}^{\prime},{\bm{p}}^{\prime}}\int_{t_{0}}^{L}{\rm
d}t_{1}\,(2\pi)^{2}\delta({\bm{p}}_{0}^{\prime}-{\bm{k}}^{\prime}-{\bm{p}}^{\prime})({\bm{k}}|{\cal
G}_{A}^{ab}(L,t_{1})|{\bm{k}}^{\prime})$ $\displaystyle\times[({\bm{p}}|{\cal
G}_{F}(L,t_{1})|{\bm{p}}^{\prime})A^{b}_{\lambda,s,s^{\prime}}({\bm{k}}^{\prime}-z{\bm{p}}_{0}^{\prime},z)\frac{1}{2E}({\bm{p}}_{0}^{\prime}|{\cal
G}_{F}(t_{1},t_{0})|{\bm{p}}_{0})]^{ij}{\cal
M}_{0s^{\prime}}^{j}({\bm{p}}_{0})\,,$ (81)

where the quark-gluon vertex is

$\displaystyle
A^{aij}_{\lambda,s,s^{\prime}}(\bm{q},z)=-\frac{2ig\mathbf{t}^{a}_{ij}}{z\sqrt{1-z}}\delta_{s^{\prime}s}\left[\delta_{\lambda
s}+(1-z)\delta_{\lambda-s}\right]\bm{q}\cdot\bm{\epsilon}_{\lambda}^{*}\,.$
(82)

Again this simplifies in the eikonal limit (13)

${\cal M}_{\lambda,s}^{ai}=\frac{1}{2E}\int_{t_{0}}^{L}{\rm
d}t_{1}\,e^{i\frac{1}{2z(1-z)E}((1-z){\bm{k}}-z{\bm{p}})^{2}t_{1}}U_{2}^{ab}(L,t_{1})[V_{1}(L,t_{1})A^{b}_{\lambda,s,s^{\prime}}V_{0}(t_{1},t_{0})]^{ij}{\cal
M}_{0s^{\prime}}^{j}({\bm{k}}+{\bm{p}})\,.$ (83)

We have denoted the adjoint Wilson line as $V_{A}({\bm{r}}_{2})\equiv U_{2}$.
Squaring the amplitude, summing/averaging over spins and colors and taking the
medium average gives

$\displaystyle\frac{\mathrm{d}I}{\mathrm{d}z\,\mathrm{d}\theta}$
$\displaystyle=\frac{\alpha_{s}}{\pi}\frac{P_{gq}(z)}{\theta}\frac{2}{N_{c}^{2}-1}2\operatorname{Re}\int_{t_{0}}^{L}\frac{\mathrm{d}t_{1}}{t_{\mathrm{f}}}\int_{t_{1}}^{L}\frac{\mathrm{d}t_{2}}{t_{\mathrm{f}}}\mathrm{e}^{-i\frac{t_{2}-t_{1}}{t_{\mathrm{f}}}}$
$\displaystyle\times\langle\left[U^{\dagger}(t_{2},L)U(L,t_{1})\right]^{\bar{b}b}\operatorname{tr}\left[V_{0}^{\dagger}(0,t_{2})t^{\bar{b}}V_{\bar{1}}^{\dagger}(t_{2},L)V_{1}(L,t_{1})t^{b}V_{0}(t_{1},0)\right]\rangle\,,$
(84)

where the relevant Altarelli-Parisi splitting function is

$P_{gq}(z)=C_{F}\frac{1+(1-z)^{2}}{z}\,.$ (85)

To continue we transform the adjoint Wilson lines into fundamental ones using
the identity (5). The resulting expression will contain many group generators
$t^{a}$, and can be simplified by using the Fierz identity (44). Finally,
completely in the fundamental representation the Wilson line structure becomes

$\displaystyle\left\langle\left[U^{\dagger}(t_{2},L)U(L,t_{1})\right]^{\bar{b}b}\operatorname{tr}\left[V_{0}^{\dagger}(0,t_{2})t^{\bar{b}}V_{\bar{1}}^{\dagger}(t_{2},L)V_{1}(L,t_{1})t^{b}V_{0}(t_{1},0)\right]\right\rangle$
$\displaystyle=\frac{1}{2}\big{\langle}\left([V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}V_{1}]_{kj}[V_{\bar{2}}^{\dagger}V_{2}]_{il}-\frac{1}{N_{c}}[V_{\bar{1}}^{\dagger}V_{1}]_{ij}\delta_{kl}\right)_{(L,t_{2})}$
$\displaystyle\times\left([V_{1}V_{2}^{\dagger}]_{jk}[V_{2}V_{0}^{\dagger}]_{li}-\frac{1}{N_{c}}[V_{1}V_{0}^{\dagger}]_{ji}\delta_{lk}\right)_{(t_{2},t_{1})}\big{\rangle}\,.$
(86)

Conservation of color then makes it possible to connect $i,l$ and $j,k$ so
when we include the proper normalization factor the whole expression turns
into (18) with (31) and (32).

###  A.3 Gluon-gluon splitting

The calculation of the emission spectrum for gluon-gluon splittings was done
in Blaizot_2013 . For completeness we will also include the main results here.
The matrix element of the process is

$\displaystyle{\cal
M}_{\lambda_{1},\lambda_{2}}^{a_{1}a_{2}}=\int_{{\bm{k}}_{0},{\bm{k}}_{0}^{\prime},{\bm{k}}_{1}^{\prime},{\bm{k}}_{2}^{\prime}}\int_{t_{0}}^{L}{\rm
d}t_{1}\,(2\pi)^{2}\delta({\bm{k}}_{0}^{\prime}-{\bm{k}}_{1}^{\prime}-{\bm{k}}_{2}^{\prime})$
$\displaystyle\times({\bm{k}}_{1}|{\cal
G}_{A}^{a_{1}b_{1}}(L,t_{1})|{\bm{k}}_{1}^{\prime})({\bm{k}}_{2}|{\cal
G}_{A}^{a_{2}b_{2}}(L,t_{1})|{\bm{k}}_{2}^{\prime})A^{b_{0}b_{1}b_{2}}_{\lambda_{0},\lambda_{1},\lambda_{2}}({\bm{k}}_{2}^{\prime}-z{\bm{k}}_{0}^{\prime},z)\frac{1}{2E}({\bm{k}}_{0}^{\prime}|{\cal
G}_{A}^{b_{0}c}(t_{1},t_{0})|{\bm{k}}_{0})$ $\displaystyle\times{\cal
M}_{0\lambda_{0}}^{c}({\bm{k}}_{0})\,,$ (87)

where the gluon-gluon vertex is

$\displaystyle
A^{b_{0}b_{1}b_{2}}_{\lambda_{0},\lambda_{1},\lambda_{2}}(\bm{q},z)=-2ig(T^{b_{0}})^{b_{1}b_{2}}\left[\frac{1}{z}({\bm{q}}\cdot\bm{\epsilon}_{\lambda_{2}}^{*})\delta_{\lambda_{0}\lambda_{1}}+\frac{1}{1-z}({\bm{q}}\cdot\bm{\epsilon}_{\lambda_{1}}^{*})\delta_{\lambda_{0}\lambda_{2}}-({\bm{q}}\cdot\bm{\epsilon}_{\lambda_{0}})\delta_{\lambda_{1}\lambda_{2}}\right]\,.$
(88)

In the eikonal approximation (13) the amplitude is

$\displaystyle{\cal M}_{\lambda_{1},\lambda_{2}}^{a_{1}a_{2}}$
$\displaystyle=\frac{1}{2E}\int_{t_{0}}^{L}{\rm
d}t_{1}\,e^{i\frac{1}{2z(1-z)E}((1-z){\bm{k}}_{2}-z{\bm{k}}_{1})^{2}t_{1}}U_{1}^{a_{1}b_{1}}(L,t_{1})U_{2}^{a_{2}b_{2}}(L,t_{1})A^{b_{0}b_{1}b_{2}}_{\lambda_{0},\lambda_{1},\lambda_{2}}U_{0}^{b_{0}c}(t_{1},t_{0})$
$\displaystyle\times{\cal M}_{0s^{\prime}}^{c}({\bm{k}}_{1}+{\bm{k}}_{2})\,.$
(89)

After taking the square of the amplitude, summing/averaging over spins and
colors and taking the medium average this becomes

$\displaystyle\frac{\mathrm{d}I}{\mathrm{d}z\,\mathrm{d}\theta}$
$\displaystyle=\frac{\alpha_{s}}{\pi}\frac{P_{gg}(z)}{\theta}\frac{2}{N_{c}(N_{c}^{2}-1)}2\operatorname{Re}\int_{t_{0}}^{L}\frac{\mathrm{d}t_{1}}{t_{\mathrm{f}}}\int_{t_{1}}^{L}\frac{\mathrm{d}t_{2}}{t_{\mathrm{f}}}\mathrm{e}^{-i\frac{t_{2}-t_{1}}{t_{\mathrm{f}}}}$
$\displaystyle\times
f^{b_{0}b_{1}b_{2}}f^{\bar{b}_{0}\bar{b}_{1}\bar{b}_{2}}\langle[U_{1}^{a_{1}d_{1}}U_{2}^{a_{2}d_{2}}U_{\bar{1}}^{{\dagger}\bar{b}_{1}a_{1}}U_{\bar{2}}^{{\dagger}\bar{b}_{2}a_{2}}]_{(L,t_{2})}[U_{1}^{d_{1}b_{1}}U_{2}^{d_{2}b_{2}}U_{0}^{{\dagger}b_{0}\bar{b}_{0}}]_{(t_{2},t_{1})}\rangle\,,$
(90)

where the relevant Altarelli-Parisi splitting function is

$P_{gg}(z)=N_{c}\left[z(1-z)+\frac{1-z}{z}+\frac{z}{1-z}\right]\,.$ (91)

Conservation of color lets us decouple the Wilson lines in the two time
intervals $(L,t_{2})$ and $(t_{2},t_{1})$

$f^{b_{0}b_{1}b_{2}}U_{1}^{d_{1}b_{1}}U_{2}^{d_{2}b_{2}}U_{0}^{{\dagger}b_{0}\bar{b}_{0}}=\frac{1}{N_{c}(N_{c}^{2}-1)}f^{d_{1}d_{2}\bar{b}_{0}}f^{d_{1}^{\prime}d_{2}^{\prime}\bar{b}_{0}^{\prime}}f^{b_{0}b_{1}b_{2}}U_{1}^{d_{1}^{\prime}b_{1}}U_{2}^{d_{2}^{\prime}b_{2}}U_{0}^{{\dagger}b_{0}\bar{b}_{0}^{\prime}}\,.$
(92)

The part in the time interval $(t_{2},t_{1})$ can be calculated explicitly
because of its simple color structure

$\displaystyle\frac{1}{N_{c}(N_{c}^{2}-1)}f^{d_{1}^{\prime}d_{2}^{\prime}\bar{b}_{0}^{\prime}}f^{b_{0}b_{1}b_{2}}U_{1}^{d_{1}^{\prime}b_{1}}U_{2}^{d_{2}^{\prime}b_{2}}U_{0}^{{\dagger}b_{0}\bar{b}_{0}^{\prime}}={\rm
e}^{-\frac{N_{c}}{2}\int_{t_{1}}^{t_{2}}{\rm
d}t\,n(t)[\sigma_{01}+\sigma_{02}+\sigma_{12}]}\,.$ (93)

What remains are the Wilson lines in time interval $(L,t_{2})$

$f^{d_{1}d_{2}\bar{b}_{0}}f^{\bar{b}_{0}\bar{b}_{1}\bar{b}_{2}}\langle[U_{1}^{a_{1}d_{1}}U_{2}^{a_{2}d_{2}}U_{\bar{1}}^{{\dagger}\bar{b}_{1}a_{1}}U_{\bar{2}}^{{\dagger}\bar{b}_{2}a_{2}}]\rangle_{(L,t_{2})}\,.$
(94)

However, these are not that easy to calculate. The procedure for calculating
Wilson line products detailed in Sec. 4 only involve fundamental Wilson lines,
so (5) is used to turn all the adjoint Wilson lines into fundamental ones.
Then one can use the definition of the structure constants
$[t^{a},t^{b}]=if^{abc}t^{c}$ and the identity (44) to get rid of all the
group generators. This was done in Blaizot_2013 , and we quote the result

$\displaystyle
f^{d_{1}d_{2}\bar{b}_{0}}f^{\bar{b}_{0}\bar{b}_{1}\bar{b}_{2}}\langle[U_{1}^{a_{1}d_{1}}U_{2}^{a_{2}d_{2}}U_{\bar{1}}^{{\dagger}\bar{b}_{1}a_{1}}U_{\bar{2}}^{{\dagger}\bar{b}_{2}a_{2}}]\rangle_{(L,t_{2})}$
$\displaystyle=\frac{1}{2}\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}]\operatorname{tr}[V_{2}V_{\bar{2}}^{\dagger}V_{\bar{1}}V_{1}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{2}^{\dagger}]-\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}V_{2}V_{\bar{2}}^{\dagger}V_{\bar{1}}V_{1}^{\dagger}V_{\bar{2}}V_{2}^{\dagger}]+\textrm{h.c.}\rangle_{(L,t_{2})}$
$\displaystyle=\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}]\operatorname{tr}[V_{2}V_{\bar{2}}^{\dagger}V_{\bar{1}}V_{1}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{2}^{\dagger}]-\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}V_{2}V_{\bar{2}}^{\dagger}V_{\bar{1}}V_{1}^{\dagger}V_{\bar{2}}V_{2}^{\dagger}]\rangle_{(L,t_{2})}\,.$
(95)

The last step is true because the medium averaged products of Wilson lines are
real. This means that in the gluon-gluon case we end up with medium averaged
products of up to eight Wilson lines. Putting it all together we get the
formula (18) with (35) and (36).

##  Appendix B Six and eight Wilson lines

###  B.1 Six lines

In Sec. 4 we developed the tools to calculate the correlators of six and eight
Wilson lines, which appeared in (31) and (35). To start we will look at the
case of six lines, which follows from (4.2) with $K=3$. The relevant
expression is $\frac{{\rm d}}{{\rm
d}t}\langle[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{1}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{2}}[V_{3}V_{\bar{3}}^{\dagger}]_{i_{3}j_{3}}\rangle$.
If this is contracted with
$\delta_{j_{1}i_{2}}\delta_{j_{2}i_{1}}\delta_{j_{3}i_{3}}$ it becomes

$\delta_{j_{1}i_{2}}\delta_{j_{2}i_{1}}\delta_{j_{3}i_{3}}\frac{{\rm d}}{{\rm
d}t}\big{\langle}[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{1}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{2}}[V_{3}V_{\bar{3}}^{\dagger}]_{i_{3}j_{3}}\big{\rangle}=\big{\langle}\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}V_{2}V_{\bar{2}}^{\dagger}]\operatorname{tr}[V_{3}V_{\bar{3}}^{\dagger}]\big{\rangle},$
(96)

which is the structure encountered in (31). To get exactly the same as in that
equation we need only change the labels
$(1,\bar{1},2,\bar{2},3,\bar{3})\to(1,2,\bar{2},\bar{1},2,\bar{2})$, which
also simplifies the system somewhat. The six different projections are
gathered into a vector

$\displaystyle\bm{C}^{\intercal}$
$\displaystyle=\big{(}C_{1\bar{2}2},C_{\bar{2}12},C_{2\bar{2}1},C_{12\bar{2}},C_{21\bar{2}},C_{\bar{2}21}\big{)}$
$\displaystyle=\Big{(}\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\operatorname{tr}[V_{2}V_{\bar{2}}^{\dagger}]\rangle,\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\operatorname{tr}[V_{2}V_{\bar{2}}^{\dagger}]\rangle,\langle\operatorname{tr}[V_{1}V_{\bar{2}}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle,$
$\displaystyle\quad\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{\bar{1}}^{\dagger}V_{2}]\rangle,\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}]\rangle,\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}V_{2}V_{\bar{2}}^{\dagger}]\rangle\Big{)}\,.$
(97)

We can write the system of differential equations as

$\frac{{\rm d}}{{\rm
d}t}\bm{C}=-\frac{1}{2}n(t)\underline{\bm{\mathbb{M}}}\bm{C}\,.$ (98)

One can get the elements of the $6\times 6$ matrix
$\underline{\bm{\mathbb{M}}}$ from (63) and (64). The 6 diagonal entries are
simply

$\displaystyle\mathbb{M}_{m_{1}m_{2}m_{3}}^{m_{1}m_{2}m_{3}}=N_{c}(\sigma_{2m_{1}}+\sigma_{\bar{1}m_{2}}+\sigma_{\bar{2}m_{3}})-\frac{1}{N_{c}}\sigma_{1\bar{1}}\,,$
(99)

where $(m_{1},m_{2},m_{3})$ now is some permutation of $(1,\bar{2},2)$. Thee
non-zero non-diagonal entries are given by

$\displaystyle\mathbb{M}_{m_{1}m_{2}m_{3}}^{m_{2}m_{1}m_{3}}$
$\displaystyle=\sigma_{2m_{2}}+\sigma_{\bar{1}m_{1}}-\sigma_{m_{2}m_{1}}-\sigma_{2\bar{1}}$
$\displaystyle\mathbb{M}_{m_{1}m_{2}m_{3}}^{m_{3}m_{2}m_{1}}$
$\displaystyle=\sigma_{2m_{3}}+\sigma_{\bar{2}m_{1}}-\sigma_{m_{3}m_{1}}-\sigma_{2\bar{2}}$
$\displaystyle\mathbb{M}_{m_{1}m_{2}m_{3}}^{m_{1}m_{3}m_{2}}$
$\displaystyle=\sigma_{\bar{1}m_{3}}+\sigma_{\bar{2}m_{2}}-\sigma_{m_{3}m_{2}}-\sigma_{\bar{1}\bar{2}}\,.$
(100)

This leads to six differential equations which can be solved numerically for
the six functions in $\bm{C}$. Interestingly this $6\times 6$ system is
reducible into two $3\times 3$ systems. The first of these systems leads to
three differential equations that actually can be solved exactly:

$\displaystyle\frac{{\rm d}}{{\rm
d}t}\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}]\rangle=$
$\displaystyle-
C_{F}n(t)\sigma_{1\bar{1}}\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}]\rangle$
$\displaystyle\frac{{\rm d}}{{\rm
d}t}\langle\operatorname{tr}[V_{1}V_{\bar{2}}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle=$
$\displaystyle-\frac{1}{2}n(t)[N_{c}(\sigma_{\bar{1}\bar{2}}+\sigma_{1\bar{2}})-\frac{1}{N_{c}}\sigma_{1\bar{1}}]\langle\operatorname{tr}[V_{1}V_{\bar{2}}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle$
$\displaystyle-\frac{1}{2}n(t)(\sigma_{1\bar{1}}-\sigma_{1\bar{2}}-\sigma_{\bar{1}\bar{2}})\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}]\rangle$
$\displaystyle\frac{{\rm d}}{{\rm
d}t}\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{2}V_{\bar{1}}^{\dagger}]\rangle=$
$\displaystyle-\frac{1}{2}n(t)[N_{c}(\sigma_{\bar{1}2}+\sigma_{12})-\frac{1}{N_{c}}\sigma_{1\bar{1}}]\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{2}V_{\bar{1}}^{\dagger}]\rangle$
$\displaystyle-\frac{1}{2}n(t)(\sigma_{1\bar{1}}-\sigma_{12}-\sigma_{\bar{1}2})\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}]\rangle$
(101)

This is a nice consistency check, as taking the system for four Wilson lines
(47) and letting $2\to\bar{2}$ reproduces the first and second of these
equations. Similarly, (47) with $\bar{2}\to 2$ reproduces the first and third.
Solving the first two gives

$\displaystyle\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}]\rangle$
$\displaystyle=N_{c}\,{\rm e}^{-C_{F}\int_{t_{2}}^{t}{\rm
d}s\,n(s)\sigma_{1\bar{1}}(s)}$
$\displaystyle\langle\operatorname{tr}[V_{1}V_{\bar{2}}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle$
$\displaystyle=(N_{c}^{2}-1)\,e^{-\frac{1}{2}\int_{t_{2}}^{t}{\rm
d}s\,n(s)(N_{c}(\sigma_{\bar{1}\bar{2}}+\sigma_{1\bar{2}})-\frac{1}{N_{c}}\sigma_{1\bar{1}})}+e^{-C_{F}\int_{t_{2}}^{t}{\rm
d}s\,n(s)\sigma_{1\bar{1}}(s)}\,.$ (102)

One can easily get
$\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{2}V_{\bar{1}}^{\dagger}]\rangle$
from the second of these equations by changing $\bar{2}\to 2$. The first
equation in (B.1) is a well known result, so it is nice that we reproduce
that.

However, this is nothing new, merely a check that the system of six Wilson
lines is consistent with the previous calculations. The remaining $3\times 3$
system contains the correlator we actually want to solve, but is also a more
complicated nonhomogeneous system. It is useful to define two vectors with the
3 unknown and 3 known functions

$\displaystyle\bm{C}_{1}^{\intercal}$
$\displaystyle=\Big{(}\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\operatorname{tr}[V_{2}V_{\bar{2}}^{\dagger}]\rangle,\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\operatorname{tr}[V_{2}V_{\bar{2}}^{\dagger}]\rangle,\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}V_{2}V_{\bar{2}}^{\dagger}]\rangle\Big{)}$
$\displaystyle\bm{C}_{2}^{\intercal}$
$\displaystyle=\Big{(}\langle\operatorname{tr}[V_{1}V_{\bar{2}}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle,\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{\bar{1}}^{\dagger}V_{2}]\rangle,\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}]\rangle\Big{)}\,.$
(103)

Then we can write the remaining system of differential equations as

$\frac{{\rm d}}{{\rm
d}t}\bm{C}_{1}=-\frac{1}{2}n(t)\left(\underline{\bm{\mathbb{M}}}_{1}\bm{C}_{1}+\underline{\bm{\mathbb{M}}}_{2}\bm{C}_{2}\right)\,.$
(104)

The $3\times 3$ matrices $\underline{\bm{\mathbb{M}}}_{1}$ and
$\underline{\bm{\mathbb{M}}}_{2}$ are subsets of the $6\times 6$ matrix
$\underline{\bm{\mathbb{M}}}$ and have the form

$\underline{\bm{\mathbb{M}}}_{1}=\begin{bmatrix}N_{c}(\sigma_{12}+\sigma_{\bar{1}\bar{2}}+\sigma_{2\bar{2}})-\frac{1}{N_{c}}\sigma_{1\bar{1}}&\sigma_{1\bar{1}}+\sigma_{2\bar{2}}-\sigma_{1\bar{2}}-\sigma_{\bar{1}2}&0\\\
\sigma_{12}+\sigma_{\bar{1}\bar{2}}-\sigma_{1\bar{2}}-\sigma_{\bar{1}2}&2(C_{F}\sigma_{1\bar{1}}+N_{c}\sigma_{2\bar{2}})&\sigma_{1\bar{2}}+\sigma_{\bar{1}2}-\sigma_{12}-\sigma_{\bar{1}\bar{2}}\\\
0&\sigma_{1\bar{1}}+\sigma_{2\bar{2}}-\sigma_{12}-\sigma_{\bar{1}\bar{2}}&N_{c}(\sigma_{2\bar{2}}+\sigma_{1\bar{2}}+\sigma_{\bar{1}2})-\frac{1}{N_{c}}\sigma_{1\bar{1}}\end{bmatrix}\,.$
(105)

$\underline{\bm{\mathbb{M}}}_{2}=\begin{bmatrix}\sigma_{1\bar{2}}-\sigma_{12}-\sigma_{2\bar{2}}&\sigma_{\bar{1}2}-\sigma_{2\bar{2}}-\sigma_{\bar{1}\bar{2}}&0\\\
0&0&-2\sigma_{2\bar{2}}\\\
\sigma_{\bar{1}\bar{2}}-\sigma_{2\bar{2}}-\sigma_{\bar{1}2}&\sigma_{12}-\sigma_{1\bar{2}}-\sigma_{2\bar{2}}&0\end{bmatrix}\,.$
(106)

This can be solved numerically for the three functions in $\bm{C}_{1}$, and
the result can be seen in figure 7.

Figure 7: The exact and large-$N_{c}$ solutions to the system of differential
equations (104).

####  B.1.1 Quark-gluon splitting in the large-$N_{c}$

As showed in section 4.3 all the functions in $\bm{C}$ can be solved exactly
in the large-$N_{c}$ limit. The two terms with highest powers of $N_{c}$,
$\langle[\operatorname{tr}[V_{1}V_{2}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{\bar{1}}^{\dagger}]\operatorname{tr}[V_{2}V_{\bar{2}}^{\dagger}]\rangle$
and
$\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\operatorname{tr}[V_{2}V_{\bar{2}}^{\dagger}]\rangle$
can be gotten directly from (71) and (4.3) respectively. Alternatively one can
count the $N_{c}$ powers in (B.1) and realize that
$\underline{\bm{\mathbb{M}}}_{1}$ and $\underline{\bm{\mathbb{M}}}_{2}$
simplify to

$\underline{\bm{\mathbb{M}}}_{1}\simeq\begin{bmatrix}N_{c}(\sigma_{12}+\sigma_{\bar{1}\bar{2}}+\sigma_{2\bar{2}})&0&0\\\
\sigma_{12}+\sigma_{\bar{1}\bar{2}}-\sigma_{1\bar{2}}-\sigma_{\bar{1}2}&N_{c}(\sigma_{1\bar{1}}+2\sigma_{2\bar{2}})&0\\\
0&\sigma_{1\bar{1}}+\sigma_{2\bar{2}}-\sigma_{12}-\sigma_{\bar{1}\bar{2}}&N_{c}(\sigma_{2\bar{2}}+\sigma_{1\bar{2}}+\sigma_{\bar{1}2})\end{bmatrix}\,.$
(107)

$\underline{\bm{\mathbb{M}}}_{2}\simeq\begin{bmatrix}0&0&0\\\ 0&0&0\\\
\sigma_{\bar{1}\bar{2}}-\sigma_{2\bar{2}}-\sigma_{\bar{1}2}&\sigma_{12}-\sigma_{1\bar{2}}-\sigma_{2\bar{2}}&0\end{bmatrix}\,.$
(108)

The solutions to the simplified differential equation leads to (3.1).

###  B.2 Eight Wilson lines

For more than six Wilson lines the matrix in (4.2) becomes so big that it is
impractical to analyze it by hand.

For eight lines it involves the $4!=24$ projections of
$\langle[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{1}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{2}}[V_{3}V_{\bar{3}}^{\dagger}]_{i_{3}j_{3}}[V_{4}V_{\bar{4}}^{\dagger}]_{i_{4}j_{4}}\rangle$,
and the matrix $\underline{\bm{\mathbb{M}}}$ has $24^{2}$ elements. The power
of our result in section 4 is here evident, as simply solving the differential
equation (4.2) for $K=4$ numerically immediately gives the result for eight
lines Wilson lines. To get the Wilson line correlators we want from (35) the
four last labels must be changed
$(3,\bar{3},4,\bar{4})\to(\bar{1},1,\bar{2},2)$ so that
$\langle[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{1}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{2}}[V_{3}V_{\bar{3}}^{\dagger}]_{i_{3}j_{3}}[V_{4}V_{\bar{4}}^{\dagger}]_{i_{4}j_{4}}\rangle\to\langle[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{1}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{2}}[V_{\bar{1}}V_{1}^{\dagger}]_{i_{3}j_{3}}[V_{\bar{2}}V_{2}^{\dagger}]_{i_{4}j_{4}}\rangle$.
The two relevant solutions are shown in figure 8.

Figure 8: The exact and large-$N_{c}$ version of
$C_{1\bar{1}2\bar{2}}=\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}]\operatorname{tr}[V_{2}V_{\bar{2}}^{\dagger}V_{\bar{1}}V_{1}^{\dagger}]\operatorname{tr}[V_{\bar{2}}V_{2}^{\dagger}]\rangle$
and
$C_{2\bar{1}\bar{2}2}=\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}V_{2}V_{\bar{2}}^{\dagger}V_{\bar{1}}V_{1}^{\dagger}V_{\bar{2}}V_{2}^{\dagger}]\rangle$.

One thing to notice in figure 8 is that for the case of eight Wilson lines
correlators, keeping only the first term in the large-$N_{c}$ limit does not
work well.

##  Appendix C Derivation of differential equation

In this appendix we will show in more detail how the differential equation
(4.2) was derived.

We start with the derivation of (4.2). To illustrate we will first show the
calculation for $K=2$, that is calculating
$\langle[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{1}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{2}}\rangle$.
This generalizes rather easily to the arbitrary $K$ case (59). Expanding the
first of these pairs like in (4.1) up to first order of $\epsilon$ and
defining ${\bm{A}}\equiv A^{a}t^{a}$ it becomes

$\displaystyle[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{1}}(t+\epsilon,t_{0})=[V_{1}V^{\dagger}_{\bar{1}}+ig\int^{t+\epsilon}_{t}{\rm
d}s({\bm{A}}_{1}(s)V_{1}V^{\dagger}_{\bar{1}}-V_{1}V^{\dagger}_{\bar{1}}{\bm{A}}_{\bar{1}}(s))$
$\displaystyle+\frac{1}{2}g^{2}\int^{t+\epsilon}_{t}{\rm
d}s\int^{t+\epsilon}_{t}{\rm
d}s^{\prime}(2{\bm{A}}_{1}(s)V_{1}V^{\dagger}_{\bar{1}}{\bm{A}}_{\bar{1}}(s^{\prime})-{\bm{A}}_{1}(s){\bm{A}}_{1}(s^{\prime})V_{1}V^{\dagger}_{\bar{1}}-V_{1}V^{\dagger}_{\bar{1}}{\bm{A}}_{\bar{1}}(s){\bm{A}}_{\bar{1}}(s^{\prime}))]_{i_{1}j_{1}}\,.$
(109)

Here all the Wilson lines on the right hand side go from $t_{0}$ to $t$. After
taking the medium average (7) and using the Fierz identity (44) the last term
becomes

$\displaystyle\frac{1}{2}g^{2}\int^{t+\epsilon}_{t}{\rm
d}s\int^{t+\epsilon}_{t}{\rm
d}s^{\prime}\langle[2{\bm{A}}_{1}(s)V_{1}V^{\dagger}_{\bar{1}}{\bm{A}}_{\bar{1}}(s^{\prime})-{\bm{A}}_{1}(s){\bm{A}}_{1}(s^{\prime})V_{1}V^{\dagger}_{\bar{1}}-V_{1}V^{\dagger}_{\bar{1}}{\bm{A}}_{\bar{1}}(s){\bm{A}}_{\bar{1}}(s^{\prime})]_{i_{1}j_{1}}\rangle$
$\displaystyle=\frac{1}{2}g^{2}n(t)\epsilon\left[\gamma_{1\bar{1}}\langle\operatorname{tr}(V_{1}V_{\bar{1}}^{\dagger})\rangle\delta_{i_{1}j_{1}}-(2C_{F}\gamma_{0}+\frac{1}{N_{c}}\gamma_{1\bar{1}})\langle
V_{1}V_{\bar{1}}^{\dagger}\rangle_{i_{1}j_{1}}\right]\,.$ (110)

Now adding the second pair of Wilson lines and taking the medium average,
while disregarding higher orders of $\epsilon$, it takes the form

$\displaystyle\langle[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{1}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{2}}\rangle(t+\epsilon,t_{0})=\langle[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{1}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{2}}\rangle(t,t_{0})$
$\displaystyle+\frac{1}{2}g^{2}\epsilon\,n(t)\left\langle\left[\gamma_{1\bar{1}}\operatorname{tr}(V_{1}V_{\bar{1}}^{\dagger})\delta_{i_{1}j_{1}}-(2C_{F}\gamma_{0}+\frac{1}{N_{c}}\gamma_{1\bar{1}})[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{1}}\right][V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{2}}\right.$
$\displaystyle+[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{1}}\left.\left[\gamma_{2\bar{2}}\operatorname{tr}(V_{2}V_{\bar{2}}^{\dagger})\delta_{i_{2}j_{2}}-(2C_{F}\gamma_{0}+\frac{1}{N_{c}}\gamma_{2\bar{2}})[V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{2}}\right]\right\rangle$
$\displaystyle-g^{2}\int^{t+\epsilon}_{t}{\rm d}s\int^{t+\epsilon}_{t}{\rm
d}s^{\prime}\big{\langle}\left[{\bm{A}}_{1}(s)V_{1}V_{\bar{1}}^{\dagger}-V_{1}V_{\bar{1}}^{\dagger}{\bm{A}}_{\bar{1}}(s)\right]_{i_{1}j_{1}}\left[{\bm{A}}_{2}(s^{\prime})V_{2}V_{\bar{2}}^{\dagger}-V_{2}V_{\bar{2}}^{\dagger}{\bm{A}}_{\bar{2}}(s^{\prime})\right]_{i_{2}j_{2}}\big{\rangle}\,.$
(111)

The last term simplifies to

$\displaystyle g^{2}\int^{t+\epsilon}_{t}{\rm d}s\int^{t+\epsilon}_{t}{\rm
d}s^{\prime}\big{\langle}\left[{\bm{A}}_{1}(s)V_{1}V_{\bar{1}}^{\dagger}-V_{1}V_{\bar{1}}^{\dagger}{\bm{A}}_{\bar{1}}(s)\right]_{i_{1}j_{1}}\left[{\bm{A}}_{2}(s^{\prime})V_{2}V_{\bar{2}}^{\dagger}-V_{2}V_{\bar{2}}^{\dagger}{\bm{A}}_{\bar{2}}(s^{\prime})\right]_{i_{2}j_{2}}\big{\rangle}$
$\displaystyle=\frac{1}{2}g^{2}n(t)\epsilon\big{\langle}\gamma_{12}[V_{1}V_{\bar{1}}^{\dagger}]_{i_{2}j_{1}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{1}j_{2}}+\gamma_{\bar{1}\bar{2}}[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{2}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{1}}$
$\displaystyle-\gamma_{1\bar{2}}\delta_{i_{1}j_{2}}[V_{2}V_{\bar{2}}^{\dagger}V_{1}V_{\bar{1}}^{\dagger}]_{i_{2}j_{1}}-\gamma_{\bar{1}2}\delta_{i_{2}j_{1}}[V_{1}V_{\bar{1}}^{\dagger}V_{2}V_{\bar{2}}^{\dagger}]_{i_{1}j_{2}}$
$\displaystyle+\frac{1}{N_{c}}(\gamma_{1\bar{2}}+\gamma_{\bar{1}2}-\gamma_{12}-\gamma_{\bar{1}\bar{2}})[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{1}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{2}}\big{\rangle}\,.$
(112)

Letting $\epsilon$ go to zero this turns into a differential equation

$\displaystyle\frac{{\rm d}}{{\rm
d}t}\langle[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{1}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{2}}\rangle$
$\displaystyle=\frac{1}{2}g^{2}n(t)\big{\langle}[\frac{1}{N_{c}}(\gamma_{12}+\gamma_{\bar{1}\bar{2}}-\gamma_{1\bar{2}}-\gamma_{\bar{1}2}-\gamma_{1\bar{1}}-\gamma_{2\bar{2}}-2(N_{c}^{2}-1)\gamma_{0})][V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{1}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{2}}$
$\displaystyle+\gamma_{1\bar{1}}\operatorname{tr}(V_{1}V_{\bar{1}}^{\dagger})\delta_{i_{1}j_{1}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{2}}+\gamma_{2\bar{2}}\operatorname{tr}(V_{2}V_{\bar{2}}^{\dagger})[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{1}}\delta_{i_{2}j_{2}}$
$\displaystyle-\gamma_{12}[V_{1}V_{\bar{1}}^{\dagger}]_{i_{2}j_{1}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{1}j_{2}}-\gamma_{\bar{1}\bar{2}}[V_{1}V_{\bar{1}}^{\dagger}]_{i_{1}j_{2}}[V_{2}V_{\bar{2}}^{\dagger}]_{i_{2}j_{1}}$
$\displaystyle+\gamma_{1\bar{2}}\delta_{i_{1}j_{2}}[V_{2}V_{\bar{2}}^{\dagger}V_{1}V_{\bar{1}}^{\dagger}]_{i_{2}j_{1}}+\gamma_{\bar{1}2}\delta_{i_{2}j_{1}}[V_{1}V_{\bar{1}}^{\dagger}V_{2}V_{\bar{2}}^{\dagger}]_{i_{1}j_{2}}\big{\rangle}\,.$
(113)

Now we only have to project out the two possible ways to connect the Wilson
lines. Contracting with $\delta_{j_{2}i_{1}}\delta_{j_{1}i_{2}}$ and
$\delta_{j_{1}i_{1}}\delta_{j_{2}i_{2}}$ gives $\frac{{\rm d}}{{\rm
d}t}\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}V_{2}V_{\bar{2}}^{\dagger}]\rangle$
and $\frac{{\rm d}}{{\rm
d}t}\langle\operatorname{tr}[V_{1}V_{\bar{1}}^{\dagger}]\operatorname{tr}[V_{2}V_{\bar{2}}^{\dagger}]\rangle$
respectively. In section 3 we wanted to calculate
$\langle\operatorname{tr}[V_{1}V_{2}^{\dagger}V_{\bar{2}}V_{\bar{1}}^{\dagger}]\rangle$,
which is similar to the above, but not exactly the same. Fortunately, our
choice of labels is just a convention, and completely arbitrary. Simply making
the three changes $\bar{1}\to 2$, $2\to\bar{2}$ and $\bar{2}\to\bar{1}$ turns
(C) into the system of differential equations in (47). The difference in this
approach compared to what we did in section 4.1 is that (C) contains both of
(4.1) and (4.1). This compact form is highly convenient when considering more
than four Wilson lines. Generalizing the steps from equation (C) to (C) to an
arbitrary number $K$ pairs of Wilson lines produces the differential equation
(4.2).

The next step is to show how to get from Eq. (4.2) to the matrix elements (63)
and (64). Any pair of Wilson lines has two free indices. Take for example the
second Wilson line pair in (4.2) which is $W_{i_{2}j_{2}}^{2}$. Start by
projecting out these two indices in all the ways possible, and at the same
time making as few assumptions as possible about the rest of the Wilson lines.
It turns out that projecting out with two Kronecker deltas gives all the
information we need. There are also only two possibilities that need to be
considered: either $W^{2}$ can connect to other Wilson lines, or it connects
to itself and becomes a trace. These two possibilities are given by projecting
with $\delta_{j_{1}i_{2}}\delta_{j_{2}i_{3}}$ and
$\delta_{j_{1}i_{3}}\delta_{j_{2}i_{2}}$, respectively. To use Wilson lines 1,
2 and 3 is arbitrary. These labels can be changed to anything else without
changing the result, so the calculation is completely general.

Using (4.2) and projecting out by the two deltas
$\delta_{j_{1}i_{2}}\delta_{j_{2}i_{3}}$ gives a differential equation for
$\langle[W^{1}W^{2}W^{3}]_{i_{1}j_{3}}W_{i_{4}j_{4}}^{4}\dots
W_{i_{K}j_{K}}^{K}\rangle$.

$\displaystyle\frac{{\rm d}}{{\rm
d}t}\langle[W^{1}W^{2}W^{3}]_{i_{1}j_{3}}W_{i_{4}j_{4}}^{4}\dots
W_{i_{K}j_{K}}^{K}\rangle$
$\displaystyle=-\frac{1}{2}g^{2}n(t)\left(N_{c}(\sigma_{\bar{1}2}+\sigma_{\bar{2}3})+\frac{1}{N_{c}}A_{K}\right)\langle[W^{1}W^{2}W^{3}]_{i_{1}j_{3}}W_{i_{4}j_{4}}^{4}\dots
W_{i_{K}j_{K}}^{K}\rangle$
$\displaystyle-\frac{1}{2}g^{2}n(t)(\sigma_{\bar{1}3}+\sigma_{2\bar{2}}-\sigma_{23}-\sigma_{\bar{1}\bar{2}})\langle\operatorname{tr}W^{2}[W^{1}W^{3}]_{i_{1}j_{3}}W_{i_{4}j_{4}}^{4}\dots
W_{i_{K}j_{K}}^{K}\rangle$ $\displaystyle+(\dots)\,.$ (114)

The $(\dots)$ in the end are terms that are not completely determined by the
projection that was made.

Next up is the case where we project out with
$\delta_{j_{1}i_{3}}\delta_{j_{2}i_{2}}$, making a differential equation for
$\langle\operatorname{tr}W^{2}[W^{1}W^{3}]_{i_{1}j_{3}}W_{i_{4}j_{4}}^{4}\dots
W_{i_{K}j_{K}}^{K}\rangle$.

$\displaystyle\frac{{\rm d}}{{\rm
d}t}\langle\operatorname{tr}W^{2}[W^{1}W^{3}]_{i_{1}j_{3}}W_{i_{4}j_{4}}^{4}\dots
W_{i_{K}j_{K}}^{K}\rangle$
$\displaystyle=-\frac{1}{2}g^{2}n(t)\left(N_{c}(\sigma_{\bar{1}3}+\sigma_{\bar{2}2})+\frac{1}{N_{c}}A_{K}\right)\langle\operatorname{tr}W^{2}[W^{1}W^{3}]_{i_{1}j_{3}}W_{i_{4}j_{4}}^{4}\dots
W_{i_{K}j_{K}}^{K}\rangle$
$\displaystyle-\frac{1}{2}g^{2}n(t)(\sigma_{\bar{2}3}+\sigma_{2\bar{1}}-\sigma_{23}-\sigma_{\bar{1}\bar{2}})\langle[W^{1}W^{2}W^{3}]_{i_{1}j_{3}}W_{i_{4}j_{4}}^{4}\dots
W_{i_{K}j_{K}}^{K}\rangle$ $\displaystyle+(\dots)\,.$ (115)

In the notation from section 4 these equations become

$\displaystyle\frac{{\rm d}}{{\rm d}t}C_{23m_{3}\dots m_{N}}$
$\displaystyle=-\frac{1}{2}g^{2}n(t)\left(N_{c}(\sigma_{\bar{1}2}+\sigma_{\bar{2}3})+\frac{1}{N_{c}}A_{K}\right)C_{23m_{3}\dots
m_{N}}$
$\displaystyle-\frac{1}{2}g^{2}n(t)(\sigma_{\bar{1}3}+\sigma_{2\bar{2}}-\sigma_{23}-\sigma_{\bar{1}\bar{2}})C_{32m_{3}\dots
m_{N}}$ $\displaystyle+(\dots)\,,$ (116) $\displaystyle\frac{{\rm d}}{{\rm
d}t}C_{32m_{3}\dots m_{N}}$
$\displaystyle=-\frac{1}{2}g^{2}n(t)\left(N_{c}(\sigma_{\bar{1}3}+\sigma_{\bar{2}2})+\frac{1}{N_{c}}A_{K}\right)C_{32m_{3}\dots
m_{N}}$
$\displaystyle-\frac{1}{2}g^{2}n(t)(\sigma_{\bar{2}3}+\sigma_{2\bar{1}}-\sigma_{23}-\sigma_{\bar{1}\bar{2}})C_{23m_{3}\dots
m_{N}}$ $\displaystyle+(\dots)\,.$ (117)

Both of these equations are consistent with the matrix elements (63) and (64).
The point is that when all the indices are projected out all the Wilson lines
will connect in one of these two ways. Either they will connect to other
Wilson lines or they will only connect to themselves. And since we have shown
that in either way the resulting expression is given by (63) and (64) it means
that these two equations are correct for all the possible combinations.

## References

  * (1) D. d’Enterria, “Jet quenching,” Landolt-Bornstein, vol. 23, p. 471, 2010\.
  * (2) A. Majumder and M. Van Leeuwen, “The Theory and Phenomenology of Perturbative QCD Based Jet Quenching,” Prog. Part. Nucl. Phys., vol. 66, pp. 41–92, 2011.
  * (3) Y. Mehtar-Tani, J. G. Milhano, and K. Tywoniuk, “Jet physics in heavy-ion collisions,” Int. J. Mod. Phys. A, vol. 28, p. 1340013, 2013.
  * (4) 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, vol. 757, pp. 102–183, 2005.
  * (5) K. Adcox et al., “Formation of dense partonic matter in relativistic nucleus-nucleus collisions at RHIC: Experimental evaluation by the PHENIX collaboration,” Nucl. Phys. A, vol. 757, pp. 184–283, 2005.
  * (6) K. Aamodt et al., “Suppression of Charged Particle Production at Large Transverse Momentum in Central Pb-Pb Collisions at $\sqrt{s_{NN}}=$ 2.76 TeV,” Phys. Lett. B, vol. 696, pp. 30–39, 2011.
  * (7) V. Khachatryan et al., “Charged-particle nuclear modification factors in PbPb and pPb collisions at $\sqrt{s_{\mathrm{N}\;\mathrm{N}}}=5.02$ TeV,” JHEP, vol. 04, p. 039, 2017.
  * (8) S. Chatrchyan et al., “Observation and studies of jet quenching in PbPb collisions at nucleon-nucleon center-of-mass energy = 2.76 TeV,” Phys. Rev. C, vol. 84, p. 024906, 2011.
  * (9) G. Aad et al., “Observation of a Centrality-Dependent Dijet Asymmetry in Lead-Lead Collisions at $\sqrt{s_{NN}}=2.77$ TeV with the ATLAS Detector at the LHC,” Phys. Rev. Lett., vol. 105, p. 252303, 2010.
  * (10) B. Abelev et al., “Measurement of charged jet suppression in Pb-Pb collisions at $\sqrt{s_{NN}}$ = 2.76 TeV,” JHEP, vol. 03, p. 013, 2014\.
  * (11) R. Baier, Y. L. Dokshitzer, S. Peigne, and D. Schiff, “Induced gluon radiation in a QCD medium,” Phys. Lett. B, vol. 345, pp. 277–286, 1995\.
  * (12) R. Baier, Y. L. Dokshitzer, A. H. Mueller, S. Peigne, and D. Schiff, “Radiative energy loss and p(T) broadening of high-energy partons in nuclei,” Nucl. Phys. B, vol. 484, pp. 265–282, 1997.
  * (13) 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, vol. 483, pp. 291–320, 1997.
  * (14) R. Baier, Y. L. Dokshitzer, A. H. Mueller, and D. Schiff, “Radiative energy loss of high-energy partons traversing an expanding QCD plasma,” Phys. Rev. C, vol. 58, pp. 1706–1713, 1998.
  * (15) B. G. Zakharov, “Fully quantum treatment of the Landau-Pomeranchuk-Migdal effect in QED and QCD,” JETP Lett., vol. 63, pp. 952–957, 1996.
  * (16) B. G. Zakharov, “Radiative energy loss of high-energy quarks in finite size nuclear matter and quark - gluon plasma,” JETP Lett., vol. 65, pp. 615–620, 1997.
  * (17) U. A. Wiedemann, “Gluon radiation off hard quarks in a nuclear environment: Opacity expansion,” Nucl. Phys. B, vol. 588, pp. 303–344, 2000.
  * (18) Y. Mehtar-Tani and K. Tywoniuk, “Radiative energy loss of neighboring subjets,” Nuclear Physics A, vol. 979, p. 165–203, Nov 2018.
  * (19) Y. Mehtar-Tani, “Gluon bremsstrahlung in finite media beyond multiple soft scattering approximation,” JHEP, vol. 07, p. 057, 2019.
  * (20) P. Caucal, E. Iancu, A. H. Mueller, and G. Soyez, “Vacuum-like jet fragmentation in a dense QCD medium,” Phys. Rev. Lett., vol. 120, p. 232001, 2018.
  * (21) P. Caucal, E. Iancu, A. H. Mueller, and G. Soyez, “Nuclear modification factors for jet fragmentation,” JHEP, vol. 10, p. 204, 2020.
  * (22) F. Domínguez, J. G. Milhano, C. A. Salgado, K. Tywoniuk, and V. Vila, “Mapping collinear in-medium parton splittings,” The European Physical Journal C, vol. 80, Jan 2020.
  * (23) J.-P. Blaizot, F. Dominguez, E. Iancu, and Y. Mehtar-Tani, “Medium-induced gluon branching,” Journal of High Energy Physics, vol. 2013, Jan 2013.
  * (24) L. Apolinário, N. Armesto, J. G. Milhano, and C. A. Salgado, “Medium-induced gluon radiation and colour decoherence beyond the soft approximation,” JHEP, vol. 02, p. 119, 2015.
  * (25) A. Kovner and U. A. Wiedemann, “Eikonal evolution and gluon radiation,” Physical Review D, vol. 64, Oct 2001.
  * (26) F. Dominguez, C. Marquet, B.-W. Xiao, and F. Yuan, “Universality of Unintegrated Gluon Distributions at small x,” Phys. Rev. D, vol. 83, p. 105005, 2011.
  * (27) P. Arnold, “Landau-Pomeranchuk-Migdal effect in sequential bremsstrahlung: From large-$N$ QCD to $N$=3 via the SU($N$) analog of Wigner 6-$j$ symbols,” Phys. Rev. D, vol. 100, no. 3, p. 034030, 2019.
  * (28) B. G. Zakharov, “Color randomization of fast gluon-gluon pairs in the quark-gluon plasma,” J. Exp. Theor. Phys., vol. 128, no. 2, pp. 243–258, 2019.
  * (29) Y. Hatta and T. Ueda, “Non-global logarithms in hadron collisions at $N_{c}$ = 3,” Nucl. Phys. B, vol. 962, p. 115273, 2021.
  * (30) J. Jalilian-Marian and Y. V. Kovchegov, “Inclusive two-gluon and valence quark-gluon production in DIS and pA,” Phys. Rev. D, vol. 70, p. 114017, 2004.  [Erratum: Phys.Rev.D 71, 079901 (2005)].
  * (31) E. Iancu and D. N. Triantafyllopoulos, “Higher-point correlations from the JIMWLK evolution,” JHEP, vol. 11, p. 105, 2011.
  * (32) T. Lappi, H. Mäntysaari, and A. Ramnath, “Next-to-leading order Balitsky-Kovchegov equation beyond large $N_{c}$,” Phys. Rev. D, vol. 102, no. 7, p. 074027, 2020.
  * (33) F. Dominguez, C. Marquet, A. M. Stasto, and B.-W. Xiao, “Universality of multiparticle production in QCD at high energies,” Phys. Rev. D, vol. 87, p. 034007, 2013.
  * (34) Z. Nagy and D. E. Soper, “Parton shower evolution with subleading color,” JHEP, vol. 06, p. 044, 2012.
  * (35) K. Hamilton, R. Medves, G. P. Salam, L. Scyboz, and G. Soyez, “Colour and logarithmic accuracy in final-state parton showers,” 11 2020.
  * (36) M. Gyulassy, P. Levai, and I. Vitev, “Reaction operator approach to nonAbelian energy loss,” Nucl. Phys. B, vol. 594, pp. 371–419, 2001.
  * (37) X.-N. Wang and M. Gyulassy, “Gluon shadowing and jet quenching in A + A collisions at s**(1/2) = 200-GeV,” Phys. Rev. Lett., vol. 68, pp. 1480–1483, 1992.
  * (38) P. Aurenche, F. Gelis, and H. Zaraket, “A Simple sum rule for the thermal gluon spectral function and applications,” JHEP, vol. 05, p. 043, 2002\.
  * (39) Y. Mehtar-Tani and K. Tywoniuk, “Improved opacity expansion for medium-induced parton splitting,” Journal of High Energy Physics, vol. 2020, Jun 2020\.
  * (40) J. a. Barata, Y. Mehtar-Tani, A. Soto-Ontoso, and K. Tywoniuk, “Revisiting transverse momentum broadening in dense QCD media,” 9 2020.
  * (41) T. Altinoluk, N. Armesto, G. Beuf, M. Martınez, and C. A. Salgado, “Next-to-eikonal corrections in the cgc: gluon production and spin asymmetries in pa collisions,” Journal of High Energy Physics, vol. 2014, Jul 2014.