Anisotropic momentum broadening in the 2+1D Glasma: analytic weak field approximation and lattice simulations

In heavy ion collisions, transverse momentum broadening quantifies the modification of a hard probe due to interactions with the quark-gluon plasma (QGP). We calculate momentum broadening in the Glasma, which is the highly non-isotropic precursor stage of the QGP right after the collision. We show that the Glasma leads to anisotropic momentum broadening: high energy partons accumulate more momentum along the beam axis than transverse to it. The physical origin of anisotropic broadening can be traced back to differences in the shapes of chromo-electric and chromo-magnetic flux tubes in the Glasma. We provide semi-analytic results for momentum broadening in the dilute Glasma and numerical results from real-time lattice simulations of the non-perturbative dense Glasma.


I. INTRODUCTION
Understanding the properties of QCD matter in heavyion collisions provides many challenges and requires a combination of multiple approaches and tools for the different stages involved. An important probe that can be sensitive to all stages involved is high-energetic jets. They originate from initial hard scatterings and probe the emerging and evolving medium throughout its evolution before they hit the detectors in sprays of particles. The reduction of the energy through medium effects is known as jet quenching [1,2]. Experimental signatures include suppression of hadron spectra at high transverse momentum (p T ), modification of back-to-back correlations, and modification of reconstructed jets in nucleus-nucleus collisions as compared to proton-proton collisions. These effects have been extensively studied in the quark-gluon plasma stage using combined models, where the perturbative QCD effects stem at leading order from collisional energy loss and medium-induced gluon radiation [3][4][5][6][7]. Recently, it has been shown that the effect of jet quenching can be very sensitive to prehydrodynamic time-scales if the nuclear modification factor R AA and high-p T elliptic flow v 2 have to be matched simultaneously [8]. Therefore, a complete treatment of jet quenching should also consider modifications from the later hadron gas stage [9] as well as from earlier stages before hydrodynamics becomes applicable.
The Glasma [10,11] is a pre-hydrodynamic stage based on the color glass condensate (CGC) effective theory. It differs from the later quark-gluon plasma phase by nonequilibrated, highly anisotropic flux tubes with longitudinal chromo-magnetic and chromo-electric fields.
A first qualitative estimate of jet quenching from the Glasma has been given by considering synchrotron-like gluon emission in the coherent color fields of the Glasma [12]. Such contributions turned out to be smaller than the radiative energy loss in the QGP phase, but this result relies on a small angle approximation which strictly only holds for high-p T hadrons. In an energy regime ω 5 GeV, which is still important for the total energy loss, the contributions from the Glasma phase may no longer be negligible. Momentum broadening and energy loss due to collisions have been solved on a lattice using the Wong-Yang-Mills equations with hard collisions among the partons [13]. The estimate provided for a jet traversing a classical Yang-Mills field indicates a significant fractional energy loss of 10% -20% over the first τ 1 fm/c. One should take into account though that [13] did not use typical initial conditions for the Glasma and did not account for the longitudinal expansion of the system. A proper time expansion has been taken into account in a recent calculation for the leading analytic contribution for heavy quarks traversing the Glasma [14]. However, the truncation used is valid only for very small times and gauge invariance is not apparent in the chosen set-up.
In this work, we study the contribution of momentum broadening of a parton in the evolving Glasma in a systematic and gauge-invariant way. We compare analytic calculations with real-time lattice simulations and verify that they agree in the dilute limit. We observe anisotropic momentum broadening with larger momentum growth along the beam axis than transverse to it. This effect is also observed in our lattice simulations for the dense case and can be traced back to initial differences in the chromo-electric and chromo-magnetic field correlators of the Glasma.
The source of the anisotropic transverse momentum broadening that we find in the Glasma differs from various origins and mechanisms that have been proposed for the quark-gluon plasma: in [15] a calculation based on kinetic theory shows that the anisotropic momentum distribution in a longitudinally expanding quark-gluon plasma leads to larger jet broadening along the beam axis. Instabilities in non-Abelian plasmas [16] and turbulent color fields [17] can also give rise to anisotropic broadening. Similar effects have been found in anisotropic supersym-metric Yang-Mills plasmas based on the gauge/gravity duality [18][19][20]. Furthermore, the flow of the medium can deform the energy and momentum distribution of a jet [21]. In contrast, the physical origin of anisotropic momentum broadening investigated in this paper is related to particular properties of Glasma flux tubes. The model for the Glasma we consider is homogeneous in the transverse plane. Hence, no net flow is possible which allows us to rule out effects due to the flow of the medium. Moreover, we show that the anisotropy is unrelated to the longitudinal expansion of the Glasma. The anisotropy is directly linked to differences in the average spatial shapes of chromo-electric and -magnetic Glasma flux tubes (more precisely, their respective spatial autocorrelation functions) which exert different Lorentz forces on moving color charges.
This paper is organized as follows. In Section II we lay out the general formalism to calculate momentum broadening for a test particle traversing through a non-Abelian background field. The formulae that we develop are manifestly gauge invariant. In Section III we apply this formalism to the Glasma background field in a weak-field approximation and in the lattice formalism. In Section IV we finally present numerical checks of weakfield results with lattice simulations in the dilute limit and provide results for the dense case.

II. BASIC FORMALISM
In this section we present the basic formalism which we use to compute transverse momentum broadening of highly relativistic test particles as they traverse the Glasma. The methods we use are largely based on the definition of momentum broadening via a lightlike Wilson loop [22][23][24] and, equivalently, the dynamics of classical colored particles [25][26][27]. We review these approaches and derive a gauge invariant expression for transverse momentum broadening of an ultrarelativistic colored particle.
The geometry of the problem is shown in fig. 1. Two high energy nuclei move along the beam axis z and collide at z = 0 at time t = 0. During the collision, hard partons from the colliding nuclei can scatter and produce high energy jets. In fig 1, a single particle from such a scattering (either a quark or a gluon) is shown, which moves along the x-axis. We assume that this particle has a high initial momentum p x and zero momentum in the two orthogonal directions, p y = p z = 0. Around the same time, the Glasma is formed from interactions among the soft partons of the nuclei. The Glasma is a purely gluonic, highly occupied state and can be described as a classical color field. Initially, the Glasma consists of chromo-electric and chromo-magnetic flux tubes aligned with the beam axis z. In fig. 1 these flux tubes are depicted as colorful cylinders. As the system evolves and the particle traverses the Glasma, it interacts with the color fields of the Glasma and accumulates momentum Figure 1. Momentum broadening of fast partons from interaction with the Glasma. In this schematic we show a high energy particle created in the collision of two nuclei escaping towards the detector along the x axis with initial momentum px. The collision of the two nuclei produce the Glasma: a purely gluonic, highly occupied anisotropic precursor state to the quark-gluon plasma. As the particle moves through and interacts with the Glasma, it accumulates transverse momentum (in the sense of orthogonal to x, i.e. py and pz). This leads to momentum broadening ∆pz along the beam axis z and broadening ∆py within the transverse plane (spanned by x and y). along its path due to the non-Abelian Lorentz force exerted on the particle. For sufficiently high initial momentum p x , the particles' trajectory is almost unaffected by this interaction, but over time the particle accumulates transverse momentum p y and p z . This leads to a momentum broadening effect depicted in fig. 1 as a red cone around the x axis.
A note of caution is necessary regarding the terms longitudinal and transverse. In the context of jet momentum broadening, the terms transverse and longitudinal are defined relative to the jet trajectory. Using this nomenclature the x axis in the situation shown in fig. 1 would be the longitudinal axis, while y and z are transverse coordinates. This should not be confused with the longitudinal and transverse directions relative to the beam axis as used in the context of Glasma calculations. In the Glasma the longitudinal axis is identified with the beam axis z, while x and y span the transverse plane. For example, transverse momenta of gluons in the Glasma refer to momenta within the (x, y) plane. Transverse momentum broadening of a jet moving along x refers to broadening of momenta in the (y, z) plane.

