Medium-induced gluon radiation and colour decoherence beyond the soft approximation

We derive the in-medium gluon radiation spectrum off a quark within the path integral formalism at finite energies, including all next-to-eikonal corrections in the propagators of quarks and gluons. Results are computed for finite formation times, including interference with vacuum amplitudes. By rewriting the medium averages in a convenient manner we present the spectrum in terms of dipole cross sections and a colour decoherence parameter with the same physical origin as that found in previous studies of the antenna radiation. This factorisation allows us to present a simple physical picture of the medium-induced radiation for any value of the formation time, that is of interest for a probabilistic implementation of the modified parton shower. Known results are recovered for the particular cases of soft radiation and eikonal quark and for the case of a very long medium, with length much larger than the average formation times for medium-induced radiation. Technical details of the computation of the relevant $n$-point functions in colour space and of the required path integrals in transverse space are provided. The final result completes the calculation of all finite energy corrections for the radiation off a quark in a QCD medium that exist in the small angle approximation and for a recoilless medium.


Introduction
The dense and hot state of matter produced in heavy ion collisions, commonly referred to as a quark-gluon plasma (QGP), is characterised by the deconfinement of quarks and gluons up to distances much larger than the size of hadrons. Thus, it affords a window of opportunity to study QCD in a regime usually not accessible perturbatively. Given the very short lifetime of the QGP, only probes generated within, as part of the overall collision, can be used to determine its properties. Of particular importance are hard probes, produced in a hard scattering, which carry information about the initial stages of the collision and can be addressed within perturbation theory. The large centre-of-mass collision energy per nucleon pair in collider experiments -the Relativistic Heavy Ion Collider (RHIC) at BNL and, above all, the Large Hadron Collider (LHC) at CERN -leads to abundant hard particle production and, consequently, to the possibility of measuring a variety of observables with high statistics. Among such hard probes, those related to the modification of jets and jet-like properties of particle production resulting from the effects imparted by the hot and dense medium to the propagation dynamics of high-energy particles -what is commonly referred to as jet quenching, see e.g. [1,2] -provide the opportunity to extract detailed information about the medium. Manifestations of jet quenching phenomena have been experimentally observed both for jet-like observables -e.g., in single inclusive particle spectra [3][4][5][6][7] and two-particle correlations [8][9][10] -and for fully reconstructed jets in PbPb collisions at the LHC [11][12][13][14][15][16][17][18][19][20][21] (see also [22,23] for related results at RHIC). In particular, the observation of a sizeable increase, with respect to the proton-proton case, of the energy asymmetry in dijet systems without modification of their azimuthal correlation and for which the missing energy is recovered at large angles away from the jet axis in the form of soft particles, appeared, at first sight, to pose serious challenges for the standard explanation of jet quenching in terms of medium-induced gluon radiation in which energy loss and broadening are intrinsically connected. While several phenomenological explanations [24][25][26][27][28][29][30][31][32][33] have since been put forward, none provides a complete quantitative description of the data nor proper theoretical justification for some of the assumptions. This fact stresses the importance of having a faithful description of the mechanisms of energy loss so that a successful comparison with data can be done and the properties of the QGP deconvoluted from the confounding effects present in heavy-ion collisions. Such a requirement is particularly relevant for Monte Carlo implementations. While analytical models are currently derived within the high-energy approximation and usually implement energy-momentum conservation a posteriori, Monte Carlo generators encode it by construction. This implies further assumptions whick lack a firm theoretical basis and often extend the use of theoretical models beyond their strict validity region [34]. In order to avoid this situation, finite energy corrections for the full radiation spectrum must be correctly (and fully) accounted for. Efforts in this direction include [35][36][37]. In this work we extend the results derived in [36], where an interpolation function between the soft and hard limits for gluon radiation off a quark was deduced, to account for an arbitrary momentum fraction to be carried by the radiated gluon.
The path-integral formalism [38] is used throughout. We allow transverse motion of all particles in the emission process, thus relaxing the common assumption that only the softest particle is allowed such movement. In an earlier work where the same constraint was removed [37], the gluon radiation off a gluon was computed assuming small formation times for the radiated gluon. By going beyond this approximation we not only present the full result of medium-induced gluon radiation at finite formation time and finite energy, but also a clear physical picture in which the relevant color coherence effects are identified. A future goal will be to convert this physical picture into a probabilistic approach suitable for Monte Carlo implementation. Let us also mention that while the region of validity of small formation times can be quite comfortable, it demands a large medium to be crossed by the radiating partons. This approximation will not hold for the ones that are emitted close to the edge of the medium.
Let us anticipate the main results of our work. We provide the most general form of a double inclusive spectrum, for gluon radiation off a fast quark in a QCD medium in the small angle approximation and for a recoilless medium. From it and in order to get expressions that can be used in practice, we make use of the Gaussian approximation for medium averages, work in the limit of large number of colours and employ the approximation of many soft scatterings, known as multiple soft or harmonic oscillator approximation. One surprising effect we find is that the kinematic part of the splitting -terms which come from the hard vertex -is modified with respect to the vacuum and cannot be factorized as an Altarelli-Parisi splitting function for the case of finite formation time. The physical origin of this effect, whose relevance for applications remains to be studied in the future, is the different histories that the quark in amplitude and conjugate amplitude follow.
Our main qualitative result is, however, the physical picture in terms of colour coherence/decoherence of the radiation process, with a factorization of dipole cross sections and a colour decoherence parameter with the same physical origin, but with some differences in the mathematical form, as the one found in the studies of the antenna radiation, ∆ med [39][40][41][42][43]. We find that colour coherence between the outgoing quark and gluon survives longer than in those calculations done in the eikonal limit for all propagators but the one of the softest particle, therefore suppressing the spectrum of radiated gluons. All these corrections are found to vanish in the limit of negligible formation time of the produced gluon, where we recover the results known in the literature, in particular those in [37].
The paper is organized as follows: In section 2, the formalism used to describe the inmedium propagation of partons is introduced, and the different contributions to the singlegluon emission spectrum for a static colour medium profile are calculated. The average over all possible colour configurations of the medium and generalisation of the decoherence parameter is carried out in section 3. The total radiation spectrum with the evaluation of the necessary path-integrals is presented in section 4, while the final conclusions are presented in section 5. Technical details are provided as appendices.

