Wilson line correlators beyond the large-$N_c$

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.


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 [1][2][3].
Experiments at RHIC [4,5] and the LHC [6][7][8][9][10] 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., [11][12][13][14][15][16][17]. 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 [18][19][20][21]. 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 → 2 splittings.
Previous studies of such processes focused mostly on a hard photon splitting into a quark-antiquark pair [18,22] 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 where the angular brackets denote an average over medium configurations. For a generic 1 → 2 process, the amplitude squared can be reduced to a product of two-, three-and fourpoint correlators [23,24]. To calculate these processes it is imperative to know the form of the Wilson line correlator appearing on the right hand side of (1.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 tr which scales like N K c . 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 ∼ 1/N 2 c 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 [25] for a technique based on diagonalization of the evolution matrix and [24,26] for an iterative procedure. The derivation of the evolution matrix culminates in Eq. (4.26). 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.1) above, the generic expansion in N c takes the following form, 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 sub-leading Nc 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 [23,24], 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 [27] and to understand hadronization after exiting the QGP [28]. 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. (2.5), which is also employed in other contexts than for a thermal medium, see, e.g., [29] 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 [25,30,31]. Furthermore, sub-leading corrections in color have also been considered in the context of high-energy QCD evolution at next-to-leading order [32]. Finally, the generic color structure of high-energy QCD events is actively studied [26,33]. 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., [34,35].
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. (1.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.

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 − , p), where p + = (p 0 + p 3 )/2 is identified with the LC energy E ≡ p + , p − = p 0 − p 3 is negligible and p is the transverse momentum. The parton interacts with the medium, which is modelled by a classical background gauge field A µ,a (t, 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 [11][12][13][14] and Zakharov [15,16], 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) [36] approximation.
It is possible to construct Feynman rules from the BDMPS-Z approach, with special in-medium propagators and vertices [18]. In this formulation a highly energetic parton travelling through the medium can be described by the propagator In this expression V R is a Wilson line in the representation R, which is given by where the symbol P enforces path ordering. A quark transforms in the fundamental representation, so the group generator is T a F ≡ t a ij . Similarly, a gluon transforms in the adjoint representation, and its group generator is T a A ≡ (T a ) bc = −if abc . The final results in this paper will mainly concern fundamental lines, which we will denote by V ≡ V F . Similarly, we will write the adjoint lines as U ≡ V A . Focusing on fundamental lines is sufficient, since one can always transform adjoint Wilson lines to fundamental ones through the identity In the absence of interactions, i.e. when the Wilson line is evaluated at g = 0, we simply get 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. (2.1). 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 . . . , and we assume that the correlator of the medium fields takes the form which corresponds to the Gaussian noise approximation. Here, n(t) is the (time-dependent) density of scattering centers in the medium and 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 The form of the function σ depends on how the medium is modelled. The two main ways of calculating this is through the Gyulassy-Wang model [37] or through Hard Thermal Loop theory [38]. These models differ mainly in how infrared screening is implemented when q ⊥ → 0. In this paper, we will however work in the harmonic oscillator approximation, which accounts for multiple soft interactions. In this case, the potential σ(r) can be cast as whereq 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 = N 2 c −1 2Nc and C A = N c , respectively. In this paper we will useq =q F unless otherwise stated. In Eq. (2.9) 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. [19,39,40].
We stress that the approximation in (2.8) 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 → 2 splitting processes in the medium, where the initial particle has energy E and the two splitting products carry, respectively, ω 1 = (1 − z)E and ω 2 = zE. This is formally equivalent to setting the energy of the mother particle, E → ∞, and considering a finite momentum sharing fraction 0 z 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 (x|G(t, t 0 )|x 0 ), in configuration space, for E (t − t 0 ) −1 gets strongly constrained to the classical path connecting the initial and final transverse positions, see Eq. (2.1), and leads to where the classical trajectory is given by . This corresponds to the product of a Wilson line, trailing the direction of the particle, times a vacuum propagator, see Eq. (2.4). Corrections to this limit can also be systematically be calculated [41]. In the mixed representation, this leads to, where n = p/E, see [18,22]. The last term in this product is simply the Fourier transform of the vacuum propagator. In detail, the 1 → 2 partonic processes we consider are: 1) γ → q +q, 2) q → q + g, 3) g → 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 → q +q, but we do not explicitly evaluate the spectrum in this case. All three processes consist of one (off-shell) particle 1 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 r 1 (t) (r1(t)) and r 2 (t) (r2(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 where n 1 ≡ p 1 zE and n 2 ≡ 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 (r 1 ) ≡ V 1 and γ(r 1 − r2) ≡ γ 12 , etc.
Finally, in the harmonic approximation, we need the square of the differences of the transverse coordinates. Using the eikonal approximation this is where we have assumed that the angle θ is small.

Emission spectra
In this section we will present the results for the in-medium emission spectra dI dzdθ for the inmedium 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 where α can be α em or α 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 [22] dI full dz dθ = dI vac dz dθ + dI med dz dθ = dI vac dz dθ (1 + F med (z, θ)) . The term F med (z, θ) contains the medium modification to the processes. For a generic medium profile the medium radiation reads 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 ⊥ = z(1 − z)Eθ is comparable to the transverse momentum accumulated in the medium Q s ∼qL [23]. Finally, the factors C (n) (t b , t a ) appearing in (3.3) are n-particle correlators that have support during time t a < t < t b , and t f = 2 z(1−z)Eθ 2 is the formation time of the process. The splitting process is illustrated in Fig. 1.
For a medium with fixed density and extension L, we haveq(t) =q Θ(L − t). In this case, the integrals over the emission times t 1 and t 2 in (3.3) can be split, so that Taking into account that the Wilson line correlators are real the medium modification term can be written [22] (3.6) 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 θ to map the regions where medium-induced corrections appear, as quantified by the factor F med (z, θ). Our focus here is to provide a test bed for evaluating precisely the multi-Wilson line correlators appearing in (3.6), and we will therefore not worry about the validity of the eikonal approximation (2.10) of the splitting products. Including non-eikonal corrections on the particle trajectories will be postponed to future work.

Derivation of the splitting functions
Photon splitting We will start with the case of a photon splitting into a quark-antiquark pair, i.e. γ → q +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 (3.1) with the QED coupling constant α em and the Altarelli-Parisi splitting function being where n f is the number of active flavors. Furthermore, the correlator C (3) reduces to an effective two-point function because the photon does not carry color charge. We have 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 C qγ (t, t 1 ) = S 12 (t, t 1 ), is generally referred to as a dipole correlator, and is known to be where r = r 1 − r 2 is the difference of transverse positions of the two Wilson lines. For a fixed separation, i.e. r = const., in the HO approximation and in a medium with constant density, it simply reads S 12 (t, t 1 ) = e − 1 4q (t−t 1 ) r 2 , where, as a reminder, we have denoted q ≡q F . However, for the kinematics we consider, see Eq. (2.13), this becomes Then, assuming that t = t 2 and t 2 − t 1 ∼ t f , we find with ω = z(1−z)E. This implies that medium modifications appear, in this term, whenever ω 3 θ 4 q. The correlator of four Wilson lines 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. (4.10). In the large-N c limit, however, it can be calculated analytically through the simplified differential equation Eq. (4.12). 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 , note the absence of explicit normalization factors at this stage. In the large-N c approximation, the first of these correlators reads simply, where τ ≡ t 2 − t 1 , and the dependence on t 1 appears as a consequence of the fixed trajectories. Here we have used the eikonal (2.12) and harmonic oscillator approximations (2.8).
Similarly, at large-N c , the second correlator is Only the latter of these correlators appears in the spectrum, cf. Eq. (3.7), but we include both for completeness. 2 Assuming the dominance of the first term in (3.13), setting t = L and assuming that L t f , we find that where we put ξ ≈ 2/3. The factor in the exponential becomes large whenever ωθ < √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 S 12 (t 2 , t 1 ) or quadrupole C2 1 (L, t 2 ). For a fixed energy ω, the two conditions are equal at the critical angle Let us also define the characteristic energies ω d = (q/θ 4 ) 1/3 and ω broad = √q L/θ. At large angles θ > θ c , the condition from the dipole starts affecting soft gluon emissions, i.e. ω d < ω broad . This reflects the length-dependence color coherence. On the one hand, the dipole, which has support only during the formation time t f 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 C   Fig. 2a, 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 aŝ q = 1.5 GeV 2 /fm, θ = 0.5 and z = 0.5. We notice the fast decay of the dipole, that goes like ∼ e −t 3 according to (3.10), compared to the exponential decay of the quadrupole, i.e. ∼ e −t , at large times. Finally, we notice that the large-N c approximation to the full quadrupole, given in Eq. (3.13), 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 [22], but not calculated explicitly. For this process, the vacuum emission spectrum is given by (3.1), with the QCD coupling constant α s and the Altarelli-Parisi splitting function P gq (z) = C F . The four-and three-point functions read The emission spectrum is composed of correlators of two, four and six Wilson lines. The three-point function can be solved exactly, see (B.7), resulting in This expression is very similar to the dipole term in Eq. (3.10) 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 1 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.13).

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 [23]. For this process, the vacuum emission spectrum is given by (3.1) with the QCD coupling constant α s and the Altarelli-Parisi splitting function In this case the 4-and 3-point functions read When cast as correlators of Wilson lines in the fundamental representation, the 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.26) or by writing it in terms of adjoint Wilson lines (A.17). In the end, the result reads Note the similarity to the previous results, see Eqs. (3.10) and (3.18). The 4-point function consists of two different correlators of eight Wilson lines. They can be calculated through the differential equation in Eq. (4.26). Interestingly, the fourpoint function C (4) gg involves a eight-point correlator, see the second term in (3.20), 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.

