Simulating jets and heavy quarks in the Glasma using the colored particle-in-cell method

We explore the impact of strong classical color fields, which occur in the earliest stages of heavy-ion collisions and are known as the Glasma, on the classical transport of hard probes, namely heavy quarks and jets. To achieve this, we simulate SU(3) color fields using classical real-time lattice gauge theory and couple them to an ensemble of test particles whose dynamics are described by Wong's equations. We provide an overview of how classical color algebras are constructed and introduce a method to generate random classical SU(3) color charges. We extensively test our numerical particle solver in the limits of infinitely massive heavy quarks and ultra-relativistic light-like jets and obtain excellent quantitative agreement with previous studies. Going towards realistic masses and initial moment, we extract longitudinal and transverse momentum broadening for heavy quarks and jets. The resulting accumulated momenta and the anisotropy of these dynamical hard probes exhibit deviations from limiting scenarios, showing that the full dynamics have a significant effect.


I. INTRODUCTION
Relativistic heavy-ion collision experiments, as conducted at the Large Hadron Collider (LHC) or the Relativistic Heavy Ion Collider (RHIC), provide the remarkable opportunity to study hadronic matter under extreme conditions with increasing statistics and precision. Immediately after the collision, the medium is characterized by large gluon occupation numbers and a highly nonlinear regime, known as the Glasma [1][2][3][4]. Particularly sensitive probes of the very early stage of the collision are heavy quarks and jets. Due to their short formation time, they experience the initial stage of the collision. By understanding the imprint of the Glasma fields on these probes, one can disentangle important information about the structure of initially produced matter, in both proton-nucleus and nucleus-nucleus collisions.
Previous approaches that focus on the effect of the Glasma on hard probes include a study on jets in the Glasma [23,24] based on a lattice discretization of the Yang-Mills equations, where the transport properties of jets are evaluated by treating them as ultra-relativistic light-like partons. More precisely, the jet momentum broadening is extracted from Glasma field correlators computed on the lattice, without explicitly solving the dynamical particle equations of motion. Another lattice study [18,25,26] with over-occupied Yang-Mills plasma instead of Glasma, evaluates the heavy quark transport coefficient from electric field correlators (assuming the heavy quarks to be infinitely massive and static) and emphasizes the emergence of plasmon mass induced oscillations. In another series [27][28][29][30][31][32][33][34], the effect of the Glasma phase on the diffusion of heavy quarks is extensively studied and compared to the standard Langevin description of heavy quark dynamics, with a recent focus on memory effects. A different approach is taken in [35][36][37], where both the Glasma fields and particle transport equations are derived using analytical frameworks. The Glasma fields are obtained in the proper time expansion and the transport of the hard probes is treated using the Fokker-Planck equations adapted to the Glasma. Complementary, it was shown that the initial stage, implemented in different frameworks, has an effect on jet quenching [38,39]. Even though these approaches vary with respect to the approximations which are used, they all converge to the same key result: the Glasma phase has a considerable effect on the transport of hard probes. Never-theless, very few of these studies have a built-in way to describe the very early stage consistently and in many cases they are constructed on approximations applicable at later stages.
In this work, we present a novel framework that simulates the full dynamics of hard particles right after the collision on top of an evolving boost-invariant SU(3) Glasma background field. This is practically achieved by developing a numerical solver for the equations of motion of particles propagating in these fields. The particles are initialized with finite masses, formation times and initial momenta. The solver is used to extract relevant quantities such as the momentum accumulated as the partons propagate in the background fields. The novelty consists in the numerical methods developed for the particle solver and the techniques used to efficiently solve both the Glasma and particle equations concurrently. In particular, we introduce a novel way to generate SU(3) classical color charges using the Haar measure. The code runs on GPUs and allows for the systematic study of the full dynamics of particles and the dependence on many parameters used for particle initialization.
There exist two relevant limiting cases in which the accumulated momentum of hard probes in Glasma may be evaluated only from Glasma lattice field correlators, without solving the particle equations of motion. These correspond to infinitely massive heavy quarks and highly energetic jets. When we consider such quarks in our particle solver, we reproduce the limiting results. The limiting case of extremely fast light-like jets is extracted using two setups, namely the classical transport framework using Wong's equations and a quantum pQCD computation. By comparing the resulting momentum broadening, we notice a discrepancy between the classical computation and the quantum one, and propose a way to resolve it. Going beyond these limiting cases, towards realistic dynamical results, we quantitatively study whether the full dynamics has a considerable effect. We extract the instantaneous transport coefficients, namely κ for heavy quarks andq for jets, and check if the large transport coefficients of hard probes in the Glasma obtained by previous studies are an artifact of the approximations used or still persists with our full numerical setup. Most remarkably, we observe that that momentum broadening along rapidity oscillates as a function of proper time, which could indicate plasmon modes in the Glasma [25]. Preliminary results obtained using our solver have been presented previously in [40].
This study is structured as follows. Section II contains an overview of the classical description of the early stage in terms of Glasma initial conditions and classical boost-invariant Yang-Mills equations. In Section III we present Wong's equations. In Section IV we describe how SU (2) and SU(3) color charges are sampled correctly. Section V describes the limiting cases for infinitely massive heavy quarks and for extremely fast light-like jets. In Section VI we show how to reconcile classical particle simulations with the calculation of momentum broaden-ing within pQCD. The results obtained with our particle solver are showcased in Section VII. Finally, Section VIII includes a summary of all the results, along with viable future extensions of our study. Detailed calculations including the correct sampling of color charges can be found in Appendices A through D.

II. GLASMA IN A NUTSHELL
Within the Color Glass Condensate framework [5][6][7], the medium produced after the collision of relativistic nuclei is a state dominated by strong classical color fields known as the Glasma [1-4, 8, 9]. The CGC is an effective theory for high energy nuclei and relies on the separation of scales between degrees of freedom with small and large longitudinal momentum fraction x. Hard (large-x) partons behave as highly Lorentz-contracted, static color sources J µ for the gauge fields A µ described by the soft (small-x) partons. At leading order in the coupling constant g, the hard and soft sectors are coupled via the Yang-Mills equations with D µ ( . . . ) ≡ ∂ µ ( . . . ) − ig A µ , . . . denoting the gauge-covariant derivative, F µν = ∂ µ A ν − ∂ ν A µ − ig A µ , A ν the field strength tensor and J µ the color current. At sufficiently high energies, we can approximate the nuclei to be propagating along the light-cone directions x ± ≡ (x 0 ± x 3 )/ √ 2. Their color currents are given by where ρ A,B represent classical color charge densities and the subscripts A and B denote the two colliding nuclei. The color charge densities are treated as stochastic variables whose statistics are determined by the probability functional W [ρ]. We take it to be given by the McLerran-Venugopalan (MV) model [41][42][43]. It considers the color charges ρ to follow Gaussian statistics, which are determined by the one-and two-point correlators where λ A,B is the average color charge per unit volume. One may extract µ 2 A,B = dx ± λ A,B (x ∓ ), which denotes the MV model parameter (in units of energy squared) and represents the variance of the color charge density fluctuations of each nucleus.
We can solve the Yang-Mills equations from Eq. (1) for the special choice of color current in Eq. (2) in the covariant gauge ∂ µ A µ cov = 0. The only non-zero components of the gauge field are given by in which ∆ ⊥ is the transverse Laplace operator. The Poisson equation can be formally solved via Fourier transformation where λ is an infrared regulator andρ cov A,B are the Fourier transformed charge densities in the covariant gauge. By performing a gauge transformation to the light-cone gauge A + lc = 0, the gauge field only has transverse components given by with the light-like Wilson line where P(. . . ) denotes the path-ordering operation.
In the ultrarelativistic limit, the nuclei are contracted to infinitesimally thin sheets. This can be expressed via where the two-dimensional charge densities obey the correlator The transverse gauge fields are given by in which θ represents the Heaviside function and involves a Wilson line depending on the transverse coordinate, obtainable as We now consider the classical collision problem where the initial conditions in the asymptotic past are provided by the color fields of the nuclei. The Glasma is described by the gauge field in the future light-cone of the collision. In the ultra-relativistic limit, the total color current generated by the two nuclei J µ = J µ A + J µ B possesses invariance under longitudinal Lorentz boosts, which implies that any observables of the Glasma must be invariant under boosts as well. An appropriate choice of coordinates is given by the Milne coordinates (τ, η) defined as with proper time τ and space-time rapidity η. By fixing the residual gauge freedom by imposing the temporal gauge condition A τ = 0 and requiring boost invariance of the gauge fields as A µ (τ, η, x ⊥ ) = A µ (τ, x ⊥ ), one may formulate initial conditions for the Glasma fields along the boundary of the future light-cone as [44] A accompanied by The conjugate momenta associated with the gauge fields are The Yang-Mills action expressed in Milne coordinates, together with boost-invariance, yields the field equations along with the Gauss constraint D i P i + ig [A η , P η ] = 0 which is fulfilled throughout the evolution. In order for these equations to preserve gauge invariance upon discretization, they need to be recast in a lattice QCD formulation.