Quasi-eikonal in-medium parton propagation
The time scale involved in the propagation of energetic partons is much smaller than the characteristic time for changes in the configuration of the medium that they traverse. This difference in time scales allows for the computation of the parton-medium interaction to be performed for a fixed, but arbitrary, medium configuration and, at a later stage, for the ensemble of medium configurations to be accounted for through an averaging procedure.
The multiple scattering of the propagating parton off medium components is mediated by the exchange of gluons with typical, purely transverse, momenta of the order of the characteristic medium scales. As a result, the otherwise eikonal trajectory of the partonthe rotation of its colour phase without degradation of its (large) longitudinal momentum -is perturbed by Brownian motion in the transverse plane. The in-medium propagation of a parton with light-cone 1 plus momentum p + from transverse position x i at time x i+ (where its colour is α i ) to transverse position x f at time x f + , with colour rotated to α f , is given by the Green's function where the Wilson line accounts for the colour rotation resulting from an arbitrary number of scatterings off the medium field A − ≡ A a − T a (T a being the colour matrix in the corresponding representation), while the free propagator encodes the random walk in the transverse plane. The Wilson line W α f α i in eq. (2.2), and consequently G α f α i in eq. (2.1), should be understood to carry colour indices in the relevant representation for the parton under consideration. In the following, fundamental colour indices, as relevant for propagating quarks, will be written in uppercase latin letters, while for the gluon the adjoint indices will be written in lowercase latin letters. For compactness, and improved readability, we introduce the shorthand notation