and the correlators read
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.

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 θ. In Fig. 2, we show how C ij (t, t 1 ), with blue, solid curves, and C 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 chooseq = 1.5 GeV 2 /fm, θ = 0.5 and z = 0.5.
While this approximation turns out work extremely well for the photon splitting, see Fig. 2a, we note that it has a more limited range of applicability for both the quark-gluon, see Fig. 2b, and gluon-gluon, see Fig. 2c, splitting processes, respectively. In all of the cases the exact value is slightly higher than the approximate one. As we derived analytically,  ij terms all decay as ∼ e −q(t−t 1 ) 3 τ 2 θ 2 f (z) , where f (z) is a process dependent regular function. The C (4) ij terms are more complicated, especially at early times where all terms contribute, but at late times the dominant contribution comes from ∼ e −q(t−t 2 )θ 2 .
The ratio of double-differential in-medium to vacuum spectrum reveals the medium modification factor F med (z, θ) = dI med /(dzdθ) dI vac /(dzdθ). We plot this factor, calculated at finite N c , for the three processes in Fig. 3. These results have been obtained for the medium parametersq = 1.5 GeV 2 /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 med (z, θ) has roughly the same characteristic shape for all three processes. The medium modifications appear at large angles θ > θ c , in between the characteristic lines ω 3 θ 3 <q and ω 2 θ 2 <qL 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 f < L. In fact, we can recast these conditions in terms of the formation time of the process, namely t f < t d and t f < t broad , where (3.25) (c) Gluon-gluon splitting. The modifications appear for the range of formation times t broad < t f < t d and θ > θ c [22]. 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 ≈ 2 larger for gluon splitting than for the photon. 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 med (z, θ)| large−Nc /F med (z, θ) 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 struc-  ture, 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 f ∼ t broad and t f ∼ t d , but at the same time t f < L. This is most clearly seen in the gluon-gluon splitting, cf. Fig. 4c. 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 dI dzdθ for three different splitting processes, and shown that the resulting expressions factorize into threeand 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 broad , and decoherence during the formation of the splitting, t d . These were identified first in [22] 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.

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.