A. Classical approach to momentum broadening
We start by sketching the basic derivation of momentum broadening in a more general context before focusing on non-Abelian background fields and eventually the Glasma. Consider an ultrarelativistic test particle traveling through a medium, for example a high momentum quark through the Glasma. It follows a straight, lightlike trajectory parametrized by where t is the time coordinate in the rest frame of the medium, the tangent vector u µ is constant and lightlike, i.e. u 2 = u µ u µ = 0, and x µ 0 defines the starting point of the particle at t = 0.
Without loss of generality we choose u µ = (1, 1, 0, 0) µ and x µ 0 = 0. In this context, x 1 is referred to as the longitudinal coordinate, while x 2 and x 3 are referred to as transverse coordinates x = (x 2 , x 3 ). The longitudinal part of the trajectory of the particle reads x 1 (t) = t, where we have set the speed of light c to one, i.e. c = 1. While the particle is moving through the medium, forces acting on the particle would deflect it from its straight trajectory. Assuming that the particle is highly energetic, the force has almost no effect on the trajectory if considered over short periods of time. However, we can still obtain an estimate on the momentum that the particle accumulates along its path.
The classical equation of motion of the particle in the frame of the medium is formally solved by where p µ = (p 0 , p , p ⊥ ) is the particle's four-momentum and F µ (x(t)) is an external force due to interactions with the medium. Here p = p x is the longitudinal momentum and p ⊥ = (p y , p z ) is the transverse momentum of the particle with respect to its trajectory. Assuming that the particle starts without transverse momentum, i.e. p ⊥ (0) = 0, we find where i ∈ {2, 3} is a transverse coordinate index. The force F i exerted on the particle, when computed from the Glasma, can be considered a random variable with zero mean but non-zero two-point function Here, . . . refers to the average over all realizations of the medium (e.g. the Glasma in our case). We stress that eqs. (5) and (6) are consequences of the Glasma picture (see section III) and should not be seen as additional assumptions. As a result, the averages of the transverse components p i (t) vanish, but their variances p 2 i (t) are given by where no sum over the index i ∈ {2, 3} is implied. The total transverse momentum is then given by The (instantaneous) jet broadening parameterq(t) is defined as the time derivative of the accumulated transverse momentum (see e.g. [26]) More generally, one can define a jet broadening parameter with respect to the transverse axis x i by taking the time derivative of a single momentum component Analogous to eq. (8) it holds that As we will see in section IV, it is reasonable to distinguish between the momentum broadening in the two transverse directions due to the special field configuration of the Glasma present immediately after the collision.

B. Momentum broadening in an Abelian background field
Before tackling the problem of non-Abelian fields, we derive p 2 i (t) for an electrically charged particle in an Abelian background field A µ . Using the same tangent vector as before, u µ = (1, 1, 0, 0) µ , the Abelian Lorentz force acting on a particle with charge q is given by where F µν = ∂ µ A ν − ∂ ν A µ is the Abelian field strength tensor associated with the background field A µ . Introducing the light cone coordinates x + and x − relative to the particle trajectory via we can write Inserting this into eq. (7) yields where, as before, no sum over the transverse index i is implied.
The result, eq. (15), obtained using the classical equations of motion for a charged particle, can also be found using the dipole approximation [22] via the expansion of a Wilson loop with lightlike edges aligned with the particle trajectory. A sketch of this Wilson loop is shown in fig. 2. Consider the rectangular Wilson loop W y+ spanned across one of the two transverse coordinates y (from y = 0 to y = L) and the lightlike axis x + (from x + = 0 to x + = L + ) given by where the transverse Wilson lines are and the lightlike Wilson lines are given by In the dipole approximation this Wilson loop can be related to the squared momentum p 2 y via for small transverse extensions of the loop L L + . To show this equivalence, we use Stokes' theorem to rewrite the integral around the closed, rectangular path as an area integral over the field strength tensor Taking the real part of W y+ and expanding for small transverse distances L L + yields Performing the medium average and rewriting the integrals along x + as an integration along the lightlike trajectory of the particle x(t), we find On the other hand, we can expand eq. (21) up to second order in L: Comparing coefficients yields This derivation can be repeated for the Wilson loop W i+ , where the transverse extent is aligned with the x i axis with i ∈ {y, z}. The general result reads where no sum over i is implied. We therefore recovered eq. (15), which we previously derived using classical particle dynamics.
More generally, if we align the transverse extension of the Wilson loop along some other transverse direction, we obtain the expression for L L + , where L denotes the length of the loop along x i . In order to directly compute p 2 i we expand the Wilson loop W i+ as a function of L and take the expectation value of the quadratic term: C. Momentum broadening in a non-Abelian background field

Derivation from the Wilson loop
The results of the previous sections can be generalized to non-Abelian fields and particles with color charge. In Yang-Mills theory with gauge group SU(N c ), Wilson loops are elements of SU(N c ) and the trace over a Wilson loop depends on the chosen representation R of SU(N c ). A straightforward generalization of eq. (28) is given by where D R is the dimension of the representation and Tr [. . . ] R denotes the trace operation in the given representation. The choice of representation of SU(N c ) decides the species of particle that we consider. For quarks we choose the fundamental representation F and for gluons we choose the adjoint representation A. We add a subscript R in p 2 i R to denote that the accumulated squared momentum depends on the representation. As before, p 2 i R is given by the second coefficient of the expansion in L: where In fact, it is not necessary to actually compute the quadratic coefficient in order to obtain p 2 i (t) . Enforcing the unitary constraint for the expansion eq. (33), yields two constraints for the Taylor coefficients Using the second relation we can write The accumulated momenta p 2 i (t) are therefore given by the linear coefficient of the Wilson loop, which drastically simplifies the calculation.
A detailed derivation of the linear coefficient W i+ can be found in appendix A. The final result for the squared accumulated momentum (see eq. (A24)) reads where we have defined the parallel transported field strength tensor and the lightlike Wilson line W + (b, a) connecting x + = a to x + = b Here, P denotes path ordering along x + . In the above expression, the gauge field A + is evaluated along the lightlike particle trajectory.
Our main result eq. (38) is the non-Abelian generalization of eq. (27). It is evident that this expression for p 2 i (t) R is gauge invariant: the Wilson lines in eq. (39) ensure that under a general gauge transformation Ω(x), the parallel transported field strength tensor transforms according to The trace over the square of the field strength tensor in eq. (38) is then guaranteed to be gauge invariant.

Derivation using classical colored particle dynamics
Having derived eq. (38) using the Wilson loop, we now show that the same result can be found using classical colored particle dynamics. The equivalent derivation is based on colored particles affected by a non-Abelian Lorentz force as given by Wong's equations where the non-Abelian force is given by Here, Q a (t) are the components of the time-dependent color charge of the particle and f abc are the antisymmetric structure constants of SU(N c ). In order to rewrite everything in terms of group and algebra elements of SU(N c ), we use where T a are the generators of SU(N c ), and the relation where R ∈ {F, A} and T F = 1/2 and T A = N c in the fundamental and adjoint representations respectively. We can write the equations of motion as The formal solution to the second equation is given by where If we choose a lightlike trajectory, i.e. x µ (t) = u µ t with u µ = (1, 1, 0, 0) µ , we recover the lightlike Wilson line where x + (t) = √ 2t. For readability, we suppress other coordinates on which A µ depends. The equations of motion for the transverse momenta read Repeating the same steps as before, we find that the accumulated squared transverse momentum reads where Q a (0) are the components of the initial color charge of the particle and no sum over i is implied. The next step, following [25,26,28], is to perform an average over Q a (0) using and (as done in [26]) divide by the dimension of the representation D R . This yields Using the relation where X and Y are arbitrary elements of the Lie algebra of SU(N c ), we can further simplify eq. (56). The result reads If we define the parallel transported field strength tensor as in eq. (39), we recover the same result as eq. (38):