Amplitudes
To compute the radiation of a gluon off an energetic quark produced in a hard process in the early stages of a heavy-ion collision, two separate contributions to the amplitude ought to be considered: the case in which the splitting occurs outside a fixed length medium (see figure 1a) and thus only the initial quark experiences medium interactions; and the complementary situation in which the splitting occurs within the medium boundaries (x 0+ , L + ) (see figure 1b) and the interaction of all partons with the medium must be accounted for 2 .
(a) q → qg splitting where only the initial particle interacts with the medium.  Taking into account the dominance of the plus component of the initial momentum ( / p p + γ − ) and the preservation of the longitudinal light-cone momentum of the radiated gluon through the propagation in the medium (k 1+ = k + ⇒ k ·k 1 = 0), the total amplitude can be written as where the out and in contributions are given respectively by We recall that the hard process, of amplitude M h , from which the quark originates is unmodified by the surrounding environment since it occurs within a time/length scale too small to be resolved by the medium. and with X 0 = (x 0+ , x 0 ) the coordinates of the quark at the beginning of the medium, X 1 = (x 1+ , x 1 ), Y 1 = (x 1+ , y 1 ) and Z 1 = (x 1+ , z 1 ) the coordinates at the emission point for the initial quark, final quark and gluon 3 respectively, and X = (L + , x) , Y = (L + , y) , Z = (L + , z), the coordinates after the final scatterings.

Emission cross section
The cross section for single gluon emission is given by the average over the ensemble of medium configurations . . . (to be carried out in section 3) with dΩ k = (2π) −3 dk dk + /(2k + ) and analogously for dΩ q . This inelastic cross sectionthe squared amplitude averaged over initial colour and summed over final spin, colour and gluon polarisation -is given by After simplification of both the Dirac algebra -for which details are given in appendix A -and the colour algebra -explicitly carried in appendix B.1 -the separate contributions on the right hand side of eq. (2.9) read 4 (2.10) 3 Although we use different emission coordinates, these are in fact the same (x1 = y1 = z1). 4 For compactness, we will use throughout the shorthands +∞ −∞ dx ≡ x and where P g←q (ζ) = C F [1 + (1 − ζ) 2 ]/ζ is the (vacuum) Altarelli-Parisi splitting function [44] with ζ the fraction of longitudinal momenta carried by the gluon (as defined in appendix A). It was necessary to introduce new coordinates for the complex conjugate amplitude, to represent the initial quark (final quark, gluon). The emission point in the complex conjugate amplitude is denoted by x 2+ . The complete set of coordinates for the three contributions in eq. (2.9) is shown in figure 2 assuming It should be noted that one of the transverse momenta from the Dirac structure of the in − out term, and all the ones in the in − in term, correspond to an internal component in the T in amplitude. Therefore, they should be written as the derivative of the initial transverse coordinate of the corresponding Green's function at the vertex: (2.13) As mentioned in section 1, remarkably, it is not possible to simplify the number of kinematic terms in eq. (2.12) as in [37] -we are only able to do it if we assume that the formation time (i.e. the difference in plus light-cone coordinates of the gluon emission point in the amplitude and in the complex conjugate amplitude, τ f orm = x 2+ − x 1+ ) are small. Therefore, the initial broadening in the amplitude and in the complex conjugate amplitude, and consequently the corresponding internal momenta are different (this difference vanishes with decreasing τ f orm ). This implies that the kinematic part of the in-medium splitting depends on medium properties through the formation time. See details in Appendix A.
We further note that the (vacuum) Altarelli-Parisi splitting function is recovered not only in the out − out term, but also in the in − out and in − in contributions once we take q = q i = −k = −k i .

Average over the medium ensemble
The squared matrix elements -eqs. (2.10), (2.11), and (2.12) -have been calculated for a fixed arbitrary medium colour configuration in which the propagation and emission    (2.12). In each plot the full arrows represent the quarks and the dashed arrows the gluons. Black arrows hold for the amplitude and red arrows for the complex conjugate amplitude. The medium is shown by a coloured region that starts at x 0+ and ends at L + . processes take place. In order to account for the ensemble of possible medium colour configurations an averaging procedure is carried out. Since all medium information is encoded in the Wilson lines W , where all colour information resides, the averaging procedure amounts to the evaluation of correlators of Wilson lines (several examples are listed in appendix C). The longitudinal locality of the interactions with the medium, a consequence of the high-energy approximation, allows for the decomposition of the longitudinal support of each squared matrix element into regions with a constant number of Wilson lines. Recalling that the colour structure and transverse momentum dynamics (i.e. the random walk in transverse plane) are factorised, see eq. (2.1), together with the convolution properties of the G's, it is possible to separate a propagator at an arbitrary point within a given longitudinal support. That is to say, that a Green's function with longitudinal support in the interval (x i+ , x f + ) can be separated at a longitudinal location