Four Wilson lines
The simplest Wilson line correlator comes from the pair production process (3.7), where there is a trace of four Wilson lines tr 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 + . Then, following [25], we expand them between t and t + to get 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ˆt where in the first step we applied the medium average (2.5). Then, keeping terms up to the first order of this becomes Using the Fierz identity this results in the differential equation , both going from times t 2 to an arbitrary time t. The grey lines at the beginning and end indicate the colour connections.
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 tr[V 1 V † 2 V2V † 1 ] . 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 tr[V 1 V † 2 ] tr[V2V † 1 ] that appeared in equation (4.5). The inclusion of this term in the differential equation (4.5) just represents the possibility for color rotation to happen at each time.
To continue one can find a complementary differential equation for tr[V 1 V † 2 ] tr[V2V † 1 ] and see if we can find a solution for the set. Going through the same procedure as above gives We now have a set of two coupled differential equations. To save space the following notation will be used C 12 . This notation warrants some more explanation. Both of these expressions are composed of the two pairs of Wilson lines, namely V 1 V † 2 and V2V † 1 . 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 C2 1 means that . 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.5) and (4.6) can be gathered into the following system, d dt where the evolution matrix takes the following form, (4.8) Here we have used Eq. (2.7) to define σ 12 = σ(r 1 − r 2 ), and introduced To proceed, we employ the the harmonic approximation (2.8). For the eikonal, straight-line trajectories, given in Eqs. (2.12), the evolution matrix becomes 10) where we have defined τ ≡ t 2 −t 1 , ξ = z 2 +(1−z) 2 and assumed that the angle between the two particles θ 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 [22] calculated the four-point function tr[V 1 V † 2 V2V † 1 ] 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 (4.7) and count the powers of N c in each term in the evolution matrix and the vector of correlators, taking into account that C 12 ∼ N 2 c and C2 1 ∼ N 1 c . In this limit we also have C F ∼ N c /2. The terms on the right-hand side of (4.7) then have the following powers of N c , . (4.11) 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 (4.11). We see that the the next-to-leading power of N c turns out to be a factor N −2 c 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 d dt Now it is evident that the differential equation for C 12 is separable and can be solved easily, which means that C2 1 also can be solved. This leads to the equations (3.12) and (3.13).
The physical picture of this differential equation is quite transparent. The correlator of (orange, solid and orange, dashed lines, respectively). We also plot only the leading, diagonal term of the large-N c approximation of C2 1 (orange, dotted line) which exhibits the correct large-time asymptotic behavior. 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 (4.8), and is associated with a factor of ∼ σ, which scales as N −1 c . Say you start in the state C 12 shown on the right in figure 5, scaling as N 2 c . If you exchange a gluon you pick up a factor N −1 c from the σ, and go to the state C2 1 , which is a single trace correlator that scales as N 1 c , so in total this transition is associated with a factor N 0 c . This is a factor N −2 c smaller compared to the starting point so it can safely be dropped in the large-N c limit. However, starting with C2 1 and going to C 12 you go from a state that scales as N 1 c to one scaling as N 2 c , but you lose a power of N c from the σ, 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.
The solutions to (4.10) and (4.12) 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 C2 1 correlator. At late times, we observe an exponential suppression, ∝ e −t , with a slope that is in good agreement with the first term of Eq. (3.13). At early times, there is an interplay between C2 1 and C 12 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.13).
The C 12 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 (4.12) and some smaller correction term, (1) , (4.13) where C = C 12 , C2 1 is a vector of the correlators in question, so that C (1) is a factor O(N −2 c ) smaller than C (0) . We can also write the matrix M in a form that isolates the large-N c terms from the finite-N c corrections, i.e. M = M (0) + M corr. , (4.14) where the first term strictly corresponds to the leading terms in the large-N c limit. In our example above, we find that 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 −2 c . The correlators at leading color, i.e. C (0) , are known. They solve the simplified set of equations dC (0) (t)/dt = − n(t) 2 M (0) C (0) (t), and are given explicitly in (3.12) and (3.13). This can now be used to calculate the color sub-leading contributions C (1) . Simply plugging this into the full differential equation (4.7) results in the following differential equation for the first correction where we have neglected terms that are even more sub-leading, i.e. resulting from M corr. C (1) . This is an nonhomogeneous version of the large-N c system of differential equations (4.12), and can also be solved exactly. As an example the first correction to C 1 12 (t) is (4.18) The first correction contains C (0) 21 (t), given in (3.13) 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 12 (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

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 (4.7), and will show how to easily calculate all the matrix elements in the K! × 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 Note that one can divide the correlators of these Wilson lines into pairs on the form [V n V † m ] injm 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 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 (4.19) to whatever is needed in the specific problem at hand. For example, choosing K = 3 and changing labels (1,1, 2,2, 3,3) → (1, 2,2,1, 2,2) gives the structure needed in (3.16), while K = 4 and changing labels (1,1, 2,2, 3,3, 4,4) → (1,1, 2,2,1, 1,2, 2) reproduces the correlators in (3.20). 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 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, 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.21) and projecting out with . . 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. 3 Although it is possible to use (4.21) 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. (4.7). Making use of the notation we described above we want to write the system of differential equations for K pairs of Wilson lines as where p 1 p 2 . . . p K also is one of the K! permutations of 12 . . . K. Starting from (4.21), one can deduce the general form of the matrix elements M p 1 p 2 ...p K m 1 m 2 ...m 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 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 (4.24) for p 1 p 2 . . . p K being any other permutation of m 1 m 2 . . . m K . This means that out of the K! 2 matrix elements, only 1 2 K!(K 2 − K + 2) are non-zero. These are given by the relatively simple formulas (4.23) and (4.24). Putting it all together this becomes (4.26) Of course, for all differential equations you need to specify some initial conditions. It is clear from the definition of the Wilson line (2.2) that V ij (t 0 , t 0 ) = δ ij . The trace of this is tr V (t 0 , t 0 )=N c . This means that the initial condition of a 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 C M , where the superscript M is meant to indicate that this scales as N M c . Then, Eq. (4.26) can be represented as (4.27) -23 - The first matrix contains the diagonal elements, written in detail in (4.23). The second matrix represents the non-diagonal elements, and σ is a block containing non-zero elements, which we get from Eq. (4.24).