D. Casimir scaling of accumulated momenta
Having derived an expression for p 2 i R given a representation R of SU(N c ), it is easy to see that the accumulated transverse momentum scales with the Casimir C R of the representation. Using eq. (37) and the fact that W (1) i+ is anti-hermitian and traceless, we can write and consequently If we look at the ratio of the accumulated momentum of a gluon and a quark, we find the well-known result For N c = 3, this factor evaluates to 9/4.

III. MOMENTUM BROADENING IN THE GLASMA
In this section we apply eq. (59) to the boost-invariant Glasma. As noted in the beginning of section II, care has to be taken with nomenclature: in the Glasma, transverse coordinates are usually x = x 1 and y = x 2 , whereas z = x 3 is the longitudinal coordinate along the beam axis. In contrast, the transverse momentum p 2 ⊥ refers to the momentum perpendicular to the particle trajectory, which we take to be aligned with one of the transverse directions in the Glasma. Here, we choose the particle to move along the x axis. Additionally, we use a different definition for light cone coordinates. They are given by where, compared to eq. (13), we use x 3 in place of x 1 .
Note that x 3 = z is the beam axis. We start with a quick review of the boost-invariant color glass condensate framework [10,[29][30][31] and the Glasma. In the CGC, the hard partons of a nucleus moving at the speed of light along the positive z axis (referred to as nucleus "A") are described by a color charge density ρ A , which is highly localized around x − = 0. The color current of such a nucleus is given by (under the assumption of covariant gauge, i.e. ∂ µ A µ = 0) where x = (x, y) is the transverse coordinate vector.
Analogously, the color current of nucleus "B", which moves in the opposite direction, is given by Associated with these currents J µ A and J µ B are the gauge fields A µ A and A µ B , which represent the soft gluons of the nuclei. In the appropriate light cone gauges A + = 0 for nucleus "A" and A − = 0 for "B", the gauge fields become purely transverse and are given by The lightlike Wilson lines along x ± are given by where ∆ T is the two-dimensional Laplace operator in the transverse plane and m is an infrared regulator. The infrared regulator is used to model the phenomenon of color confinement in the nucleus. It is common practice to choose the value of m on the order of Λ QCD ≈ 200 GeV in order to enforce color neutrality on transverse length scales of ∼ 1 fm. In the ultrarelativistic limit, where nuclei are assumed to be infinitesimally thin along x ∓ , the transverse gauge fields simplify to where θ(x) denotes the Heaviside step function and the asymptotic Wilson lines are given by To specify the charge densities ρ (A,B) , a model has to be chosen. In this work, we focus on the McLerran-Venugopalan (MV) model [32,33]. For a nucleus moving along x + (nucleus "A"), the random charge density ρ A is assumed to be Gaussian and given by [34] where g is the Yang-Mills coupling constant and µ A is the MV model parameter given in units of energy. The function λ A (x − ) defines the shape of the nucleus along x − and is highly peaked around The expectation value . . . refers to the functional integration over all realizations of ρ, as specified by the probability functional W MV [ρ], i.e.
where the functional W MV [ρ] is given by with and the normalization constant Z is determined by We note that the MV model is both homogeneous and isotropic in the transverse plane as evident from eq. (72). Consequently, nuclei described by the MV model have no finite radius.
For a nucleus moving along x − (nucleus "B"), the charge density correlator eq. (72) is defined analogously with x + in place of x − .
The classical Glasma field, which is created in the collision of the two nuclei, is obtained from solving the classical Yang-Mills equations in the future light cone, which is spanned by the Milne coordinates: proper time τ = √ t 2 − z 2 ≥ 0 and space-time rapidity η = 1/2 ln ((t + z)/(t − z)). In general, there are no analytical solutions for the gauge field in the future light cone for arbitrary τ > 0. It is however possible to deduce that gauge-invariant observables in the future light cone need to be invariant under boosts along the beam axis. This boost symmetry can be extended to gauge fields (without loss of generality) if rapidity dependent gauge transformations are disregarded. The fields in the future light cone therefore only depend on τ and x, i.e. A µ (x) = A µ (τ, x). Furthermore, using the assumption of boost invariance, it is possible to solve the Yang-Mills equations for the fields at the boundary of the future light cone at τ = 0. The solution in temporal gauge A τ = 0, is given by [35] A and Taking the above expressions for A i and A η at τ = 0 as initial data, we need to solve the temporal gauge Yang-Mills equations in Milne coordinates for τ > 0: Here, the first two equations define the canonical momenta P i and P η . The initial conditions in these variables are given by Given a solution to the Yang-Mills equations, we can evaluate the accumulated transverse momentum. Note that x + in eq. (59) refer to x + = (τ + x)/ √ 2 as defined in eq. (13). The accumulated momenta in y and z are given by where f y and f z are defined as Each of the terms in parenthesis are understood to be evaluated along the lightlike trajectory x µ (τ ) = (τ, τ, 0, 0) µ . In temporal gauge, the Wilson lines reduce to where A x (τ ) is evaluated at x µ (τ ). In terms of chromoelectric and -magnetic field strengths, the functions f y and f z are given by The expectation values in eqs. (89) and (90) refer to the medium average introduced in sec. II. In the case of the Glasma, this expectation value corresponds to an average over the charge densities of both nuclei, i.e.
Due to homogeneity and isotropy of the MV model, the final results for p 2 y (τ ) and p 2 z (τ ) do not depend on the starting position of the particle. This would not be the case for more complicated initial conditions, which, for instance, account for the finite radii of nuclei. More generally, if µ 2 A in the charge density correlator eq. (72) is promoted to a function µ 2 A (x) of the transverse coordinates x, the results for p 2 y (τ ) and p 2 z (τ ) will depend on the initial transverse coordinate of the particle.
It should also be noted that the MV model initial conditions naturally fulfill eqs. (5) and (6), since no single color component is preferred over any other. Therefore, no additional assumptions are necessary in order for eqs. (89) and (90) to hold.
There are two main problems with evaluating eqs. (89) and (90): firstly, it is difficult to obtain fully nonperturbative solutions to the Yang-Mills equations. Secondly, even if a solution A µ is known, eqs. (89) and (90) are highly non-linear functionals of the gauge field due to the presence of the lightlike Wilson lines. In order to make progress, some approximations are necessary.
In the rest of this section we discuss two methods we can use to approximately compute p 2 ⊥ . The first method is to assume that the fields of the nuclei are weak and therefore perturbative analytic solutions for the Glasma field can be obtained from a linearization of the Yang-Mills equations. The solution to the linearized equations can be plugged into eqs. (89) and (90), where contributions from the Wilson lines are negligible at leading order. The second method is to solve the non-perturbative Yang-Mills equations numerically on a lattice for the case of strong fields. This is done using real-time lattice gauge theory. Our expression for p 2 ⊥ can be evaluated numerically using the solution obtained from the lattice calculation. Although the weak field approximation is not as phenomenologically interesting as the fully non-perturbative strong field case, it provides a non-trivial check of the numerics of the lattice calculation. Additionally, the weak field approximation leads to a transparent interpretation of momentum broadening from the early stages of heavy ion collisions in terms of the properties of Glasma flux tubes.

A. Weak field approximation
In this section we first summarize the weak field approximation, which has been discussed more carefully in [36]. Our main goal is to use the results of [36] to find semi-analytical expressions for eqs. (89) and (90). A more detailed version of the calculation can be found in [37].

Weak field limit of the Glasma
The weak field approximation is based on the assumption that the charge density of the nucleus is small and that the Coulomb field of the nucleus is therefore small as well. This allows us to treat the Glasma initial conditions perturbatively and the further time evolution of the Glasma, once it has been created, as effectively Abelian. We start by deriving the weak field approximation to the Glasma initial conditions. Using the fact that ρ is small, we expand the Wilson lines in eq. (70) in powers of gρ: where the potentials φ are related to the charge densities of the incoming nuclei via In the second line P x ∓ denotes path ordering along x ∓ . The transverse gauge fields of the nuclei up to quadratic order in ρ are given by Inserting this expression into eqs. (79) and (80), we can derive the weak field Glasma initial conditions: In terms of field strengths, the Glasma at τ = 0 corresponds to purely longitudinal chromo-electric and -magnetic fields [10,38], E z = P η and B z = −F xy . In the weak field limit, the initial values for P η and B z read Note that at leading order in the weak field expansion, the fields E z and B z do not contain quadratic φ (2) terms, which implies that path ordering effects do not contribute at this order of the perturbative series. Provided that the initial field strengths are to be considered small, the further time evolution of the Glasma is approximately Abelian, see e.g. [36] and [39]. A general solution to the linearized (Abelian) Yang-Mills equations (83) -(86), which is consistent with the Glasma initial conditions P η = 0, A i = 0, P i = A η = 0 at τ = 0, is given in terms of the Green's functions: whereẼ z (0, k) andB z (0, k) are the Fourier components of the initial longitudinal field strengths and J α (x) are Bessel functions of the first kind. We also use the shorthand k (. . . ) for momentum integrals: The time evolution of the transverse field strengths B y = D x A η /τ and E y = P y /τ in eqs. (94) and (95) is given by

Momentum broadening in the dilute Glasma
Having summarized the weak field approximation of the Glasma, we now apply it to compute the accumulated momenta of fast partons as derived in section II. In the Abelian approximation, the Wilson lines can be dropped from eqs. (94) and (95). Inserting the above relations for the time dependent field strengths and evaluating them on the trajectory, we find where we define The accumulated momenta can be written as Note that the only difference in p 2 y and p 2 z is due to the different initial electric and magnetic autocorrelation functions. The time evolution of the fields and the movement of the colored test particle through the Glasma is hidden in γ(τ, k) and is the same for both p 2 y and p 2 z . Since the correlators of E z and B z are explicitly evaluated at τ = 0, the integrals over proper time only act on the functions γ(τ, k). Using eqs. (103) and (104) and performing the averages using the MV model, we find that the initial correlators in momentum space are given by Here, we use p × k to denote the two-dimensional cross product in the transverse plane. Using the fact that γ(τ, −k) = γ * (τ, k), we find Additionally, since m is the only dimensionful scale that occurs in the integrals, we can re-parameterize both expressions and find with the dimensionless function with dimensionless proper timeτ = mτ . Using the above expressions we can determine the accumulated transverse momenta as a function ofτ by (numerically) solving the integrals in d (E,B) (τ ). Numerical results are presented in section IV.

B. Lattice approximation for strong fields
In this section we discuss our approach to compute the accumulated transverse momenta of fast partons in the Glasma using real-time lattice gauge theory methods. This allows us to numerically determine p 2 (y,z) (τ ) for non-perturbatively large Glasma fields, which are produced in relativistic heavy-ion collisions.

Real-time lattice gauge theory description of the Glasma
The real-time lattice gauge theory approach is a finite difference method to numerically solve the Yang-Mills equations and is usually formulated the following way: we discretize the transverse plane on which the gauge field A µ is defined as a regular, square lattice of size N T × N T with transverse lattice spacing a T . The transverse length of the lattice is given by L T = N T a T . The discrete time step is denoted as ∆τ and the discrete proper times are given by τ n = n∆τ . For convenience, we choose ∆τ = a T /n τ , where n τ ≥ 2 is an even integer. Requiring that n τ is an even integer simplifies the lattice formulation of transverse momentum broadening in terms of test particles (see section III B 2). The lower bound of n τ = 2 guarantees numerical stability of the finite difference scheme. As we use the MV model as our model for nuclei, we can employ periodic boundary conditions for the transverse lattice.
Instead of the usual pair of fields and conjugate momenta, (A µ , P µ ), lattice gauge theory is formulated in terms of gauge links or link variables, which are the Wilson lines (usually in the fundamental representation) connecting the lattice sites of the regular lattice. The transverse gauge field components A i (τ, x) are replaced by the gauge links U x,î (τ n ), which connect the lattice site x and the neighboring site x +î withî denoting the unit vector along the coordinate axis x i . For small lattice spacing a T the gauge links are related to the usual gauge fields via Gauge links starting at x and going in the opposite direction ofî are denoted by The conjugate momenta P i (τ, x) are replaced by P i x (τ ) and are defined in the continuous time limit ∆τ → 0 as As with the gauge links U x,î , P i x are defined at mid-point of the edge connecting x and x +î. The rapidity component of the gauge field A η (τ, x) is replaced by A x,η (τ ), which is defined at the lattice sites. The conjugate momentum P η x (τ ) is defined via Using this set of variables, the Yang-Mills equations (83) -(86) are replaced by a leapfrog scheme for finite time steps ∆τ > 0 (see e.g. [40,41]): where we define the anti-hermitian, traceless part of a matrix X in the fundamental representation as (131) The so-called plaquette variables U x,îĵ are 1 × 1 Wilson loops on the lattice and given by Note that the conjugate momenta are defined at fractional time steps τ n+1/2 , which allows the leapfrog scheme to be accurate up to quadratic order in ∆τ . The discrete analogues of the gauge covariant derivative D µ are given by the forward and backward differences: where all terms on the left are evaluated at proper time τ n . The second derivative D 2 i is defined as where no sum over i is implied. The Glasma initial conditions have to be discretized on the lattice as well. The details of the discretized initial conditions have been worked out in detail in [40]. The starting point is to approximate the charge density correlator of the MV model (see eq. (72)) on a threedimensional lattice. For simplicity we focus on a nucleus moving along x + without loss of generality. We assume a rectangular shape λ(x − ) along x − and split the nucleus into N s "color sheets" along the longitudinal direction, which allows us to account for path ordering along x − [42]. The discretized correlator then reads ρ a x,m ρ b y,n = where the indices m and n refer to the different color sheets. The values ρ a x,m of the lattice color charge density are uncorrelated, normally distributed random numbers. After having generated a particular random configuration of color charges, we solve the discretized Poisson equation in the transverse plane for each individual color sheet: The above equation can be readily solved in terms of Fourier coefficients. The solution reads where the lattice momentumk 2 T is given bỹ and m is the infrared regulator. Finally, the lightlike Wilson line can be approximated as a path ordered product of the individual sheets [42]: where T a F are the generators in the fundamental representation of SU(N c ). The transverse gauge links of the nucleus are given by The above outlined procedure is performed for each nucleus separately. The discrete analogues of eqs. (79) and (80), which provide the initial Glasma fields right after the collision at τ = 0, are given by [40] U and Note that eq. (141) has to be solved numerically, except for N c = 2, where an analytic solution is known [40].

Lattice approximation for transverse momentum broadening of fast colored particles
Having defined the initial conditions and the discrete field equations, the next step is to find discretizations of eqs. (94) and (95). In order to obtain a consistent procedure we take inspiration from the nearest grid point (NGP) scheme of the colored particle-in-cell (CPIC) method (see e.g. [43][44][45]). We imagine a particle with color charge Q moving across the two-dimensional transverse lattice with trajectory x(τ n ). As the particles' position changes with time, the color charge Q of the particle contributes to the lattice charge density ρ x (τ n ) at the lattice site x which is nearest to x(τ n ). Each time the nearest grid point of x(τ n ) changes from x to e.g.
x +î, the color charge Q must be color rotated using the Wilson line connecting x and x +î, i.e. the gauge link U x,î (τ n ). In this example, the color rotation of the charge reads Color rotation of the charges is necessary in order to guarantee local gauge-covariant charge conservation.
In this work, we only consider fixed, lightlike trajectories along the x-axis, i.e. x µ (τ n ) = (τ n , τ n , 0, 0) µ + x µ 0 , where x 0 denotes a starting lattice site at τ 0 = 0 in the transverse plane. At time steps t n = nn τ ∆τ = na T the particle position ends up at one of the sites of the transverse lattice, denoted by x n . The nearest grid point changes at timest n = (n + 1/2)a T at which point the charge undergoes color rotation with the gauge link connecting the appropriate edge on the lattice. The consecutive color rotations of the particles' charge amount to the lightlike Wilson lines in eqs. (94) and (95). We therefore use the approximation The field strengths E y , E z , B y and B z in eqs. (94) and (95) acting on the particle are most easily computed at the lattice sites x n where the particle ends up at times t n . The evaluation of the non-Abelian Lorentz force and the Wilson line are therefore not done at the same time. Rather, we alternate between computing field strengths at t n and the Wilson line att n . We use the following approximations for the field strengths, which are accurate up to quadratic order in both time step ∆τ and lattice spacing a T : + U xn,−ŷ (t n ) P y xn−ŷ (t n + ∆τ 2 ) Finally, we put all the above expressions together to obtain the lattice approximation of eqs. (94) and (95): with the field strengths evaluated at t n and the Wilson line evaluated att n .
In order to compute the accumulated transverse momenta, we have to integrate over τ , which we approximate using a simple sum The accumulated momenta on the lattice then read The expectation value in the above expression is a simple average over the random color charge densities used in the initial conditions. The statistical quality of the average can be improved by additionally averaging over the starting positions x 0 of the particle. This is possible because the MV model used in this work is homogeneous in the transverse plane, which implies that the results must not depend on x 0 if one averages over enough configurations. The additional average over x 0 simply accelerates the convergence of the expectation value, thus reducing the computational cost of numerically determining eq. (154).
In practice, we actually evaluate eq. (154) only for quarks, since the variables in our code are all in the fundamental representation. However, as shown in section II D, a simple rescaling of our results yields the accumulated momenta for gluons.  Figure 3. Accumulated transverse momenta in the dilute Glasma: weak field approximation vs. lattice approximation. We performed simulations using both SU(2) (crosses and circles) and SU(3) (triangles and squares) to show that our lattice results scale with DA = N 2 c − 1. The results for the weak field approximation are shown as red and blue lines. It is evident that our lattice simulations agree with the weak field approximation. While p 2 z (mτ ) increases monotonically with mτ , p 2 y (mτ ) converges to a finite value after mτ > 1. As a result, the two momentum components exhibit a large anisotropy which grows with mτ .

IV. RESULTS AND DISCUSSION
In the first part of this section we study the dilute Glasma, where we are able to perform both semi-analytic calculations and pure lattice calculations. We find that fast partons in the dilute Glasma exhibit large anisotropic momentum broadening, which can be traced back to infrared behaviour of the initial chromo-electric and -magnetic fields in the Glasma. In the second part we study the dense Glasma by applying the lattice simulations to more realistic initial conditions. This allows us to estimate the effect of the Glasma on quark and gluon jets. We find that also in the case of the dense Glasma the effect of anisotropic momentum broadening persists.
The methods to compute the accumulated transverse momenta presented in the previous sections have been implemented as part of curraun [46], an open-source code for simulating the Glasma in 2+1D using GPUs written in Python and Numba. Jupyter notebooks are provided with the code which can be used to replicate the results presented in this section.

A. The dilute Glasma
We compute the accumulated transverse momentum in the weak field approximation eq. (121) by numerically solving the integrals in eq. (122). We focus on collisions of identical nuclei and set Q A = Q B = Q. The dimensionless transverse momenta for quarks are then given by We compute these integrals numerically using SciPy [47]. The same quantities can be computed using our lattice simulation. We use a transverse lattice with N T = 1024 cells per side and a time step of ∆τ = a T /n τ with n τ = 16. In order to approach the weak field limit we need to choose our infrared regulator m to be much larger than Q = g 2 µ. For the results presented here we used m = 100 Q. On the lattice we have to make sure that there is no interference from the cutoffs due to the lattice approximation, i.e. we make sure that L −1 is the infrared cutoff due to the finite size of the lattice and a −1 T is the ultraviolet cutoff due to the finite size of each lattice cell. In our simulations we chose L T = N T a T = 20m −1 to satisfy this condition. For the discretization of the MV model we use N s = 50 color sheets per nucleus. The results of eq. (154) have been averaged over 20 collision events. In order to show that p 2 (y,z) /(N 2 c − 1) is indeed independent of N c , we performed our simulations with both SU(2) and SU(3). The results of both our lattice calculation and the weak field approximation are shown in fig. 3, which show good agreement between our purely numerical lattice calculation and semi-analytical weak field approximation. We take this as evidence for a correct implementation of our lattice code.
The first important observation in fig. 3 is that although the accumulated momentum p 2 z along the beam axis keeps increasing with time, the y components starts settling after τ m −1 to some constant value. The result is a large anisotropy p 2 z / p 2 y which increases monotonically with time τ . Consequently, the widening momentum cone around the particle trajectory becomes highly elliptical with much more pronounced broadening in z or, equivalently, rapidity η. This anisotropy is rather interesting in the light of eq. (155). The only difference between p 2 y and p 2 z is the type of correlator that is integrated over: either the initial magnetic correlator c B (|k|/m, 1) or the electric correlator c E (|k|/m, 1) (see eqs. (117) and (118)). The term in eq. (155) which takes care of the time evolution of the system is hidden in G(mτ, k T /m) and is the same for both momentum components. The interpretation is that in the dilute Glasma, a high energy parton receives momentum kicks in the z direction only from the chromo-electric Glasma flux tubes, whereas the momentum kicks in the y direction are due to chromo-magnetic flux tubes. The physical reason for anisotropic momentum broadening is therefore not the longitudinal expansion of the Glasma as one might naively assume. The anisotropy emerges due  To be more precise: although the boost-invariant Glasma is anisotropic in the sense that it is a longitudinally expanding system and that at τ = 0 the field strengths all point in the direction of the beam, the anisotropy observed in the momentum broadening is of a different physical origin and is related to the initial correlations in the Glasma. In order to investigate the reason for this anisotropy more closely, we plot the autocorrelation functions of E z and B z at τ = 0 in both momentum and coordinate space. In the lattice approximation, we compute the gauge invariant field correlator, i.e.
at τ = 0, where U x→y (τ ) is the Wilson line connecting the lattice site x and y. We include the Wilson line to guarantee gauge invariance, even though this only yields negligible corrections of higher order in the weak field limit. The momentum space correlator on the lattice is obtained by Fourier transforming the spatial correlator. The magnetic correlator is computed analogously. The results are shown in fig. 4.
First focusing on the correlators in momentum space, we see that chromo-electric and -magnetic fields behave drastically differently in the infrared region |k| < m. While the electric correlator admits a positive, non-zero value, the magnetic correlator tends to zero: On the other hand, for large momenta |k| m both correlators agree: The momentum broadening anisotropy of fast partons seems to be strongly affected by the infrared modes of the correlators. In the case of other observables that are less sensitive to the behaviour in the infrared, the mismatch in the two correlators for |k| < m can be generally neglected. For example, the energy density contributions of E z and B z at τ = 0 are given by where Λ m is an ultraviolet cutoff. In the above expressions, the factor of k from the integration measure suppresses the infrared modes of c E (|k|, m) and c B (|k|, m). Additionally, the integrand only falls off slowly with increasing |k| and therefore the main contributions to the (ultraviolet divergent) energy density come from high momenta k, where c E (|k|, m) ≈ c B (|k|, m) and therefore ε E (τ = 0) ε B (τ = 0) for Λ m.
The fact that both electric and magnetic fields roughly yield the same initial energy density, combined with the Abelian time evolution of the Glasma, would suggest that there is, on average, no qualitative difference between a chromo-electric and a chromo-magnetic Glasma flux tube. However, this is only true at short distances. The different behaviour at larger length scales can be seen in the Fourier transformed correlators from the weak field approximation and our lattice simulation shown in fig. 4 (b). For small distances mr 1 both correlators and the large k approximation of eq. (159) agree. At intermediate distances mr ≈ 1 however, the behaviour of the electric and magnetic correlators differs: the electric field correlator shows positive correlation for all mr, but the magnetic correlator has characteristic anti-correlated regions around mr ≈ 1. Intuitively, this anti-correlation means that when looking at two points x 1 and x 2 in the transverse plane separated at a distance of |x 1 − x 2 | ≈ m −1 , it is very probable that the two magnetic field values B z (x 1 ) and B z (x 2 ) will have opposite signs. In contrast, the signs of E z (x 1 ) and E z (x 2 ) are more likely to be the same. The phenomenon of anti-correlated regions of the magnetic field at τ = 0 has also been observed in the dense Glasma [48].
In the dilute Glasma a physical interpretation of the zero mode of the correlators is given by the electric and magnetic flux in the transverse plane. We define the magnetic flux Φ B (R) at τ = 0 over a circular area of radius R with origin x 0 = 0 as Although the expectation value of Φ B vanishes trivially, the expectation value of the magnetic flux with respect to the field at the origin, i.e.
does not. In the limit of R → ∞ we can use the fact that and find The vanishing of the zero mode implies that the total magnetic flux in the transverse plane vanishes. In contrast, the electric flux is non-zero: In the context of transverse momentum broadening of fast partons, the anti-correlation of the magnetic field (or the vanishing of the magnetic flux) has immediate consequences on the momentum broadening anisotropy. Consider a fast colored particle immediately after the collision that traverses the Glasma, which consists of electric and magnetic flux tubes. A fast particle that gets a particular momentum kick in the y direction from a chromo-magnetic flux tube will, with high probability, get an anti-correlated momentum kick in the opposite direction after it has travelled a characteristic distance of m −1 , reducing the accumulated momentum p y of the particle. On the other hand, a particle that moves through an electric flux tube will not receive an anti-correlated kick in z and therefore the accumulated momentum p z will be (on average) higher than p y . This phenomenon can be seen clearly in fig. 3, where it is shown that p 2 y (τ ) indeed flattens after a characteristic time of m −1 , corresponding to distances of m −1 at the speed of light. In contrast, the other component p 2 z (τ ) does not exhibit such a decelerating effect as the electric field is always positively correlated.
In summary, we find that anisotropic momentum broadening of fast partons moving through the dilute Glasma is closely linked to qualitative differences in the typical spatial shapes of chromo-electric and -magnetic Glasma flux tubes. Up until now we have only studied the dilute Glasma, which is of less phenomenological interest than the fully non-perturbative dense Glasma. The dilute Glasma, which neglects all non-linear interactions of the Yang-Mills field, should only be seen as a toy model for the dense gluonic matter created in relativistic heavy ion collisions. The dilute case is unphysical in the sense that the effective classical description in the CGC formalism is only valid for large fields and weak coupling. Moreover, the governing momentum scale in the weak field approximation is the infrared regulator m, which has to be introduced by hand in order to obtain finite results. The fact that all results depend strongly on this artificial energy scale is troubling. One might expect that the effect of anisotropic momentum broadening might be nothing but an artifact of the infrared regulation. Fortunately, lattice gauge theory allows us to study the Glasma also in the case of strong fields, where the non-linear dynamics of the Yang-Mills fields seem to cure the infrared divergences that occur in the dilute limit and results should only weakly depend on the infrared regulator m. In particular, we will see that the governing momentum scale will be given by the saturation momentum Q s in the sense that our results for the accumulated momenta in the dense regime look qualitatively similar to fig. 3 if one replaces m with Q s .

B. The dense Glasma: relativistic heavy ion collisions
After having thoroughly studied momentum broadening in the dilute Glasma, we now focus on the dense Glasma, which is created in the collision of two relativistic heavy nuclei. As before, the color charge densities are generated according to the MV model. A dense Glasma is produced when the ratio of the infrared regulator m and the saturation momentum Q s (which is roughly g 2 µ) becomes small: m/(g 2 µ) 1. On the other hand, the dilute Glasma corresponds to m/(g 2 µ) 1. To access the regime of non-perturbatively strong initial fields we perform lattice simulations as described in section III B. In order to render our simulations more physical, we work at fixed saturation momentum Q s . This done by choosing g 2 µ for some fixed Q s according the numerical results found in [49]. A summary of the numerical parameters used in our simulations is shown in table I.
As in the case of the dilute Glasma, we work with square lattices of size N T = 1024 and a time step of ∆τ = a T /16, where a T is the transverse lattice spacing. For g 2 µL = 100, we achieve a lattice resolution of g 2 µa T ≈ 0.1. This is sufficient to resolve Glasma flux tubes which roughly have a diameter of Q −1 s ≈ (g 2 µ) −1 [10]. For the infrared regulator m we try out a few different values: m/(g 2 µ) ∈ {0, 0.05, 0.1, 0.2}. In the cases of m ≥ 0.05g 2 µ, the system size L is large enough to resolve multiple color neutral domains with diameter ∼ m −1 . In the case of m/(g 2 µ) = 0 we simply eliminate the zero mode of the charge density (see eq. (137)) to implement color neutrality on the size of the system L [50]. This cor- Table I. Parameters used in simulations of the dense Glasma. The numerical values for the ratios of m, Qs and g 2 µ were originally obtained by [49]. The provided values correspond to g 2 µL = 100 with Ns = 50 color sheets per nucleus. For purposes of illustration, on the right side of the table we show  example values for an intermediate saturation momentum  For gluons the momenta should be multiplied with C A /C F = 9/4. There are a few observations we make: momentum broadening at early times τ 5/Q s behaves remarkably similar to the dilute Glasma calculation, see fig. 3. The broadening within the transverse plane p 2 y rises sharply until τ ≈ Q −1 s , where it starts to flatten out. Broadening along the beam axis p 2 z is not affected by this and depends strongly on the infrared regulator m. We observe that the momentum broadening anisotropy starts to develop around τ ≈ Q −1 s and keeps growing until τ ≈ 10/Q s . The anisotropy shows strong dependence on the infrared regulator m and shrinks for smaller values of m. At later times τ 10/Q s the large growth of p 2 z stops and is reversed. Surprisingly, the momenta along the beam axis are not broadened anymore but start to "focus". Consequently, the anisotropy shrinks and in the case of m = 0 almost vanishes entirely. One should keep in mind that the surprising behaviour of shrinking or almost vanishing momentum broadening anisotropy happens at rather late times. However, the classical approach is likely not valid after τ ≈ 10/Q s . It is therefore doubtful if this particular effect is of phenomenological interest. On the other hand, the early-time momentum broadening of fig. 5 (a) is computed in the regime where the gluon occupation number of the Glasma is still large and a classical description can be applied. Therefore, we expect a similar behaviour of p 2 y and p 2 z also for more realistic CGC models of nuclei such as IP-Glasma, which, even though it accounts for non-homogeneity of the color charge density of nuclei due to varying nucleon positions, is at its core based on the MV model used in this work. In fig. 5 we also plot the total accumulated momentum p 2 ⊥ = p 2 y + p 2 z as a function of Q s τ . We observe that after τ ≈ 10/Q s a quark will have picked up roughly p 2 ⊥ ≈ Q 2 s . This value seems appropriate as the typical momentum of a gluon in the Glasma is roughly Q s [51].
Looking at the infrared dependence of the momentum  broadening anisotropy in fig. 5 (b) one might wonder if the large anisotropy around Q s τ ≈ 5 − 10 still persists if the infrared regulator m is reduced further. Recall that in the case of m/(g 2 µ) = 0 the infrared divergence of the MV model is regulated by the system size L −1 . It is therefore possible to investigate the infrared dependence of p 2 z / p 2 y for m = 0 by studying the limit of large systems where g 2 µL → ∞. These results are shown in fig. 6. We simulate the Glasma using the gauge group SU(2) for fixed lattice spacing g 2 µa T and varying lattice size from N T = 100 to N T = 3200 (corresponding to the range of g 2 µL = 100g 2 µa T to g 2 µL = 3200g 2 µa T ). The results in fig. 6 have been obtained by averaging over 200 initial conditions for each value of g 2 µa T and g 2 µL. Although our results reveal a dependence of the anisotropy on g 2 µL around the typical values used in the Glasma (g 2 µL = 100 as used for all other results in this work), the anisotropy never completely vanishes, i.e. we find p 2 z / p 2 y > 1.5 also for large values of g 2 µL. It is important to note that we evaluate p 2 z / p 2 y at fixed proper time τ 0 = 6/(g 2 µ), where the anisotropy is large according to fig. 5. At asymptotically large τ , the anisotropy almost vanishes. We conclude that while the accumulated momenta along the beam axis p 2 z are sensitive to the infrared regulator m (or L −1 in the case of m = 0), the momentum broadening anisotropy does not appear to be an artifact related to the way infrared modes are treated in the Glasma initial conditions. However, for m → 0 there is only a short time window in which the anisotropy can become large.
In order to further compare our non-perturbative results to the dilute Glasma and to better understand the time dependence of p 2 y and p 2 z , we plot the gaugeinvariant field strength autocorrelation functions at τ = 0 from our lattice simulation and their Hankel transformations in fig. 7. As in the dilute Glasma case, we observe regions of negative correlation around r ≈ Q −1 s in fig. 7 (b). The same behaviour has been observed in [48]. We presume that this behaviour is linked to the flattening of p 2 y around τ ≈ Q −1 s in fig. 5. It is interesting to note that in the dense Glasma the regions of anti-correlation are much more pronounced in comparison to the dilute limit. This is particularly apparent in the momentum space representation in fig. 7 (a), where the zero mode of g 2 µaT = 1/2 g 2 µaT = 1/4 g 2 µaT = 1/8 g 2 µaT = 1/16 Figure 6. Momentum broadening anisotropy p 2 z / p 2 y at τ0 = 6/(g 2 µ) as a function of g 2 µL for different transverse lattice spacings g 2 µaT from Glasma simulations using the gauge group SU (2). The dashed vertical line corresponds to the value g 2 µL = 100 used in the rest of this work. The anisotropy shows a dependence on g 2 µL for g 2 µL < 500, but this dependence becomes weaker for larger values of g 2 µL.
Additionally, the anisotropy is not strongly affected by the lattice resolution g 2 µaT . For each data point we performed an average over 200 random Glasma initial conditions. the magnetic field correlator is shown to be negative. As the ratio of the infrared regulator m and g 2 µ increases, the amplitude of the zero mode approaches zero from below. In the limit of m/(g 2 µ) → ∞, which corresponds to the dilute limit, the zero mode vanishes as seen in fig. 4 (b). Surprisingly, the electric field correlator also develops a region of anti-correlation for m/(g 2 µ) 0.1. This qualitative change from positive correlation for all distances r to regions with positive and negative correlation in the color-electric field perhaps explains the strong dependence of p 2 z on the infrared regulator m, see fig. 5. It is clear that in the dense Glasma, as well as in the dilute case, the initial chromo-electric and -magnetic fields have very different properties. More intuitively, fig. 7 (b) suggests that the shapes of average electric and magnetic flux tubes differ. The average chromo-magnetic flux tube exhibits a strongly pronounced "ring" of anti-correlation around its center where the sign of the magnetic field differs from the field at the center. In contrast, for a typical chromo-electric flux tube this effect is much less pronounced.
Although we are not able to establish a simple relation between the accumulated momenta and the momentum space correlators as in the dilute case, see eq. (155), we attribute the emergence of an anisotropy in transverse momentum broadening to the different correlation behaviour of the initial fields or rather the different shapes of chromo-electric and -magnetic flux tubes. Conversely, our calculation suggests that the observation of an anisotropy in the accumulated momenta of fast partons passing through the Glasma could be linked to differences in the properties of chromo-electric and -magnetic Glasma flux tubes.

V. CONCLUSIONS AND OUTLOOK
In this paper we computed the transverse momentum broadening of highly relativistic partons as they traverse the boost-invariant Glasma created directly after the collision of two high energy nuclei. Using a test particle approximation, we derived simple, gauge invariant expressions for the accumulated momenta of colored particles in a non-Abelian background field using both the Wilson loop definition and classical colored particle dynamics. This allowed us to study two interesting limits of the Glasma: the case of a dilute Glasma which can be described using the weak field approximation and the dense Glasma for which we have performed classical real-time lattice simulations.
In the case of the dilute Glasma we found that the accumulated transverse momenta can be directly computed from the correlation functions of the initial chromoelectric and -magnetic fields at τ = 0. This leads to an intuitive picture of momentum broadening of fast particles in terms of the Glasma flux tubes: we observe that momentum broadening within the transverse plane is caused by the Lorentz force from chromo-magnetic flux tubes while momentum broadening along the beam axis is due to chromo-electric flux tubes. Using the MV model as an initial condition for the Glasma, we found a large momentum anisotropy with much stronger broadening of momenta along the beam axis compared to broadening within the transverse plane. This anisotropy is caused by differences in the correlations of the initial fields: chromomagnetic Glasma flux tubes exhibit a large region of anti-correlation which strongly suppresses the growth of momentum broadening within the transverse plane after a characteristic time scale. This stands in contrast to chromo-electric flux tubes, which show no signs of such anti-correlated regions and consequently fast partons can accumulate more momentum along the beam axis. Remarkably, we also found that this momentum broadening anisotropy is not due to the longitudinal expansion of the system but stems purely from differences in the properties of the two types of Glasma flux tubes. Our real-time lattice simulation was also applied to the dilute case and found to be in very good agreement with the analytic results.
Finally, we studied the case of a dense Glasma using real-time lattice gauge theory which is phenomenologically more relevant for high energy nuclear collisions. We find that fast quarks accumulate up to p 2 ⊥ ≈ Q 2 s of total transverse momentum within proper times of τ 0 ≈ 6 − 10 Q −1 s . Here, τ 0 roughly corresponds to the switching time from the Glasma description to a hydrodynamical stage [52,53]. Translating this into appropriate units we can use saturation momenta in the range of Q s = 1 GeV to Q s = 2 GeV. Highly relativistic quarks acquire total accumulated transverse mo-  menta of p 2 ⊥ ≈ 1 GeV 2 to p 2 ⊥ ≈ 4 GeV 2 after traversing the Glasma for proper times of τ 0 ≈ 1.2 − 2 fm/c and τ 0 ≈ 0.6 − 1 fm/c respectively. For highly energetic gluons moving through the Glasma, the accumulated (squared) momenta are to be multiplied by a factor of C A /C F = 9/4. We highlight that our main results shown in fig. 5 allow us to estimate p 2 ⊥ for a large range of values for Q s , proper times τ and different values of the infrared regulator m.
Similar to the dilute Glasma, we find that momentum broadening in the dense Glasma is anisotropic with larger broadening along the beam axis. Here, we see that in particular the accumulated momentum along the beam axis p 2 z depends strongly on the infrared regulator m used in the initial conditions. In order to rule out that this anisotropy is simply an artifact of infrared regulation, we studied the limit of very large systems with vanishing infrared regulator and found that the anisotropy still persists.
We also studied the field strength correlation functions in the dense Glasma and found similar behaviour as in the dilute case. Chromo-magnetic flux tubes at τ = 0 exhibit pronounced regions of strong anticorrelation around r ≈ Q −1 s . Chromo-electric flux tubes on the other hand show much smaller anti-correlation around a larger length scale of r ≈ 2 Q −1 s . Presumably, these qualitative differences in the two types of Glasma flux tubes are causing the momentum broadening anisotropy. By reversing this argument we conjecture that highly energetic partons are able to probe differences in the properties of chromo-electric and -magnetic fields of the Glasma. Anisotropic momentum broadening of fast partons with stronger growth of momentum along the beam axis has also been found in the quarkgluon plasma with momentum anisotropy (see [15] and [18][19][20]). Our findings suggest that anisotropic broadening already occurs before the quark-gluon plasma phase. More work is necessary to see if these similar effects from two different stages of heavy ion collisions can be disentangled. It should be noted that the simple relation between the momentum broadening anisotropy and the correlators of initial fields has been derived only in the weak field approximation. Nevertheless, our simulations in the non-perturbative regime seem to suggest a similar relationship.
The methods presented in this paper can be extended in a number of interesting ways. There are two main directions that would improve our results: by using more realistic initial conditions for the Glasma and by using more accurate approximations for the interaction between the Glasma field and the partons moving through it.
The MV model used in this work is lacking a number of features that more sophisticated nuclear models like IP-Glasma [54,55] could provide to obtain more accurate estimations of jet momentum broadening from the pre-equilibrium stage. For example, the color charge density in the MV model is homogeneous on average while IP-Glasma accounts for variations due to the spatial distribution of nucleons within a nucleus. Additionally, the MV model as it was used in this work can only be applied to central collisions (i.e. zero impact parameter) and jets created in the center of the Glasma as there is no notion of a finite nuclear radius. Clearly, a relativistic parton created near the border of the Glasma, where the field strengths are much weaker compared to the center of the Glasma, would acquire less momentum from its interaction with the Glasma. Similarly, in off-central collisions the fields in the Glasma are much weaker leading to less momentum broadening overall. The MV model allows us to average over the starting position of the jet in the Glasma, but realistic initial conditions would require a more careful treatment of where and how particles are created in the collisions of high energy nuclei. It is also likely that near the border of the Glasma the transverse flow of the medium modifies the momentum broadening anisotropy (see [21] for details of such an effect in a hydrodynamical setting).
The calculation presented here is performed within the framework of a purely classical Yang-Mills field theory and does not account for corrections to the CGC initial conditions due to quantum fluctuations. To remedy this it would be necessary to use initial conditions evolved using the JIMWLK equations [56,57]. The JIMWLK evolution allows an effectively classical description of color sources including quantum fluctuations, which can then be used as Glasma initial conditions in place of the original MV model initial conditions. Earlier studies [51] showed that the gluon spectrum in the Glasma becomes harder (i.e. is shifted towards the ultraviolet) due to the JIMWLK evolution, leading to a larger average momentum per gluon in the Glasma. This could increase the amount of momentum accumulated by a jet as it moves through and interacts with the Glasma. Moreover, the JIMWLK evolution also allows for a rapidity dependent description of the Glasma [58], which would extend the purely boost invariant approach of this work.
Another way to remove the assumption of boost invariance is to account for the finite longitudinal width of the colliding nuclei along the beam axis. Numerical approaches to simulate this genuinely 3+1 dimensional scenario have been developed by the authors previously [59][60][61]. The methods for computing jet momentum broadening developed in this paper could be applied to the 3+1D case in a straightforward manner.
All of the above noted extensions are focused on changing the non-Abelian background field, but it is also important to improve the approximations used for the interaction of the particle and the background field. In this paper we approximated highly energetic partons as colored test particles moving at the speed of light, i.e. at infinitely high initial momentum. As a consequence of this approximation, the trajectory of the test particle does not change even though it acquires momentum along its path. Our results for momentum broadening therefore do not depend on initial momentum and should be interpreted as a limiting case of a more general calculation for finite initial momenta. The most obvious disadvantage of this approximation is that it does not account for energy loss of particles traversing the Glasma. A partial solution to this is to include perturbative, non-eikonal corrections by allowing for small deviations from the straight line trajectory. More generally, a simulation of dynamic test particles using the colored particle-in-cell (CPIC) method [43][44][45] could be used to numerically determine the full test particle trajectories. Going beyond the test particle approximation, it would also be interesting to include back reaction from the Glasma: as a colored particle moves through the Glasma, the Coulomb field of the particle interacts with the Yang-Mills field of the Glasma. As a consequence, the Glasma around the colored particle becomes polarized. This polarization leads to an induced chromo-electric field which would cause energy loss of the particle (see e.g. [15,62]). The effect of parton energy loss from polarization effects in the Glasma could be computed using linear response [63,64] or also using a fully dynamical CPIC setup.

VI. ACKNOWLEDGEMENT
We thank Carlos A. Salgado for suggesting to look into this topic and Kirill Boguslavski, Tuomas Lappi and Jarkko Peuron for helpful discussions. We also thank Anton Rebhan and Alexander Soloviev for comments regarding the manuscript. This work has been supported by the Austrian Science Fund FWF No. P32446-N27 and No. P28352. The Titan V GPU used for this research was donated by the NVIDIA Corporation. In this appendix we give a detailed derivation of the accumulated squared momentum p 2 i (t) R in eq. (38) using the Wilson loop definition in a non-Abelian background field. Analogous to the derivation in the Abelian case (see section II B), we consider a rectangular Wilson loop with the transverse extent aligned with the y axis and the lightlike extent aligned with x + (see fig. 2). The Wilson loop is given where the lightlike Wilson lines are Here, P denotes the path ordering symbol and g is the Yang-Mills coupling constant. The path ordering symbol orders terms according to their occurrence along the path. Quantities that come "later" in the path are to the left of quantities that come "before". For readability we suppress the dependence of A µ on x − and z. The definition of the Wilson loop W z+ with the transverse extent in the z direction is completely analogous. In order to perform the Taylor expansion in L we have to expand three of the four Wilson lines in eq. (A1) up to linear order in L to obtain W The expansion of the lightlike Wilson line W +,1 is slightly more involved. The linear coefficient reads where W + (b, a) is a lightlike Wilson line from x + = a to x + = b given by Putting everything together, the linear coefficient of W Next, we define the parallel transported field strength as F y+ (x + , 0) = W + (0, x + )F y+ (x + , 0)W + (x + , 0), (A20) and rewrite the expression for the linear coefficient: The accumulated squared momentum (see eq. (37)) where we have re-parametrized the integrals over x + as integrals over t along the particle trajectory x(t). We can perform the same derivation for the other transverse coordinate z. The general result reads .