Separation into regions
Using equation (3.1), it is possible to identify, depending on the number of Wilson lines propagating simultaneously, a single in-medium region for the out − out contribution, two regions for the in−out and three different regions for the in−in term, as shown in figure 3.
For the out − out term, eq. (2.10), no separation is necessary as the only configuration that propagates through the medium is that involving two Wilson lines (for the initial quark in both amplitude and complex conjugate amplitude). Omitting kinematical terms, where the longitudinal support for the average is explicitly shown as a subscript.
For the in − out contribution we identify two separate regions (see figure 3b): from the production point x 0+ to the gluon emission point x 1+ , where there are just two fundamental Wilson lines, and from x 1+ to the end of the medium L + where an additional adjoint Wilson line is present. To account for this separation, the Wilson line that comes from the complex conjugate amplitude, represented as a red arrow in figure 4, has to be separated as in equation (3.1) with the introduction of additional transverse coordinates. The new transverse coordinates are represented in red, while the remaining transverse structure in black. Explicit calculations are shown in appendix B.2, but hereon, we will restrict ourselves to the large-N limit. Omitting kinematical terms, the result reads  : Schematic representation of the transverse structure for the amplitude T in T † out . The transverse ending points are written in grey, while the transverse coordinates to be integrated out in black and red. The red coordinates correspond to the propagators that were divided into different regions (note that z 1 = y 1 = x 1 ).
Proceeding analogously for the in − in term, 3 regions can be defined (figure 3c). A schematic representation of the transverse coordinate structure of this expression is provided in figure 5. The complete result is derived in appendix B.2 but, in the large-N limit 5 , reads The transverse ending points are represented in grey, while the transverse coordinates to be integrated out in black and red. The red coordinates correspond to the propagators that were divided into different regions (note that z 1 = y 1 = x 1 and z 2 = y 2 = x 2 ).
Focusing now only on the coloured part of equations (3.2), (3.3) and (3.4), we expand the Wilson line up to second order in the medium fields. In particular, the result for the 2-point function is (eq. (C.9))

6)
5 Note that we assume x1+ < x2+ here. Since the result for x1+ > x2+ is the complex conjugate of the one for x1+ < x2+, to get the final result we will multiply by 2 and take the real part.
where n(ξ + ) is the longitudinal density of scattering centers in the medium and σ the cross section for scattering of the dipole formed by the particles located at transverse coordinates x and y with the medium: Using the results derived in appendix C.2, region II of the in − in contribution can be factorised at large N into the independent average of two dipoles formed by the final quark and gluon, while region III (see appendix C.3) into a dipole times an independent quadrupole. Thus (3.8)

Dipole approximation
For an opaque media, the dipole cross section can be approximated by its small distance component [45,46], whereq, the transport coefficient, characterises the typical squared transverse momentum that the particle acquires, per mean free path λ, from the interaction with the medium. This approximation, alternatively referred to as multiple soft scattering approximation or dipole approximation, is valid for small transverse distances r. Although the medium is expanding, we will perform, for simplicity, the calculations for an homogeneous static medium for whichq is a constant 6 . The result for the 4-point correlation function is explicitly derived in appendix C.2 in the Gaussian approximation for field correlators. For large N and in region III, equation (C.30) can be written (see the notation in appendices C.1 and C.2, and figure 5 for the coordinates) as (3.10) 6 An expanding medium can be accounted for by a change of variables [38,47,48].
with the colour decoherence parameter given by  With this factorisation (3.10), a friendly interpretation of the colour structure for every contribution to the total production cross section can be obtained, see figure 6. Region I corresponds to the Brownian motion of the initial quark before it emits. In region II, corresponding to the gluon formation time τ f orm , the final quark and gluon are colour correlated -this is easy to see as the spectrum is proportional to two traces of Wilson lines in this region 7 . Region III, on the other hand, admits a nice interpretation as the independent broadening of the quark and the produced gluon times a decoherence parameter controlling the probability that the gluon decorrelates in colour from the quark. Only in this case of decorrelation the spectrum exists, the opposite possibility being suppressed.
Let us elaborate more on this important point now. The ∆ coh parameter defined in (3.11) that controls the decoherence of colour sources contains the same physics of colour decoherence as the the one, named ∆ med , previously found in the antenna radiation [39,40,42,43] 8 . Note that the exponent in ∆ coh given in (3.12) is proportional to the distance between the quark and the gluon. The medium-induced spectrum is suppressed for the case ∆ coh → 0, that is when the quark and the gluon remain in a colour coherent 7 Remember that the trace of two Wilson lines is proportional to the S-matrix of the process, rather than to the scattering amplitude, hence to the probability of no interaction -this is a standard result, see e.g. [38]. 8 Indeed, our ∆ coh would provide corrections to the prefactor in the direct contributions to the antenna spectrum.
state. In the opposite limit, ∆ coh → 1, the quark and the gluon lose their coherence and appear as two independent particles. In this sense, a medium-induced radiation in which the quark and the gluon remain in a colour coherent state is exponentially suppressed. For the particular case in which x 2+ L + a suppression factor of the form τ f orm /L + appears as shown previously in [37] 9 . Equation (3.11) shows clearly that this approximation deteriorates more and more when the position of the splitting vertex gets closer and closer to the end of the medium.
One additional comment is in order: both in the BDMPS limit with strictly eikonal quark lines (r 3 =r 3 ) or in the limit of hard gluon emissions [36] with strictly eikonal gluon lines (w 3 =w 3 ), ∆ coh = 1 and the independent broadening of the q and g, given by the traces in (3.10), happens instantaneously with no colour interference between quark and gluon in region III.