Wilson line correlators in the large-N c limit
The system of differential equations given by Eq. (4.26) is to our knowledge not possible to solve analytically in the case where σ 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. (4.27).
Since σ ∼ N −1 c , the diagonal matrix elements go as ∼ N 0 c + N −2 c , while the non-diagonal ones go as ∼ N −1 c . Multiplying in the vector on the end and representing every term by its N c scaling this becomes (4.28) 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 −2 c compared to the leading term. Translating this back to the form in Eq. (4.27), this becomes d dt . . .
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 ∼ N M c . 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 ∼ N M +1 . 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 (4.29).
It is clear from this discussion that the system of differential equations (4.26) simplifies, in the large-N c limit, to (4.30) 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 M +1 c 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.30) the differential equation for the correlator scaling as N K c is Since this is exactly solvable, 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  This can also be written in terms of dipoles, namely S m kk (t, s) . (4.36) 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.36) is a recursive relation ("bootstrap") that stops when you reach the term with K traces, which is given in terms of dipoles in (4.34). Since the only Wilson line structure that appears in both (4.34) and (4.36) is dipoles, it means all the correlator can be written in terms of dipoles. In Ref. [33] 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 −2 c relative to the leading terms, following the steps in Eqs. (4.13) and (4.14), and below. We have nevertheless not pursued this program further in this work.

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.26). This system can be solved numerically. We showed that in the large-N c limit the resulting simplified system of differential equations, Eq. (4.30), can be solved exactly. We also provided a general way to compute color subleading corrections, suppressed by N −2 c relative to the leading terms. This was done in detail for the four-point correlator, in Eqs. (4.13) and (4.14), 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 → 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 θ, but the differences become bigger as θ 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 −2 c . 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 protonproton collisions. Finally, we plan on extending the formulation to account for non-eikonal corrections to the particle trajectories.