A. Numerical implementation
The Yang-Mills equations of the Glasma may be solved numerically. The work presented here is based on an approach that employs classical real-time lattice gauge theory [8,45]. In order to assure gauge invariance of the field equations from Eqs. (16) upon discretization, one may proceed as follows: the Minkowski space is discretized on a hypercubic lattice replacing the gauge fields with gauge links, which are Wilson lines connecting neighboring points on this lattice. A particularity of the boost-invariant collision scenario is that one needs to employ this procedure only in the transverse plane. This is due to the fact that A η (τ, x ⊥ ) acts as a scalar under η-independent gauge transformations, and thus the η direction is left continuous with respect to a lattice discretization.
The transverse plane, taken as a square of length L and accompanied by periodic boundary conditions for the fields, is discretized in N 2 points in which the fields are assigned values at various proper times. In the continuum limit, assuming small lattice spacings a = L/N , a gauge link connecting the lattice point located at x ⊥ and the neighboring point along a directionî, whereî is the unit vector along x i , is given by Links in opposite directions can be expressed through the . These gauge links are then used to construct a plaquette variable as Uîĵ(τ, The Yang-Mills action can be approximated using link and plaquette variables along with the conjugate momenta Varying the discretized action yields discretized equations of motion where contains the forward D F i and backward D B i gauge-covariant finite differences on the lattice, ( . . . ) ah denotes the anti-Hermitian traceless part of a matrix, and represents the parallel transported scalar field. These equations are accompanied by the Gauss constraint and are solved numerically by employing the leapfrog algorithm. In this numerical method, the conjugate momenta are evaluated at half-integer time steps, whereas the rest of the fields are computed at integer proper times.
Finally, the MV model initial conditions must be discretized as well. The naive use of Eq. (8) leads to loss of randomness in the infinitesimal direction x ± along which the nucleus propagates, due to the non-trivial path-ordering of the involved Wilson lines. Nevertheless, by sticking together infinitesimally thin sheets of color charge and regularizing the correlator as [46] ρ a m ( where m, n ∈ {1, 2, . . . N s } denotes the index of the sheet, and with N s the number of such color sheets, this issue is resolved. Numerically, color charges are generated by sampling random numbers distributed according to a Gaussian with zero mean and variance chosen to obey Eq. (20). Once the color charges are provided, the solutions of Eq. (4) now expressed for each color sheet as ∆ ⊥ α a n ( x ⊥ ) = −ρ a n ( x ⊥ ) may be obtained using Fast Fourier Transformation (FFT), where the infrared and ultraviolet cut-offs are λ and Λ. Furthermore, the Wilson lines are constructed as products computed α n ≡ α a n T a . Subsequently, the transverse gauge links are computed from these discretized Wilson lines and the initial Glasma conditions given in Eq. (13) are also numerically discretized. Once all these steps are completed, the Glasma fields are numerically solved using our numerical simulation routines 1 .

III. PARTONS IMMERSED IN GLASMA
The dynamics of particles propagating in classical Yang-Mills fields is given by Wong's equations [47] which describe how the positions and momenta of the particles evolve in time, while their charges rotate in color space [48]. In the laboratory frame they read as where i ∈ {x, y, z} and a = {1, 2, . . . , D A } with the dimension of the adjoint representation D A = N 2 c − 1. The energy is given by E = p 2 + m 2 with m being the mass of the particle, p ≡ (p x , p y , p z ), and f abc the structure constants for the SU(N c ) group.
The dynamic equation for the energy p 0 = E is given by with E i ≡ F 0i denoting the color-electric field. This relation states that the energy of a moving particle changes due to the work exerted by the color-electric field upon it, where p = γm v, with γ = E/m the Lorentz factor and v the laboratory frame velocity.
Wong's equations may be recast into a covariant form, with quantities computed along the worldline of the particle dx µ dτ = p µ m , where (dτ ) 2 = g µν dx µ dx ν denotes the relativistic proper time and in which E d/dt = m d/dτ is employed and D/dτ is the covariant derivative taken along the particle worldline. Further, we make use of the Lie-algebravalued color charges Q = Q a T a to write Q a F µν,a = Tr [QF µν ] /T R , where T R is the representation-dependent Dynkin index defined through Tr T a T b = T R δ ab , with T a ∈ SU(N c ). The above equations simplify to The equation governing the evolution of the color charge may formally be solved by where Q is rotated with Wilson lines. These Wilson lines involve the path-ordered exponential computed along the trajectory of the particle and are given by The last relation may be derived by making use of the parallel transport equation for a Wilson line d dτ The use of Wilson lines in the evolution of the color charge automatically conserves the quadratic Q a Q a ≡ q 2 (R) (28) and cubic classical Casimirs 2 where d abc are the symmetric structure constants and the values q 2 (R), q 3 (R) depend on the chosen representation R for the color charge Q. We go into detail about how these invariants are fixed in Section IV and Appendix B.
In this work we approximate hard partons as test particles, which means that we neglect any back reaction of the partons onto the Glasma.

A. Dynamics of particles in Glasma
Let us express Eqs. (23) and (25) in Milne coordinates and choose the background fields to be those of the boostinvariant Glasma. A detailed derivation can be found in Appendix A 1. The coordinate Wong equations are given by where x µ ∈ {x, y, η} and (p τ ) 2 = (p x ) 2 +(p y ) 2 +τ 2 (p η ) 2 + m 2 . The Wong equations for momenta read as and are accompanied by the dynamic equation for the temporal component where the color-electric and -magnetic fields are determined from the field strength tensor via As for the proper time evolution of the color charge, the Wilson line involved in the color rotation, see Eq. (26), can be expressed as a path-ordered integral along the worldline As used in the Glasma framework, we employ the temporal gauge A τ = 0 and the gauge field is taken to be independent of space-time rapidity η, which simplifies the Wilson lines.
As colored particles pass through the Glasma, the momentum of the particles p µ changes according to Wong's equations. The main observable we focus on, which represents a measure of the accumulated momentum, is the momentum broadening δp µ defined as Here τ form denotes the formation time at which the particle is introduced into the system and p µ (τ form ) is the initial momentum of the particle. The momentum broadening thus reflects how much momentum is accumulated through interactions with the Glasma background field compared to the initial momentum of the particle.

B. Numerical implementation
The positions and momenta of the partons are initialized using a toy model setup. Namely, the initial coordinates of the particles, chosen at formation time τ form ≥ 0, are randomly distributed in the transverse plane x(τ form ), y (τ form ) ∈ [0, L] with η(τ form ) = 0. The particles initially only have transverse momenta p T = (p x ) 2 + (p y ) 2 at formation time, with fixed p T (τ form ) and p η (τ form ) = 0 for heavy quarks, or an initial p x (τ form ) and p y,η (τ form ) = 0 for jets propagating along the x-axis. As will become evident in Sec. V, we choose jets with initial momenta along x-direction in order to compare with previous studies having the same particle setup. Color charges are randomly sampled using Darboux variables for SU (2) or using the Haar measure for SU (3), and their associated classical Casimirs are fixed according to Eqs. (44a) and (44b). A complete description of how these classical color charges are constructed is given in Section IV and Appendix B.
Numerically, the Milne proper time evolution for positions and momenta from Eqs. (30) along with (31) and also (32) for the temporal constraint is solved with Euler's method. An example of the numerical solutions of these equations for particles propagating in Glasma fields is depicted in Fig. 1. The Glasma electric and magnetic fields from Eq. (33), which appear in Wong's momenta equations given in Eq. (23), have to be approximated on the lattice. This is because the electric fields reside on gauge links, the magnetic ones on plaquettes and we need to interpolate in order to get their value on a lattice site. Appropriate approximations that are accurate up to quadratic order in the lattice and time spacing are given by where i = x, y, and similarly When evaluating these expressions, the positions of the particles are approximated with the nearest grid point (NGP) on the transverse lattice of the Glasma. Thus, the color-electromagnetic fields that are exerted on the partons are computed at The numerical solution for the rotation of the color charge is more involved and relies on the same NGP approximation. It is inspired by colored particle-incell (CPIC) methods used in the context of particles in CYM plasmas [49][50][51][52]. In the CPIC method, the color charge of a particle is rotated with gauge links only when the NGP on the underlying simulation lattice of the background YM fields changes. It should be emphasized that the Glasma fields are discretized only in the transverse plane because of boost invariance, and the rapidity direction is left continuous. Thus, one needs to adapt the CPIC method to the Glasma lattice discretization with gauge links only in the transverse plane. Numerically, one may approximate the Wilson line from Eq. (34), namely U(τ i , τ f ) at a given proper time τ f as being comprised of subsequent products of "short" Wilson lines U(τ n−1 , τ n ) as These short Wilson lines may be reduced to  Here, U xn,î (τ n ) is a transverse gauge link along the directionî with i = x, y evaluated at position x n , while U xn,η (τ n ) represents a Wilson line along theη direction, which can be computed from A η via the matrix exponential. It should be noted that this approximation is only valid for small time steps δτ n = τ n − τ n−1 .
We have made use of the fact that the displacement in rapidity δη n is numerically small and it follows that dx i A i , δη n A η 0 such that higher order terms arising from the Baker-Campbell-Hausdorff formula are suppressed. This numerical color rotation is depicted in Fig. 2.
Alternatively, one may directly solve where the transverse gauge fields are numerically extracted from the gauge links using matrix logarithms igaA y x, y + a 2 = ln(Uŷ(x, y)).
We checked that these two distinct methods for solving the evolution of the color charge, either from Eqs. (25) or (40), are consistent with each other in the limit of small time steps and yield similar final results for momentum broadening. Nevertheless, the advantage of performing numerical color rotations with Wilson lines as in Eq. (25) lies in ensuring that the color charge remains in the Lie algebra, i.e. Q ∈ su(N c ), and that the Casimir invariants are exactly conserved. The Casimirs, Eqs. (28) and (29), remain unchanged throughout the evolution: once the values of the Casimirs are fixed at formation time, color rotations with Wilson lines will not affect them.

IV. CLASSICAL COLOR CHARGES
In the previous sections we have outlined how to numerically solve the field and particle equations on a lattice. The question remains how to choose the classical color charges Q in the ensemble of partons. Here, we largely follow the seminal works on classical non-Abelian transport theory [20][21][22]53]. There are three aspects to consider: first, what values to assign to the classical Casimir invariants of the color charges from Eqs. (28) and (29); second, how to distribute the charges in color space (the particular color charge a parton assumes after a random hard scattering is a priori unknown, hence additional considerations are required in order to construct their distribution); third, how fixing one affects the other. We address these by treating the color charge components Q a as stochastic variables with fixed values of q 2 = Q a Q a and q 3 = d abc Q a Q b Q c . We emphasize that this is a choice and in our framework, where the Casimirs remain constant throughout the evolution, see Sec. A 2, the obvious choice is to fix the Casimirs. In analogy with the trace relations for operator-valued elements of the  Figure 2. Diagram with the color rotation performed during a numerical time step from τn−1 to τn. The electric and magnetic Glasma fields reside on lattice points in the transverse plane x(τn) ≡ xn, while a particle may move at any location in the transverse plane. The particle position is approximated with the NGP on the lattice x(τn) → NGP(τn) and when the transverse coordinates of the NGP change, one performs a color rotation with the corresponding transverse gauge link, in this case Ux. Along the rapidity direction, a Wilson line Uη is computed via the matrix exponential and used in the color rotation, which here is simply given by U(τn−1, τn) = Ux(τn)Uη(τn).

su(N c ) color algebra
Tr Q a = 0, we choose to have color charges randomly distributed according to one-, two-and three-point functions where the representation-dependent coefficients T R and A R are given by The Casimir invariants from Eqs. (28) and (29) constrain what values the color charges Q a can take. It may be shown that the ansatz for the two-and three-point functions from Eqs. (42b) and (42c) fixes the quadratic and cubic classical Casimirs, defined in Eqs. (28) and (29), to We point out that assigning the labels "fundamental" and "adjoint" to the classical Casimirs is inspired by the corresponding quantum representations and is inherited from the choice we made in Eqs. (42). A detailed derivation of the classical Casimirs from Eqs. (44a) and (44b) is given in Appendix B 2. Here, we provide a sketch of the derivation. If we take the two-point function from Eq. (42b), choose the color component a = b and perform a sum over it, we get the classical quadratic Casimir q 2 (R) = D A T R . Similarly, we start from the threepoint function in Eq. (42c), multiply by d abc and sum over all color indices. Eventually, with the normalization we chose in Eqs. (42), the quadratic and cubic classical Casimirs may be recast in the following form where D R denotes the dimension of the representation, namely and C 2,3 (R) are the group-theoretical quadratic and cubic Casimirs, given here for the fundamental and adjoint representations as Combining Eq. (45) with the definitions in Eqs. (46) and (47)  With the chosen normalization from Eqs. (42), the resulting classical Casimirs differ from the grouptheoretical ones, see Eq. (45). This immediately raises the question whether the factor of D R in Eq. (45) can be absorbed in the normalization of the color charges from Eqs. (42), such that the classical Casimirs of Q a automatically match the group-theoretical ones q 2,3 → q 2,3 /D R = C 2,3 . We found that this is not always possible, i.e. the classical Casimirs do not coincide with the quantum ones for any gauge group or representation. In particular, when we consider quarks in SU(3), we were not able to find a color charge vector Q a which satisfies Eqs. (28) and (29) with q 2,3 (F ) = C 2,3 (F ). Although we do not have a formal proof, we believe that there exist no solutions to the color charge constraints for these particular values of q 2 and q 3 . In Sec. VI we discuss in more detail how this difference between classical and quantum color charges affects some particular expectation values, for example the momentum broadening δp 2 defined in Eq. (35), and how to address it.
The initial random classical color charges of the partons at formation time must satisfy the above relations in order to describe the physics of heavy quarks, and jets of quarks and gluons. In the following two subsections we show how random color charges satisfying the above n-point functions can be numerically realized for SU (2) and SU(3).

A. SU(2) classical color charges
For generating SU(2) classical color charges, we rely on the Darboux variables parametrization [22,54]. One may generically construct the classical limit of any semisimple Lie algebra [55]. This is done by starting from the defining commutation relations where f abc denote the structure constants and { Q a } is the set of operator-valued generators. Taking the reverse of the quantum limit, the classical correspondent is given by the Poisson bracket If one interprets the generators as classical variables depending on the symplectic structure of the underlying manifold through some phase-space coordinates (φ i , ξ i ), the Poisson brackets may be expressed as The pair of conjugate variables obey the canonical Poisson bracket relations {φ i , ξ j } = δ ij and are called Darboux variables. For SU(2), whose generators {Q a } with a ∈ {1, 2, 3} obey Eq. (50) with f abc = abc , one identifies a single pair {φ, ξ} and the subsequent phasespace evolution is restricted to conserve the quadratic Casimir from Eq. (28).
Simply distributing the color charges uniformly on a three-dimensional sphere of fixed radius J 2 ensures that Eq. (50) is satisfied and that the Casimir is fixed by q 2 = J 2 . The SU(2) color charges are sampled according to the parametrization where φ ∈ [0, 2π) and ξ ∈ [−J, J] are uniformly distributed random numbers.

B. SU(3) classical color charges
Similar to SU(2), one may construct a parametrization for classical SU(3) color charges in terms of the Darboux variables [54]. Unfortunately, any parametrization of classical color charges only covers a portion of the underlying manifold of SU(3) [55], leading to ill-defined one-, two-and three-point functions that will differ from the expected ones given in Eq. (42). For this reason, we rely on a different method to sample them, namely through the Haar measure of SU (3). The main idea is that the integration over color charge configurations may be mapped to integration over the underlying manifold of the group. This is done by first constructing an initial color vector Q 0 = Q a 0 T a such that the quadratic and cubic Casimirs Q a 0 Q a 0 and d abc Q a 0 Q b 0 Q c 0 satisfy Eqs. (44a) and (44b). The exact choice of Q 0 is arbitrary, as long as the Casimir invariants q 2,3 (R) match the desired values. Once the initial color vector is fixed, random color charges are generated by performing color rotations as with a random special unitary matrix U ∈ SU(3) distributed according to the Haar measure. From this color vector Q, color components are given by projecting onto the generators T a where we introduced the adjoint representation matrix By construction, the quadratic and cubic Casimirs are invariant to these color rotations.
We have to verify whether replacing the integration over classical SU(3) color charges with that over the SU(3) group elements as dQ → dU yields equivalent results. For this purpose, it suffices to check that the n-point functions of the color charges computed with the Haar measure along with the classical Casimirs fixed by Eqs. (44a) and (44b), exactly match the n-point functions of the classical colors charges from Eq. (42), namely This can be explicitly checked by carrying out the required integrals for SU (3). A detailed calculation can be found in Appendix C. In particular, we show that the n-point functions become independent of the initial color charge Q 0 , except for the values of the two Casimirs q 2 and q 3 . In this way, the generation of classical color charges for SU (3) can be replaced by sampling over the Haar measure.

V. LIMITING CASES
In general, the dynamics of colored particles passing through Yang-Mills background fields are non-trivial and can only be solved numerically, e.g. with the methods introduced in earlier sections. However, there are certain limiting cases where the dynamics become trivial and observables such as the momentum broadening defined in Eq. (35) can be reduced to simple functionals of the background fields. These limiting cases are those of infinitely massive heavy quarks and highly energetic jets. In both cases, the particle trajectories are trivial in the sense that the particles are not deflected by the forces acting on them.
These cases are of interest since there exist numerous studies which rely on the infinitely massive heavy quark approximation, with momentum broadening and diffusion coefficient κ extracted from electric fields correlators computed on the lattice [16,19,25,56], and the highly energetic jet scenario, with accumulated momentum and transport coefficientq related to light-like Wilson loops [11,17,23,24,57]. Moreover, they represent valuable numerical checks for our particle solver which, in these limiting cases, should give similar momentum broadenings as those extracted solely from Glasma fields.
Using the formal solution for the evolution of the color charge from Eq. (25), one may recast Wong's equation in the Milne frame from Eq. (24) in the following form where denotes the parallel transported color Lorentz force with and the initial color charge vector is expressed as (57), the momentum broadening from Eq. (35) may be written as where no sum over µ is implied. Using the twopoint function of the color charges chosen according to Eq. (42b), along with the Fierz identity expressed as valid for traceless N c × N c complex matrices, we arrive at the formal solution for the momentum broadening A. Infinitely massive heavy quarks An infinitely massive heavy quark is static and remains at rest in the Milne frame. Due to temporal gauge, all temporal Wilson lines are unity U(τ, τ ) = 1. Therefore, no parallel transport is required, thus F µ = F µ according to Eq. (58). Furthermore, in the infinite mass limit m → ∞, the temporal component of the four-momentum simply behaves as p τ → ∞. Thus, the Lorentz force contains only contributions from the electric fields The momentum broadening for static particles thus reduces to an integral over electric field correlators where no sum over i is implied and the fields are evaluated at some fixed transverse coordinate. This expression can be evaluated purely from color-electric fields, without the need to solve the dynamical particle equations of motion.

B. Highly energetic light-like jets
The case of a highly energetic jet moving through the Glasma has already been studied in [23,24] and here only the final results are quoted, as the derivation is analogous to the case of static particles. The momentum broadening of a light-like parton traveling along the x-axis in Glasma fields is given by since for jets we assume τ form = 0. The various components of the Lorentz force are evaluated using the Glasma color electric and magnetic fields as These color field components have to be parallel transported according to using a Wilson line constructed along x as The integration over initial color charges is replaced by the Haar measure over SU(N c ), with particular choices for the Casimir invariants, as in Eqs. (44a) and (44b), and the functional integration over the background field is an average over the Glasma initial conditions encoded in a probability functional W [A µ ]. In order to reproduce the correct physics using classical calculations, these classical expectation values should match quantum expectation values computed for example in pQCD in the limit where the classical approximation is appropriate. Due to the overoccupied nature of the gluon field in the early stages of the collision, this approximation is valid for the background field, but in a strict sense fails when we approximate quarks and gluons as classical color charges. The problem is that quarks and gluons are low dimensional representations of the color algebra, whereas classical color charges are obtained in the limit of high dimensional representations (for example, in the case of SU(3) quarks, we found no classical color charges whose Casimirs coincide with the group-theoretical ones, see Sec. IV and Appendix B). However, as we shall see below, the classical framework can nevertheless reproduce quantum expectation values for certain observables of interest, such as the momentum broadening defined in Eq. (35), by constructing a meaningful quantity whose classical expectation value correctly gets mapped to its quantum version.
To establish this relationship, we focus on the case of a light-like parton moving along the x + -axis in the eikonal approximation, i.e. in the case where the trajectory is fixed. Within pQCD and holography, momentum broadening may be related to particular Wilson loops [11,14,57,58]. In particular, momentum broadening orthogonal to the trajectory may be evaluated from a rectangular Wilson loop with one side parallel to the trajectory (light-like extent L) and the other side chosen to be spatial and orthogonal to x + (transverse extent L ⊥ ). In the small transverse extent limit L ⊥ → 0 one finds where p 2 i R is the momentum broadening of a parton in representation R along the transverse directionî. The expectation value is taken over an ensemble of background fields. It is relevant to notice the factor 1/D R in front, which ensures that the identity holds for L ⊥ → 0, when the Wilson loop reduces to W i+ (R) → 1 D R . This factor will play an important role when matching with the classical computation. Performing a Taylor expansion in L ⊥ , where the Wilson loop is written as and inspecting the second-order coefficient yields the momentum broadening [23] From this relation, one expects Re Tr [ .
In the case of classical background fields such as the Glasma, the second order coefficient can be written in terms of the field strength tensor via where F i+ (τ ) ≡ F i+ x(τ ) denotes the parallel transported field strength tensor which is given by and contains the light-like Wilson line On the other hand, we can compute the same expectation value within the classical particle framework. For a light-like trajectory x + = √ 2t, see also Eq. (65), the momentum broadening is given by Moreover, as may be inferred from Eq. (60), one expects (76) suggests that the classical computation may be mapped to the quantum one by considering or equivalently, the classical expectation value coincides with the quantum one after a division by the dimension of the representation, since D R = q 2 (R)/C 2 (R) according to the classical Casimirs from Eqs. (45). An analogous calculation can be performed for infinitely massive partons, which involves a time-like Wilson loop instead of a light-like loop. Repeating the same steps, we arrive at the same factor of D R to match the classical to the quantum expectation value. We note that this division is also performed in [14,59]. It should not be surprising that calculations based on classical colored particles do not entirely match pQCD calculations since, already on a formal level, there is an important difference between the high-dimensional classical and the low-dimensional quantum representations (see Appendix B), namely they are not labeled by the same Casimir invariants. Our choice given in Eq. (45) shows that the classical Casimirs q 2,3 (R) are the dimension of the representation D R times the group-theoretical ones C 2,3 (R), which comes from how we choose to distribute the classical color charges, see Eq. (42). Thus, the origin of the difference between the classical and the quantum expectation values may be traced back to the statistical properties of the ensemble of classical color charges. Specifically, δp 2 is directly related to the classical two-point function Q a Q b , see Eq. (62), and thus, the quadratic classical Casimir δp 2 classical ∝ q 2 . On the other hand, the quantum correspondent satisfies δp 2 quantum ∝ C 2 . Therefore, when mapping classical to quantum expectation values, the meaningful quantity to compare is actually δp 2 /C 2 , where C 2 denotes either the classical or quantum quadratic Casimir, as stated in Eq. (77).
Moreover, a similar argument which leads to Eq. (77) also works for δp 3 , namely δp 3 /C 3 is the correct quantity to map from classical to quantum, where C 3 denotes the classic or group-theoretical cubic Casimir. Nevertheless, it fails for δp 4 or higher-order moments. As shown in [60], where the computation of the averages over the classical color charges is performed for SU(2), such a matching fails for the four-point function of the gauge field when compared to the 1-loop quantum effective action. It is only in the limit of high dimensional representations, where such a matching is exact. More concretely, the previous arguments generalize to δp n classic ∝ Q a1 . . . Q an . In analogy with Eq. (42), we choose where T (a T b) denotes the symmetric part of T a T b . As previously shown, such relations are satisfied for n = 1, 2, 3 with Q a obeying the classical Casimir constraints in Eqs. (45), but are violated for n ≥ 4. For consistency with pQCD calculations, the classical framework is thus limited to the quadratic and cubic moments of the momenta in a strict sense. There is another important property of highly energetic jets that suggests that the matching condition in Eq. (77) is appropriate, namely Casimir scaling. The ratio of the accumulated momentum of the adjoint and fundamental representation must yield the ratio as was already noted for Eq. (72). This scaling with the ratios of the Casimir invariants is observed in many systems. For example, the Casimir scaling of the transverse momentum broadening coefficientq is a result inherited from pQCD computations of partons in weakly-coupled QGP [12,13] or in weakly-coupled N = 4 SYM [61] and holds in the eikonal limit as in Eq. (73). Moreover, it is also a direct consequence of Wong's equations and the properties of the color charges provided that we use the matching condition in Eq. (77). To see this, we first rewrite Eq. (76) into This expression states that the classical accumulated momentum for a colored parton in representation R is simply proportional to the representation-dependent factor T R . Consequently, since T A /T F = q 2 (A)/q 2 (F ), the classical accumulated momenta behave as which resembles the Casimir scaling of Eq. (79) but in terms of the classical Casimir from Eqs. (28). The division by D R of the classical momentum broadening, see Eq. (77), restores it to a Casimir scaling with grouptheoretical Casimirs for a given SU(N c ) group. Such a relation can equivalently be seen from the mapping proposed in Eq. (77). We checked, with our particle solver, that the grouptheoretical Casimir scaling for momentum broadenings divided by D R is satisfied throughout the evolution, see Fig. 10 and the discussion in Appendix D.

VII. RESULTS
In this section we apply the previously developed numerical methods to study momentum broadening of heavy quarks and jets in the early Glasma stage of heavyion collisions. To gain trust in our methods, we first compare our particle simulations to the limiting cases of infinitely heavy quarks and infinitely energetic jets. In Subsection VII C we proceed with realistic simulations of dynamical heavy quarks, such as charm and beauty. Similarly, in Subsection VII D realistic jet momentum broadenings are extracted and the jet transport coefficient is computed.

A. Choice of parameters
For the Glasma, the saturation momentum is chosen as Q s = 2 GeV, while the MV model parameter is fixed through g 2 µ ≈ 0.8 Q s for N s = 50 color sheets and the IR regulator as m = 0.1 g 2 µ, according to [62], and the UV regulator as Λ = 10 GeV. The coupling constant is evaluated from the running coupling constant as g 2 = 4πα s (Q s ) computed at a given saturation momentum with N f = 3 and Λ QCD = 200 MeV, which yields g ≈ 2.07. The rest of the numerical parameters of the Glasma are set as follows: the length of the simulation domain in the transverse plane is L = 10 fm, the number of lattice points is N = 512 for heavy quarks or N = 1024 for jets. The time step ∆τ , which is used in the leapfrog scheme for the Glasma fields, is given in terms of the transverse lattice spacing a ⊥ = L/N : for heavy quarks we use ∆τ = a ⊥ /8 and for jets we use ∆τ = a ⊥ /16. The numerical code for this work is an extension of an earlier Glasma code used in [23,24] and is hosted publicly [63]. As initial conditions for classical particles, we rely on a toy model initialization of positions and momenta. Namely, all partons are randomly distributed in the transverse plane, at mid-rapidity, and have a fixed initial transverse momentum. For heavy quarks, their formation time is given by τ form ≈ 1/(2m HQ ), with m charm = 1.27 GeV and m beauty = 4.18 GeV [64]. All jets are formed instantaneously, at the same time as the Glasma fields. A single Glasma event contains N tp = 10 5 test particles and most of the results are obtained for N events = 30 Glasma events, although convergence was reached for fewer events. The transverse simulation region has periodic boundary conditions for the particles, whereas the rapidity direction is left continuous. We emphasize that the nuclei we simulate are not finite in size (their realistic geometry is not taken into account) and occupy a square lattice in the transverse plane, with periodic boundary conditions. Moreover, expecting that the Glasma picture holds up to τ 0.3 fm/c, the details of the geometry and the transverse expansion are expected to be less relevant in the extraction of less sensitive quantities, for example the momentum broadening.

B. Comparison with limiting cases
In the limit of infinite particle mass m → ∞ (infinitely heavy quarks) or infinite spatial momentum p x → ∞ (infinitely energetic jets), the dynamics of the particles become trivial and the accumulated momenta reduce to Eqs. (64) and (65) respectively. In order to validate our simulations, we compare results from the limiting cases (which were already used in [23,24]) to simulations with classical particles in these particular limits. Our results are shown in Fig. 3. For heavy quarks, we show the transverse δp 2 T = δp 2 x + δp 2 y and the longitudinal δp 2 L momentum components. We note that "longitudinal" refers to the component along the beam axis, whereas "transverse" denotes the orthogonal direction. For jets, we show all three independent components δp 2 i with i ∈ {x, y, z}. In both cases, the longitudinal momentum broadening increases much faster than the transverse one at early times and reaches a maximum around Q s τ ≈ 10. Peculiarly, the longitudinal component for heavy quarks undergoes multiple damped oscillations before settling to a constant value at late times. In contrast, the longitudinal component for light-like jets has only a single pronounced peak and then quickly saturates. For both heavy quarks and jets, the transverse momentum broadening components increase rapidly at early times, due to strong coherent fields, but become essentially constant within Q s τ 10. Similar phenomena have been observed in [25], where heavy quark diffusion was studied in overoccupied gluonic systems without expansion. In these systems, the  As is evident from our data, both approaches yield the same results to a large degree. The slight numerical difference is due to the fact that for the limiting case result in terms of field correlators, we discretize over time the integrals in Eqs. (64) and (65) in steps of the transverse lattice spacing a ⊥ . For our particle solver we typically use much smaller time steps ∆τ a ⊥ leading to a slightly more accurate result. More concisely, we numerically checked that reducing the lattice spacing used in the particle solver (in order to make it "less accurate") lead to a better agreement with the limiting case result.

C. Heavy quark momentum broadening
Having established that our particle simulations correctly reproduce limiting cases, we can now focus on more realistic simulations of heavy quarks with finite masses and finite formation times. Moreover, we can use our simulations to extract the heavy quark transport coefficient, which we define as This is the instantaneous heavy quark coefficient and may be interpreted as a diffusion coefficient in the limit of large proper times, namely κ diffusion = lim τ →∞ κ inst (τ ).
Our results for beauty quarks with vanishing initial transverse momentum are shown in Fig. 4, where we plot the accumulated momenta and their time derivatives. As in the case of infinitely heavy quarks, the longitudinal momentum broadening component δp 2 L increases more rapidly than the transverse component δp 2 T at early times. Even though not shown here, we checked that the longitudinal and transverse momentum broadenings have the same behaviors at larger proper times τ 2 fm/c, as already noticed in Fig. 3 for static quarks. The first peak of the oscillations in δp 2 L happens at around δτ = 0.8 fm/c, giving rise to a temporary negative heavy quark diffusion coefficient κ L . In contrast, the transverse component approaches a constant value after δτ ≈ 0.5fm/c. A qualitatively similar picture emerges for charm quarks.
In general, the accumulation of momentum of heavy quarks depends not only on their mass (and thus formation time), but also their initial transverse momentum p T . Figure 5 shows the numerical results for beauty and charm quarks for various values of the initial p T ∈ {0, 2, 5, 10} GeV 3 . For comparison, we include the static quark limit as a dashed curve. Since the Glasma affects the heavy quarks in an anisotropic manner, we also plot the heavy quark anisotropy coefficient, which we define Beauty quarks, due to their early formation time, experience the initial strong and coherent Glasma fields more than charm quarks. For this reason, their momentum broadening is generally larger than that of charm quarks.
On average, beauty quarks acquire 30-50% more momentum than charm quarks. A similar observation was also emphasized in [33].
Focusing on the heavy quark anisotropy, we find that as the initial p T increases, δp 2 L decreases and δp 2 T increases. Consequently, the corresponding anisotropy δp 2 L / δp 2 T becomes smaller. Compared to the static quark accumulated momentum (dashed lines), beauty quarks with zero initial p T have an increase in δp 2 L of 50% and charm quarks 50-80% throughout the proper time evolution. For the maximum initial p T taken in our simulations, δp 2 L for dynamic quarks differs from that for static quarks by 50% for beauty and 30-70% for charm quarks, whereas δp 2 T increases only by 20-30% compared to static quarks. These will have a complementary effect on the anisotropy. Namely, δp 2 L / δp 2 T for beauty or charm quarks is 20-40% larger or smaller, depending on the initial p T , than that of infinitely massive heavy quarks formed at the same formation time. The anisotropy is higher for small initial p T heavy quarks and lower for quite large initial p T . The anisotropy is more pronounced for the zero initial p T heavy quarks. Therefore, there are slight differences between static quarks that "see" only the Glasma electric fields, and quarks initialized with vanishing momentum but allowed to move in the Glasma.
Understanding the dynamics of heavy quarks in the Glasma in terms of the electric and magnetic color fields is generally not trivial, but some of their properties may be inferred from particular features of the background Glasma fields. Initially at τ = 0 fm/c, the Glasma consists of correlated domains of longitudinal color-electric and -magnetic flux tubes with a typical size of ≈ 1/Q 2 s . This shortly lived initial phase is probed by heavy quarks with very high mass m HQ due to their early formation time τ form = 1/(2m HQ ). These heavy quarks are accelerated due to the strong longitudinal color-electric fields of the Glasma, leading to the rapid increase of the longitudinal momentum broadening component as seen in Figs. 4 and 5. If heavy quarks have a non-negligible initial transverse momentum, there is additional transverse acceleration due to longitudinal color-magnetic flux tubes. This effect, albeit small, is seen in Fig. 5 for both beauty (left panel) and charm (right panel) quarks, where the transverse momentum broadening component increases with the initial transverse momentum. Remarkably, the opposite occurs for the longitudinal component: larger initial transverse momentum leads to reduced longitudinal broadening, but we have found no simple explanation in terms of the field structure of the initial Glasma for this effect.
Immediately after their initial formation, the flux tubes start to expand in the transverse plane, which generates transverse color-electric and -magnetic field components. These transverse electric fields lead to a slightly delayed increase of the transverse momentum broadening component. At the same time, the highly correlated regions within the Glasma are lost, and the longitudinal acceleration becomes less efficient. Eventually, the Glasma transitions to the free-streaming regime at around τ free ≈ 1/Q s ≈ 0.1 fm/c, after which the fields become more dilute and the mean energy density falls off as 1/τ . As seen in Fig. 4(b), the heavy quark diffusion coefficient κ has already peaked by then and falls off quickly. The formation time of heavy quarks has a large influence on the accumulated momenta in the Glasma stage. As can be seen from Fig. 5 (right panel), charm quarks accumulate less momentum because they "skip", at least in part, the initially highly correlated phase of the Glasma at τ Q −1 s . Even though in our current setup we initialize heavy quarks homogeneously in the transverse plane, it is more likely that partons are formed inside the Glasma flux tubes, where the energy density is larger, and thus the particle production is more favorable. Thus, for illustrative purposes, we look at trajectories of almost static or dynamic beauty and charm quarks, initialized in the "center" of such a Glasma correlation domain, where the energy density reaches its maximum value.
The results are shown in Fig. 6, where the different colors of the trajectory lines correspond to various values of initial p T ∈ {0, 2, 5} GeV and the background shows the energy density at the formation time of the corresponding heavy quark. Almost static quarks with very high mass barely move during the evolution and thus essentially remain where they were originally produced at formation time. On the other hand, quarks with realistic masses are able to move further and probe larger spatial regions of the Glasma. Moreover, the quark mass determines when the particles are being introduced into the system and what regime of the evolution the partons are able to "see". For example, as shown in Fig. 6, charm quarks are produced close to the transition to the free-streaming regime, where the color flux tubes already started to expand. Slow heavy quarks spend more time in the correlation domains before they expand, whereas fast quarks escape them more quickly, and thus lose the correlation faster. Even though the picture of heavy quarks probing the Glasma correlation domains as illustrated in Fig. 6 describes an over-simplified scenario, it still offers a valuable qualitative understanding. Moreover, within the approximations we use for particle initialization, it offers hints that beauty quarks might be more viable probes of the Glasma than charm quarks.

D. Jet momentum broadening
In recent years, jets in the Glasma have been investigated using classical simulations [23,24] and the small τ expansion [36,37]. In all of these works, the initial energy of the jet has been assumed to be very large, such that the trajectory can be approximated as essentially lightlike. Since we account for particle dynamics via Wong's equations, we can use our particle solver to go beyond the light-like jet case and consider the effect of finite initial momentum along the propagation axis and different jet masses. For simplicity, we choose the jets to be initialized with finite p x values.
Similarly to the heavy quark transport coefficient κ, we distinguish between various components of the jet transport coefficientq. We define the instantaneous jet broad-  ening coefficientq with i ∈ {x, y, z}. This is different from the collisional energy loss dE/dx. Since the jet propagates along the x-axis, we introduce the transverseq T ≡q y and longitudinalq L ≡q z jet transport coefficients.
In addition, we are interested in deviations from the light-like jet scenario by considering a finite jet mass m and an initial p x . We find that jet momentum broadening essentially only depends on the ratio p x /m. In the limit of p x /m → ∞, we can compare to the limiting case given by Eq. (65). Similarly to the heavy quark anisotropy, we also introduce a measure of how the Glasma anisotropy affects the jets by defining the ratio Our numerical results for jets are shown in Figs. 7 and 8. Figure 7 shows similar behavior as in the case of heavy quarks. After reaching a maximum at roughly τ ≈ 0.8 fm/c, the longitudinal component δp 2 z starts to decrease at late times τ 1 fm/c. The same early-time behavior is observed for heavy quarks, see Fig. 4(a). Nevertheless, at later proper times, the jets do not appear to undergo multiple oscillations. This was noticed in the limiting cases shown in Fig. 3 and we have checked that it is still present in realistic heavy quark and jets simulations. The other components, δp 2 y and δp 2 x , show a steady monotonic increase at late times. We also plot the jet broadening coefficientq i as the time derivative of the momentum broadening in Fig. 7(b). Similar to the case of heavy quarks, there is a strong peak at very early stages with τ < 0.1 fm/c and a quick decay afterwards. Due to the decrease of δp 2 z at later times (τ 0.6 fm/c), the longitudinal componentq z becomes negative.
Results for jets with various values of p x /m are shown in Fig. 8, where we also plot the momentum broadening anisotropy. The values of the initial jet momentum are chosen such that p x > 5 GeV. We also include the results for lightlike jets. As expected, one recovers the highly energetic jet limit by choosing a sufficiently large value for p x /m in the particle solver. Compared to heavy quarks, there is little difference in the results when accounting for finite masses and momenta. Increasing p x /m leads to a slight decrease in the longitudinal component δp 2 (at most 15%). The transverse component is affected in the opposite way: δp 2 y increases with p x /m (at most 25%). Remarkably, while the momenta are not strongly affected, the anisotropy is enhanced (up to 40-60%) for less relativistic jets with p x /m ≈ 1 as can be seen from the lower panel in Fig. 8. Figure 9 depicts jet trajectories overlaid on top of the initial energy density of the Glasma for various initial momenta p T ∈ {10, 20, 50} GeV. Here, instead of fixing p x as the initial direction, we choose the direction of the initial transverse momentum randomly. Unlike slow heavy quarks (see Fig. 6), jets propagate on straight lines, due to their high initial momentum. Similar to Fig. 8, the initial value for p T only weakly affects the jet trajectories.

VIII. SUMMARY AND OUTLOOK
We have investigated the impact of the early stages of heavy-ion collisions, namely the Glasma, on hard probes such as heavy quarks and jets. To accomplish this, we approximate these hard probes as classical colored particles and simulate their dynamics using Wong's equations on top of the non-Abelian background field of the boostinvariant Glasma. This work can be understood as an extension of earlier studies on highly energetic jets [23,24] and heavy quarks [27][28][29][30][31][32][33] in the pre-equilibrium medium, which were limited in different ways. Simulations of jets in the Glasma were based on the ultrarelativistic limit, i.e. the jets were assumed to be lightlike. Thus, these simulations only apply to jets at extremely high ener-gies. On the other hand, studies of heavy quarks in the Glasma relied on using SU(2) [29][30][31][32] instead of SU(3) as the gauge group, which can only provide a qualitative picture. Thus, to improve upon these earlier studies, we have developed a fully non-perturbative simulation of classical particles with SU(3) color charges based on Wong's equations. The background field in which the charges are moving in is provided by classical real-time simulations of the Glasma. As such, we have realized a unified numerical setup where the effects of the Glasma on both heavy quarks and jets can be studied quantitatively.
To measure the impact of the Glasma on hard probes, we focused on the momentum broadening components δp 2 i (τ ) , which describe how much momentum is accumulated by heavy quarks and jets as they pass through the medium. This observable is particularly interesting, because it can be related to transport coefficients such as the heavy quark diffusion coefficient κ and the jet momentum broadening coefficientq. Additionally, we studied anisotropy ratios of different components of δp 2 i (τ ) . As a consistency check for our simulations, we have performed non-trivial numerical checks of our code by comparing to certain limiting cases, where the dynamics of the hard probes become trivial. These cases are heavy quarks with infinite mass (static quarks) and jets at very high energies (lightlike jets), where the particle trajectories are fixed and the eikonal approximation applies. Consequently, it is possible to compute momentum broadening directly from Wilson loops of the background field, which provides a benchmark result that our particle simulations must be able to reproduce. By taking these limits in the particle solver and performing extensive numerical checks, we have verified that our numerical solutions to Wong's equations are indeed consistent with the calculation from Wilson loops.
Going towards more realistic settings, we then considered the effects of finite mass and initial momentum of the hard probes. In particular, we performed simulations for beauty and charm quarks. In both cases, we notice deviations from the static quark limit. We found that there is strong initial acceleration at early times which results in a strongly time-dependent diffusion coefficient κ, with a characteristic peak at early times τ 0.1 fm/c and a subsequent quick decay. This behavior differs from the standard Langevin or Boltzmann approaches, in which the momentum broadening grows slowly, is generally smaller and does not exhibit a peak [32]. Following [30], it is of future interest to investigate the impact of such large broadening induced by the Glasma on observables such as elliptic flow or nuclear modification factors in both proton-nucleus or nucleus-nucleus collisions. Moreover, our calculations showed that beauty quarks, even though they are heavier, accumulate more momentum compared to charm quarks. This is due to their larger mass, which allows them to be formed slightly earlier in the evolution of the Glasma, where the color fields are particularly strong. Regardless of quark species, there is a sizable momentum broadening anisotropy with δp 2 L > δp 2 T , i.e. more accumulation along the beam axis compared to the transverse plane at early times. Curiously, this effect is reversed at late times for charm quarks, where δp 2 L < δp 2 T . Most remarkably, we observed that the longitudinal component δp 2 L oscillates as a function of time. It is possible that this effect could be traced back to the existence of plasmon modes in the Glasma. The plasmon modes are a collective feature of the Glasma color fields themselves. They could further be transmitted to the particles propagating in these fields, thus causing oscillations in their accumulated momenta. Moreover, plasmon frequency oscillations were already observed in a study involving a Yang-Mills plasma with large occupation numbers [25]. The emergence of such oscillations only in the longitudinal direction and not in the transverse plane is intriguing and requires further investigation. Thus, a possible extension of the current work would be to determine the Glasma plasmon frequency using methods similar to [65][66][67].
We have performed analogous calculations for jets with finite mass and finite initial momenta. Similar to heavy quarks and also confirming previous studies [24], we found that the jet momentum broadening coefficientq is highly peaked at early times τ 0.1 fm/c. There is a rapid increase of both longitudinal and transverse components at early times, and a sizable momentum broadening anisotropy at later times with δp 2 L > δp 2 T . For less relativistic jets with |p| ∼ m, this anisotropy is more pronounced compared to the ultrarelativistic limit. In contrast to heavy quarks, the effects of finite masses and initial momentum are quantitatively less important. Remarkably, there is a notable absence of oscillatory behavior in the longitudinal (beam axis) component. Instead, δp 2 L exhibits a single peak around τ ≈ 0.8 fm/c. It would be interesting if this behavior could also be understood in terms of the excitation spectrum of the Glasma.
Besides determining the origin of the oscillations of δp 2 L , there are multiple other ways to extend our current work. Concerning the Glasma itself, a possible extension is to consider more complicated initial conditions beyond the McLerran-Venugopalan model used here. In particular, it would interesting to see the effects of more realistic transverse structure (such as in the IP-Glasma model [68,69]) or hot spots [70][71][72]. Another extension, related to the longitudinal structure of the colliding nuclei, could be to go beyond the boost-invariant approximation and consider the full 3+1 dimensional structure of the Glasma, either due to finite extent along the beam axis [73][74][75][76] or due to the JIMWLK evolution [77][78][79]. Although generalizing our numerical setup to 3+1 dimensions is in principle trivial, a large amount of computational resources would be required to carry out such simulations. In practice, this generalization might still be possible through the weak field approximation [80], which exhibits significantly reduced computational costs compared to lattice simulations at the expense of neglecting non-perturbative effects.
Regarding the dynamics of the hard probes, an immediate improvement would be the inclusion of the color current generated by the color charges as they propagate through the Glasma. This would induce a back reaction of the hard particles onto the Glasma. It has already been demonstrated in [32] that including the color current of heavy quarks does not significantly modify momentum broadening, spectra, or nuclear modification factor at early times. However, one would expect the back reaction to be more significant for jets, in particular regarding (classical) gluon radiation and energy loss. Unfortunately, fast moving charged particles in lattice simulations are plagued by the numerical Cherenkov instability which is not tractable in the current setup without significant changes to the numerical scheme [81].
Another interesting aspect, unrelated to classical particle simulations, would be a more detailed study of large temporal and lightlike Wilson loops in the Glasma. Beyond just the lowest moments δp 2 , the Wilson loops encode information about the probability P (p ⊥ ) that a hard parton picks up transverse momentum p ⊥ during its evolution [57,58]. The Wilson loop formulation therefore allows for the extraction of the collision kernel for momentum broadening. Such a quantity was computed in the context of anisotropic plasmas within a kinetic theory approach [82] or using perturbative computations [13] or non-perturbative lattice techniques [83]. Computing the collision kernel in the Glasma, which is an anisotropic and out-of-equilibrium medium, is an exciting prospect.
Lastly, there are additional observables which describe the effect of the Glasma on heavy quarks and jets, namely two-particle correlations that may be significantly affected by the large momentum broadening. In principle, these are possible observables within the available setup, which could be extended by off-central collisions, more sophisticated nuclear models, and more realistic ways of initializing particles in our simulation. We plan to include such features in our code and study the angular correlations of quark-antiquark pairs and how they are affected by the early stages of heavy-ion collisions. In this part of the appendix we collect some derivations and technical details regarding Wong's equations.

Wong's equations in Milne coordinates
Here, we provide an explicit derivation of Wong's equations in the Milne frame. We start from the covariant form given by Eq. (23).
The coordinate vector of the Milne frame is x µ = (τ, x, y, η) with Milne proper time τ and longitudinal space-time rapidity η. The coordinate change from the laboratory to the Milne frame is described by The inverse transformations are t = τ cosh η and z = τ sinh η. The components of the metric are g µν = diag(1, −1, −1, −τ 2 ). Consequently, the only nonvanishing Christoffel symbols are The Christoffel symbols of the second kind are related to the first-kind Christoffel symbols through which in Milne coordinates read The Christoffel symbols are used to relate the covariant derivative along the worldline of a particle, denoted by D/dτ , to the usual derivative d/dτ . For the four-velocity u µ of a particle, this relationship is given by Note that τ denotes the proper time in the rest frame of the particle, which should not be confused with Milne proper time τ . The transformations of the four-velocity components are u τ = cosh η u t − sinh η u z along with u η = −(sinh η u t + cosh η u z )/τ . The inverse transformations are u t = cosh η u τ + sinh η τ u η and u z = sinh η u τ + cosh η τ u η . The next step is to express derivatives with respect to τ in terms of τ -derivatives. In particular, we use m d/dτ = p τ d/dτ . This allows us to write the τ -evolution of the particle coordinates from Eq. (23) as where the temporal component p τ is given by with p 2 T = (p x ) 2 +(p y ) 2 . The Milne proper time evolution of the momenta is given by where the covariant derivatives can be written as The last set of equations follows from Eq. (A4) and p µ = mu µ . Finally, we write the components of the field strength tensor in terms of color-electric and -magnetic fields. Using the relations the spatial components of the momentum equations become The temporal component is given by

Conservation of Casimirs by Wong's equations
Wong's equations from Eq. (23) preserve the classical Casimir values of the color charge Q a . The quadratic and cubic Casimirs are given by Taking the time derivative of the quadratic Casimir and using the anti-symmetry of the structure constants immediately yields Using the identity (see e.g. [84]) valid for the fundamental generators the cubic Casimir can be written as where Q = Q a T a . The conservation of q 3 then follows from where we have used Eq. (24) in the second line. The last line vanishes for any Q and A µ due to the cyclic property of the trace.
The generators are orthonormal in the sense that where the Dynkin index T R depends on the chosen representation R. For the fundamental and adjoint representations it is given by The representations of the algebra su(2) are uniquely labeled by the value of the quadratic Casimir with D R the dimension of the representation. For the su(3) algebra, representations are additionally labeled by the cubic Casimir where d abc denote the symmetric structure constants (see Tab. II). For SU (2), d abc can be taken as zero. The values of the quadratic and cubic Casimirs are (see [84]) Classical color charges Q a can be understood as the limit of high dimensional representations of su(N c ). By analogy, it is then possible to also define Casimir invariants similar to the finite dimensional case in Eqs. (B5) and (B6). Thus, the quadratic and cubic Casimirs of Q a are given by and similarly The color charges Q a are real-valued numbers and there are D A = N 2 c − 1 components for a given SU(N c ) group.
The Casimir invariants in Eqs. (B8) and (B9) may be viewed as constraints on the set of possible color charges. Evidently, the manifold of admissible color charge vectors depends on the number of colors. For example, in SU(2) only the quadratic Casimir invariant applies because d abc = 0. Thus, due to D A = 3, the classical color charges are three-dimensional vectors constrained to a 2-sphere with radius q 2 . The manifold of SU(2) color charges is therefore two-dimensional. For SU (3), the manifold becomes more complicated: Firstly, due to D A = 8 for SU (3), the quadratic invariant constrains the possible choices of color charges to a 7-sphere with radius q 2 . Secondly, the cubic Casimir invariant further constrains the color charge manifold to a six-dimensional submanifold of R 8 . It turns out that not all choices of q 2 and q 3 are admissible. In contrast to SU(2), where the color charge manifold exists for any value q 2 > 0, the SU(3) color charge manifold only exists for certain values of q 2 and q 3 . We address this at the end of Appendix B 2.

Integration measure and n-point functions
The generators of SU(N c ) satisfy the following trace relations Tr where the anomaly coefficient A R is given by Similarly to Eq. (B10), one may impose that the averages performed over classical color charge configurations must satisfy where, in the three-point function, the imaginary part was discarded since the classical color charges are real valued and are also symmetric under color indices, while f abc is anti-symmetric. The integration measure dQ used in the definition of the n-point functions is constrained by the Casimir invariants. For SU(2) it reads and similarly for SU(3) The normalization constant c R is chosen such that the color charge distributions are normalized to unity Once a normalization for the integration measure is chosen, the values for the classical Casimirs q 2 and q 3 are fixed according to Eqs. (B12). Contracting the twopoint function in Eq. (B12b) with δ ab immediately yields In an analogous manner, contracting the three-point function of Eq. (B12c) with d abc yields which is analogous to Eq. (B16). Thus, according to Eqs. (B2) and (B7), the classical Casimirs are given by Similar relationships for the Casimirs of the classical color charges were established in [20][21][22]53]. In [85,86] within a transport theory, it is argued that a classical description of the color charges holds only for large representations. From this perspective, it is not suitable to assign Casimirs only for a single classical quark or gluon, but rather to work with higher-order representations and introduce them for the whole ensemble. It is within this context that the choices in Eqs. (B16) and (B17) for the classical Casimirs q 2 and q 3 are the product of the dimension of the representation D R and the group-theoretical Casimirs C 2 and C 3 .
Different choices for the two-and three-point functions would yield other classical quadratic and cubic Casimirs. For example, an elegant choice would be to set q 2 = C 2 and q 3 = C 3 , i.e. match the classical and the grouptheoretical Casimirs directly. This would allows us to avoid the matching procedure of Section VI, where we showed that observables such as momentum broadening δp 2 must be divided by a factor of D R in order to reproduce a similar calculation in perturbative QCD. However, we found that this choice is generally not admissible. Using numerical solution methods, we were not able to find a single color charge vector Q a for SU(3) quarks that satisfies both q 2 = C 2 and q 3 = C 3 in the fundamental representation. On the other hand, we found possible solutions in the case of SU(3) gluons and for both SU(2) quarks and gluons. It appears that in the case of SU(3) quarks, the color charge manifold embedded in R 8 and constrained by the two Casimir invariants does not exist. On the other hand, the choices in Eqs. (B16) and (B17) are admissible in the sense that there are valid solution vectors Q a for both quarks and gluons in SU(2) and SU(3). We require these solution vectors because our color charge sampling method (see Section IV B) is based on randomly rotating initial color charge vectors. For SU(3) gluons, which reside in the adjoint representation, the initial color vector is fixed by q 2 (A) = 24 and q 3 (A) = 0 while SU(3) quarks in fundamental representation are labeled by q 2 (F ) = 4 and q 3 (F ) = 10/3. In our numerical simulations, we use the initial color vectors where we use the notation #» Q ≡ (Q 1 , . . . , Q 8 ). The case of SU(2) is much simpler, because only the quadratic Casimir invariant applies. Possible color charges Q a are vectors on 2-sphere with radius q 2 .
To this end, we need to perform the integral over matrix elements of U ∈ SU(N c ) matrices, as can be seen by using Eq. (54) to further write and similarly for the two-and three-point functions. This matching requires the evaluation of certain integrals over SU(N c ). In what follows, we calculate them symbolically 4 , in the fundamental representation, and quote the main steps in the derivation.

Revising some relevant SU(Nc) integrals
According to [88], integrals over U ∈ SU(N c ) of the type where U ij denotes the matrix elements of U (in the fundamental representation), may be evaluated from the generating function We are particularly interested in Z 1,1 , Z 2,2 and Z 3,3 and their derivatives with respect to J and K. For each of these cases, let us see how we may generate our integrals of interest and how to evaluate them, using already computed expressions from [88][89][90], for N c = 3.
(C13) c. Integral involving three pairs of conjugate matrix elements Taking p = n = 3 gives the generating function Similarly, we want the LHS of Eq. (56b) to be expressible as in Eq. (42b). Besides a factor of 1/T 2 F from Eq. (54), the RHS is proportional to in which the integral involving SU(N c ) matrices may be evaluated from Eq. (C13). Inserting this relation back into Eq. (C24) yields where in the last step Eqs. (B10a) and (B10b) were used. This finally gives the RHS of Eq. (56b) as Using the classical Casimir of Eq. (B8) gives Q a 0 Q a 0 = D F C 2 (F ) = (N 2 c − 1)T F , which finally yields the desired two-point function exactly the two-point function from Eq. (42b).

c. Matching the three-point functions
The LHS of Eq. (56c) should be given by Eq. (42c), while the RHS contains the term dU Tr T a U T a U † Tr multiplied by 1/T 3 F = 8 coming from Eq. (54). When the integral over SU(N c ) is evaluated according to Eq. (C21) and using Eq. (B10a) along with Eq. (B10c), one may show that the only non-vanishing terms are given by (C29) where we used the fact that d abc is symmetric with respect to a, b, c, whereas f abc is anti-symmetric. This then yields Our publicly available simulation code [63] was subject to various numerical checks. All of them were done for very massive quarks with vanishing transverse momentum, but hold irrespective of particle type or initial momentum.
Firstly, we verified that the Casimirs of the classical color charges are perfectly conserved throughout the evolution, up to numerical precision. This requirement is automatically satisfied by the numerical method used to solve the evolution of the color charge, namely the color rotation with Wilson lines constructed on the lattice, as in Eq. (38).
Secondly, we checked that the Casimir scaling expressed in Eq. (81) is satisfied at all proper times. For this purpose, we initialized classical color charges in both fundamental and adjoint representations, for SU (2) and SU(3), by appropriately choosing the corresponding Casimir invariants from Eq. (44). We emphasize that classical momentum broadening components are mapped to quantum ones according to Eq. (77). Using  Additionally, we compared SU(2) to SU(3) momentum broadening components, see Fig. 10 (lower panel ). We find that, in the very early stages, there is no simple way to scale the results of SU(2) to SU (3). At later proper times, when the Glasma enters the free-streaming regime and the fields are dilute, there exists a simple scaling from SU(2) to SU(3) momentum broadening, but the scaling factor differs for adjoint and fundamental representations.

Temporal constraint
Numerically, the temporal component of the momentum is computed from Eq. (A6) but p τ may also be extracted using Eq. (32). We denote by δp τ the difference between these two distinct ways of computing the temporal component and represent it in Fig. 11. The difference is on the percent level and can be understood as an artifact of the lattice discretization.