Vacuum gluon spectrum
In the previous section, all possible colour medium configurations were taken into account by averaging over the whole ensemble of possible colour profiles. The results -eqs. (3.2), (3.3) and (3.8) -are for some general trajectories in the transverse plane, r = r(ξ), that change with the propagation time ξ. In order to account for the transverse broadening of the propagating particles, it is necessary to integrate over all possible trajectories that each particle may undergo. Particular examples are computed in appendix D using a semiclassical approximation. This method, that is able to provide an exact solution for some cases (see appendix D.1), assumes that the dominant contribution to the path integral is mainly given by the classical trajectory plus local fluctuations at the end points of the trajectory. The classical path is determined taking into account the kinetic term from the propagator G 0 , eq. (2.3), and the potential term from the effective action of the medium scattering centres encoded in the several n-point correlation function.
Using these results from appendix D.1, in particular (D.28), and after Fourier transforming from coordinate space to momentum space, the out − out term can be finally written as 9 This can be seen in the following way: the coherent piece, called non-factorisable in [37], can be written , with the exponential in m22 (in m11) giving the independent broadening of quark and gluon (the colour coherence propagation of the quark-gluon system). The latter suppresses exponentially τ > x2+ +τ f orm , while m12 can be expressed as the product of gradients of the former. Such product of gradients is proportional to 1/(L+ − τ ). where the integral over b comes from the fact that we are using plane waves, and withq F = C Fq and ∆ξ 1 = L + − x 0+ . This expression, normalised to one, simply provides the momentum broadening of the initial quark that propagates through a medium. The corresponding number of emitted gluons is where σ el is the total elastic cross section. Using the same assumptions as in section 2.2, it is possible to calculate the elastic channel that is schematically represented in figure 7. The result reads Finally, one finds that where all undetermined factors cancel. Eq. (4.5) admits a simple interpretation as the typical vacuum emission for an off-shell quark which, however, has experienced broadening while traversing the medium. This fact can be more clearly seen by taking the limitq F L + → 0 in which the vacuum spectrum is recovered. In this limit, the integration over the quark phase space can be performed, and fixing the initial transverse direction p 0 = 0 one indeed obtains the expected expression with the corresponding Altarelli-Parisi splitting function 10