A Calculation of spectrums
Here we will show the calculations leading up to the for the emission spectra dI dzdθ . The Feynman rules from [18] 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 [22] but we will restate some of the results. The amplitude is where the photon-quark vertex is given by The initial hard process is represented by the amplitude M 0 . After using the eikonal approximation (2.11) this becomes (up to some phase that cancels when we take the square) We have used the more compact notation to write 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 (3.3) with (3.7) and (3.8).

A.2 Quark-gluon splitting
The amplitude was calculated in [18] and is where the quark-gluon vertex is Again this simplifies in the eikonal limit (2.11) We have denoted the adjoint Wilson line as V A (r 2 ) ≡ U 2 . Squaring the amplitude, summing/averaging over spins and colors and taking the medium average gives where the relevant Altarelli-Parisi splitting function is To continue we transform the adjoint Wilson lines into fundamental ones using the identity (2.3). The resulting expression will contain many group generators t a , and can be simplified by using the Fierz identity (4.4). Finally, completely in the fundamental representation the Wilson line structure becomes .
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 (3.3) with (3.16) and (3.17).

A.3 Gluon-gluon splitting
The calculation of the emission spectrum for gluon-gluon splittings was done in [23]. For completeness we will also include the main results here. The matrix element of the process is where the gluon-gluon vertex is (A.12) In the eikonal approximation (2.11) the amplitude is where (m 1 , m 2 , m 3 ) now is some permutation of (1,2, 2). Thee non-zero non-diagonal entries are given by The 3 × 3 matrices M 1 and M 2 are subsets of the 6 × 6 matrix M and have the form (B.10) This can be solved numerically for the three functions in C 1 , and the result can be seen in figure 7. The solutions to the simplified differential equation leads to (3.19).