Medium emission spectrum
The remaining two terms (in − in and in − out) can be identified with the medium contribution 11 . Neglecting undetermined factors coming from the use of plane waves, the result reads where Σ i are the results of the path integrals for each region i = I, II, III (see appendix D for definitions, in particular eqs. (D.28), (D.37) and (D.43)). Note that the factor C F (explicit in the T in T † out contribution) should be approximated by its large-N limit, N/2. Performing the integration over x 0 andx 0 , Σ 1 can be simplified and eq. (4.7) results in (4.8) - 16 -In this work, we provide a complete calculation of the medium-induced single gluon radiation off a quark, in the regime where partons undergo multiple soft scatterings with the medium. The kinematic setup was extended beyond the eikonal limit by associating Brownian perturbations in the transverse plane to all propagating particles, thus extending our previous results [36] obtained in the limit of hard emitted gluons. We consider multiple scatterings of the rescattered partons and not only single scatterings as in [35]. We go beyond the work [37] by considering a finite size medium and, thus, interference effects with the vacuum radiation. This allowed us to recover the vacuum gluon radiation spectrum, with a complete factorisation of the corresponding splitting function. Our computation includes all finite-energy corrections that exist in the small angle approximation (i.e. the emission and deflections angles are small) and for static scattering centres (i.e. we do not consider recoil). We also provide technical details of the computation of the relevant npoint functions in colour space in the Gaussian approximation, and of the required path integrals in transverse space.
The final results, that include the evaluation of several n−point correlation functions and resolution up to four path-integrals, are presented in the large-N limit and using the Gaussian approximation for field correlators, but without constraints on the gluon formation time. The resulting spectrum is, therefore, a generalisation of previous works, describing the radiation spectrum off a non-eikonal quark. We find that the kinematic factors of the in-medium splitting show a dependency on medium properties that vanishes in the limit of negligible formation time, recovering in this situation the results in [37] 12 .
We confirm the finding in [37] that parton branching, in the limit of a very opaque medium, can be understood as a factorisation of single gluon emissions, where the total radiation spectrum is just an incoherent sum of each independent emitter but suppressed by the interferences with the vacuum radiation -as already derived in previous calculations of soft gluon emissions. As the parton shower continues its development, the medium starts to become less opaque and coherence effects between the final particles, fully included in our calculation, become important. In this case, the emission spectrum is additionally suppressed with respect to the factorised regime. This suppression is controlled by a noneikonal decoherence parameter, ∆ coh , that accounts for the broadening of the particles and contains the same physics of colour decoherence as in the previous results in the antenna [39][40][41][42][43]. This fact implies that most of the energy lost by a parton must occur earlier in its development than under the assumption that the fully factorised regime holds during the full development of the shower.

A Dirac algebra
In order to perform the Dirac algebra in equations (2.10), (2.11) and (2.12), we recall the completeness relation for Dirac spinors, where the massless case will be considered. Furthermore, in the light-cone gauge A + = 0, the gluon polarizations sum to where η = (0, 1, 0). At this point we make the kinematics explicit: that of the in − out term as and the in − in case as 1 4 spin,polū The origin of the kinematical combination (ζq i − (1 − ζ) k i ) seen in the above equation is the following: it is always possible to find a reference frame where the outgoing transverse momenta are opposite. For that, we consider a rotation operation around the x axis and three generic vectors for massless partons, in the original frame, given by: with p = k + q. The corresponding rotation matrix reads where the last equality comes from considering the eikonal approximation, p ⊥ , k ⊥ , q ⊥ E, with p 2 ⊥ = p 2 x + p 2 y and the same for k ⊥ and q ⊥ (|p| = p ⊥ ). Therefore p z = E, k z = ζE and q z = (1 − ζ)E. Designating p , k and q the 3-vectors in the rotated reference frame, one gets Choosing the reference frame where p = 0, With this rotation -identical for the amplitude and complex conjugate amplitudeit is possible to simplify one of the mentioned kinematical combinations in equation (A.6), but another one remains unless one assumes, as done in [37], that q 1 + k 1 = q 2 + k 2 . The difference between these quantities is proportional to the formation time, τ f , of the emitted gluon, that in [37] is taken to be negligible.

B.1 Before region separation
To perform the colour algebra of equations (2.10), (2.11) and (2.12), we explicitly separate the path integral part of the propagators from its Wilson line part (i.e. the transverse Brownian motion from the colour rotation): rewrite the adjoint Wilson lines in terms of fundamental ones (see e.g. [52]) via and make use of the Fierz identity Using these expressions, we get that the colour contribution to the out − out piece reads that the in − out one reads and that the in − in one reads (B.6)

B.2 After region separation
To simplify further the propagator structure of eq. (2.11), we use the convolution property of the Green's function (eq. (3.1)) to write the propagator of the initial quark in the complex conjugate amplitude as where the coordinates are explicitly shown in figure 4. Moreover, due to the locality of the medium averages, the only possible contraction of two fundamental indices, A and B, when only two Wilson lines are present in a given local interval, is given by .

(B.8)
Thus, from eq. (B.5), one can write As for the structure of eq. (2.12), assuming that x 1+ < x 2+ , we need to split the final quark and gluon Green's functions in the amplitude as follows (coordinates are shown in figure 5): Using eqs. (B.7) and (B.10), eq. (B.6) can be written (B.11)

C.1 Two-point correlation function
The simplest average to be computed is the one involving two fundamental Wilson lines. That is, Expanding the Wilson line up to second order in the medium field (see [52]), we find Using this relation, we get

(C.3)
A diagrammatic interpretation of these terms can be found in [38]. Fourier transforming the fields and using the fact that all scattering centres are Lorentz contracted in a plane located at x + , we get with |a − (q)| 2 being a general screened potential (usually taken as a Yukawa potential). In order to perform the average over all possible colour configuration we have to integrate over the transverse and longitudinal coordinates of the scattering centres, (x i+ , x i ). Doing so, we find where we have introduced the longitudinal density of scattering centres Thus, The dipole cross section is identified as The result can be re-exponented -with account of the ordering in the x + coordinate -due to the fact that there is only one possibility for colour state in this average. We finally find and analogously for the adjoint colour representation 13 ,

C.2 Four-point correlation function
The structure that we want to calculate is present in eqs. (3.3) and (3.4) 14 . We write (C.11) To perform the medium average we expand the Wilson lines up to second order in the medium fields (see e.g. [52]), like done in eq. (C.2), but only for the first + light-cone position τ . Doing so, we can write where V (L + , τ ; x) denotes the Wilson line from the position τ to the final extension of the medium and the correlator between two medium fields. We use the Gaussian approximation in which all information is contained in the two-point function. In the following, we will not write explicitly the + coordinates.
Then the operator T αβµν jkno in (C.11) reads (C.14) In the following, we will not work with B's but with v's that are related to the dipole cross section: We define the following vectors: that are not orthogonal: so the scalar product is defined through the matrix One can prove that Therefore we can write the operator in the following matrix form: .

(C.21)
Now we need to act repeatedly with the operator M = T − I on u 1 (as we have expanded close to this vector that is the one at initial times, (C.12) and (C.14)), and, in order to close the traces, we must project onto u 2 by doing u 2 GM n T u 1 (the + coordinate increases from right to left i.e. u 1 → u 2 ). We write where the matrix elements can be read from (C.21). At leading N , we get and analogously for m n−1−i 11 (τ,x + ) , and m ij (τ ) indicating the integrand in (C.15) evaluated at a given x + = τ . Using the notation as analogous to a path-ordered exponential between x + and y + , and taking into account that · · · , (C. 26) we can make the sum over n to get ∞ n=0 u 2 GM n T u 1 = N e N m 11 that is the expression that one would get by making the infinitesimal expansion (C.12) at late times (i.e. close to u 2 ) as we will do for the six-point function in the next sub appendix (note that m 12 = m 21 + m 11 − m 22 ). Let us note that this result is perfectly compatible with those obtained in the strict eikonal limit (with fixed transverse coordinates all along the trajectories of the colour charges) either using the same method [55] or diagonalising the matrix (C.22) as e.g. in [52]. Indeed, considering that the m ij (τ ) do not depend on τ , one gets the result at large N that reads .

(C.31)
Finally, these results can be applied to the colour structure in region II. In this case, the structure that we want to simplify is the same as in eqs. (C.32) Here, we must consider the repeated action of M on u 1 , with the change of notation and the simplifications produced by the fact that two coordinates are equal. The result reads (C. 33) At large N , the result reads: . (C.34)