B.2 Eight Wilson lines
For more than six Wilson lines the matrix in (4.26) becomes so big that it is impractical to analyze it by hand. For eight lines it involves the 4! = 24 projections of [ ] i 4 j 4 , and the matrix M has 24 2 elements. The power of our result in section 4 is here evident, as simply solving the differential equation (4.26) for K = 4 numerically immediately gives the result for eight lines Wilson lines. To get the Wilson line correlators we want from (3.20) the four last labels must be changed (3,3, 4,4) → (1, 1,2, 2) so that The two relevant solutions are shown in figure 8.
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.
Letting go to zero this turns into a differential equation Now we only have to project out the two possible ways to connect the Wilson lines. Contracting with δ j 2 i 1 δ j 1 i 2 and δ j 1 i 1 δ j 2 i 2 gives d dt tr[ respectively. In section 3 we wanted to calculate tr[V 1 V † 2 V2V † 1 ] , 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 changes1 → 2, 2 →2 and2 →1 turns (C.5) into the system of differential equations in (4.7). The difference in this approach compared to what we did in section 4.1 is that (C.5) contains both of (4.5) and (4.6). This compact form is highly convenient when considering more than four Wilson lines. Generalizing the steps from equation (C.1) to (C.5) to an arbitrary number K pairs of Wilson lines produces the differential equation (4.21).
The next step is to show how to get from Eq. (4.21) to the matrix elements (4.23) and (4.24). Any pair of Wilson lines has two free indices. Take for example the second Wilson line pair in (4.21) which is W 2 i 2 j 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 δ j 1 i 2 δ j 2 i 3 and δ j 1 i 3 δ 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.21) and projecting out by the two deltas δ j 1 i 2 δ j 2 i 3 gives a differential equation The (. . . ) 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 δ j 1 i 3 δ j 2 i 2 , making a differential equation (C.7) In the notation from section 4 these equations become Both of these equations are consistent with the matrix elements (4.23) and (4.24). 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 (4.23) and (4.24) it means that these two equations are correct for all the possible combinations.