C.3 Six-point correlation function
For the purposes of the main body of this work, it suffices to evaluate the six-point correlation function appearing in (3.4) in the large N limit. Following the procedure already adopted for the four-point case, we write and define the colour contractions u 1 = δ jk δ no δ rs , u 2 = δ js δ rk δ no , u 3 = δ jo δ nk δ rs , u 4 = δ jo δ rk δ ns , with scalar products given through the matrix Using the same technique discussed in appendix C.2 but now making the infinitesimal expansion (C.12) at late times, the matrix operator that expresses the medium averages in this basis in the Gaussian approximation reads Then, to perform the medium average in region III in (3.4), we compute n insertions of the matrix M = T − I 6×6 (see e.g. [56]), project onto the corresponding vectors and sum over n i.e. 40) with the + coordinate increasing from left to right i.e. u 6 → u 1 in this case (see the comment below (C.30)).
Successive insertions of M can be understood as either propagating a given colour structure u i (in the form of an insertion of diagonal matrix elements m ii ) or as swapping colour structure from u i into u j (with off-diagonal matrix element m ji ). As each diagonal matrix element (we will discuss their 1/N corrections later on) carries a factor of N and the projection of the leftmost colour structure (given by the leftmost of swaps) onto u 6 contributes with N p , p = 1, 2, 3, according to (C.37), a term in (C.40) with n insertions of which s are colour swaps will carry an overall power of N given by N n−s+p which we write as N p−s N n .
The leading order result is obtained from terms with no swaps, s = 0, (the diagonal propagation of the colour structure u 1 ) for which the final projection is u 6 ·u 1 = N 2 (p = 2), and terms with one swap, s = 1, m 61 from u 1 to u 6 for which the final projection yields u 6 · u 6 = N 3 (p = 3): .
It is straightforward to see that the subleading corrections m ii in the diagonal matrix elements can be recovered through the substitution m ii → m ii + m ii /N 2 .

D Path integrals in the multiple soft scattering approximation
The path integrals that appear can be solved analytically for very few examples. Here we use the semi-classical approximation [57,58], which consists in considering the path integral as the classical action. As a consequence, the trajectory of the particle inside the medium is considered to be the classical path in the exponent, while fluctuations are considered in the norm.
The semi-classical method provides an exact solution for some cases e.g. the free particle or the harmonic oscillator that, in order to clarify subsequent computations, we elaborate as a first step.

D.1 Semi-classical method and some general examples
The trajectory of free particles is described by the following propagator: withṙ = dr(ξ)/dξ. In the semi-classical approximation, we identify the above expression with G 0 (y + , y; where the classical action R cl is defined as 15 Note the common factor v(w3 −w3) in m11 and m66 and, in order to see the exact equivalence with (C.30), make the substitutions m11 → m22 and m66 → m11, and keep in mind that m61 = m21+m11−m22 = m12. and the Lagrangian is Using the Euler-Lagrange equations, we find so the classical trajectory r(ξ) is given by a straight line: The solution for the path integral (2.3) is [57,58] where D is the number of dimensions, i, j ∈ {1, . . . , D} and R cl is evaluated at the classical path r = l(ξ) and integrated over ξ ∈ [x + , y + ]. The derivatives are taken on the initial and final transverse coordinates. In our case, D = 2 and the result of the action is (D.9) So, one finally finds Another example is the harmonic oscillator: K osc (y + , y; x + , x|p + ) = 1 N Tr G(y + , y; x + , x|p + )W † (y + , x + , 0) = r(y + )=y where it was used the multiple soft scattering approximation (eq. (3.9)). The corresponding Lagrangian is: that results in the following equation of motion: r + Ω 2 r = 0, (D. 13) with imaginary frequency: The solution is Plugin this solution into the classical action, we find It is easy to show that in this case: Thus, the solution for the harmonic oscillator is the well-known result

D.2 Region I: Two path integrals
The region I is shared by the tree contributions to the total spectrum, eq.(2.5), only changing the ending coordinates (see figure 3 and eqs. (3.2), (3.3) and (3.8)): where ξ 1 ∈ [x 0+ , x 1+ ]. In this expression, both s 1 ands 1 enter the kinematical and potential terms. Performing the following change of variables: we can write the potential term with a dependency on one single variable. Moreover, using Fujikawa's method [59], the change of variables in path integrals comes with the inverse Jacobian of the transformation matrix: Therefore, omitting the time dependency, This allow us to find the equations of Euler-Lagrange that constrain u 1 to be the same as the free particle (eq. (D.6)): As for the evolution of the norm, since one of the path integrals comes from a complex conjugate propagator: where ∂ (f ) and ∂ (i) are the derivatives with respect to all final and initial coordinates (formed by the 2-vectors (u 1f , v 1f ) and (u 1i , v 1i ) ). Putting all results together, where ∆l 1 = l 1f − l 1i , and ∆v 1 = v 1f − v 1i .