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

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Simulating jets and heavy quarks in the glasma using the colored particle-in-cell method © Published by the American Physical Society. Funded by SCOAP3. Published version Avramescu, Dana; Băran, Virgil; Greco, Vincenzo; Ipp, Andreas; Müller, David; Ruggieri, Marco Avramescu, D., Băran, V., Greco, V., Ipp, A., Müller, D., & Ruggieri, M. (2023). Simulating jets and heavy quarks in the glasma using the colored particle-in-cell method. Physical Review D, 107, Article 114021. https://doi.org/10.1103/PhysRevD.107.114021 2023 Simulating jets and heavy quarks in the glasma using the colored particle-in-cell method Dana Avramescu ,1,2,* Virgil Băran,3,†Vincenzo Greco ,4,5,‡Andreas Ipp ,6,§ David Müller ,6,∥and Marco Ruggieri 4,¶ 1Department of Physics, P.O. Box 35, 40014 University of Jyväskylä, Finland 2Helsinki Institute of Physics, P.O. Box 64, 00014 University of Helsinki, Finland 3Faculty of Physics, University of Bucharest, Atomiştilor 405, Măgurele, Romania 4Department of Physics and Astronomy, University of Catania, Via S. Sofia 64, I-95123 Catania, Italy 5INFN-Laboratori Nazionali del Sud, Via S. Sofia 62, I-95123 Catania, Italy 6Institute for Theoretical Physics, TU Wien, Wiedner Hauptstraße 8, A-1040 Vienna, Austria (Received 11 April 2023; accepted 15 May 2023; published 15 June 2023) 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 ultrarelativistic lightlike jets and obtain excellent quantitative agreement with previous studies. Going towards realistic masses and initial momenta, 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. DOI: 10.1103/PhysRevD.107.114021 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–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. The glasma is described using a wider framework entitled Color Glass Condensate (CGC) [5–7] which is formulated at the high-energy limit of quantum chromodynamics (QCD). The field equations for the color fields of the gluons are solved numerically using methods from lattice QCD [8–10]. To describe the properties of hard probes from high-energy nuclear collisions, numerous approaches based on perturbative QCD (pQCD) techniques [11–14], lattice computations [15–19], or non- Abelian Yang-Mills transport theories [20–22] have been used. These probes are produced immediately after the collision and may be affected by the entire evolution of the resulting Quark Gluon Plasma (QGP). 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 ultrarelativistic lightlike 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 *Corresponding author. [email protected] †[email protected] ‡[email protected] §[email protected]t ∥[email protected]c.at ¶[email protected]nict.it Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. Funded by SCOAP3. PHYSICAL REVIEW D 107, 114021 (2023) 2470-0010=2023=107(11)=114021(30) 114021-1 Published by the American Physical Society (assuming the heavy quarks to be infinitely massive and static) and emphasizes the emergence of plasmon mass induced oscillations. In another series [27–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–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. Nevertheless, 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 lightlike 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 and ˆ qfor 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 boostinvariant Yang-Mills equations. In Sec. III we present Wong’s equations. In Sec. 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 lightlike jets. In Sec. VI we show how to reconcile classical particle simulations with the calculation of momentum broadening within pQCD. The results obtained with our particle solver are showcased in Sec. VII. Finally, Sec. 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 Appendixes Athrough D. II. GLASMA IN A NUTSHELL Within the Color Glass Condensate framework [5–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 DμFμν ¼Jν;ð1Þ with Dμð…Þ≡∂μð…Þ−ig½Aμ;…denoting the gaugecovariant 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≡ðx0x3Þ=ffiffiffi 2 p. Their color currents are given by Jμ A;B ¼δμρA;Bðx∓; x⊥Þ;ð2Þ where ρA;B represent classical color charge densities and the subscripts Aand Bdenote the two colliding nuclei. The color charge densities are treated as stochastic variables whose DANA AVRAMESCU et al. PHYS. REV. D 107, 114021 (2023) 114021-2 statistics are determined by the probability functional W½ρ. We take it to be given by the McLerran-Venugopalan (MV) model [41–43]. It considers the color charges ρto follow Gaussian statistics, which are determined by the one- and two-point correlators hρaðx∓; x⊥ÞiA;B ¼0; hρaðx∓; x⊥Þρbðy∓; y⊥ÞiA;B ¼g2λA;Bðx∓Þδab ×δðx∓−y∓Þδð2Þð x⊥− y⊥Þ; ð3Þ where λA;B is the average color charge per unit volume. One may extract μ2 A;B ¼Rdxλ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 nonzero components of the gauge field are given by A covðx∓; x⊥Þ≡αA;Bðx∓; x⊥Þ where αA;B obeys a Poisson equation restricted to the transverse plane Δ⊥αA;Bðx∓; x⊥Þ¼−ρcov A;Bðx∓; x⊥Þ;ð4Þ in which Δ⊥is the transverse Laplace operator. The Poisson equation can be formally solved via Fourier transformation αA;Bðx∓; x⊥Þ¼Zd2 k⊥ ˜ ρcov A;Bðx∓;  k⊥Þ  k2 ⊥þλ2expð−i  k⊥· x⊥Þ;ð5Þ 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 Ai A;Bðx∓; x⊥Þ¼i gVðx∓; x⊥Þ∂iV†ðx∓; x⊥Þ;ð6Þ with the lightlike Wilson line V† A;Bðx∓; x⊥Þ¼PexpigZx∓ −∞ dy∓αA;Bðy∓; x⊥Þ;ð7Þ where Pð…Þdenotes the path-ordering operation. In the ultrarelativistic limit, the nuclei are contracted to infinitesimally thin sheets. This can be expressed via Jμ A;B ¼δμδðx∓ÞρA;Bðx⊥Þwhere the two-dimensional charge densities obey the correlator hρað x⊥Þρbð y⊥ÞiA;B ¼g2μ2 A;Bδabδð2Þð x⊥− y⊥Þ:ð8Þ The transverse gauge fields are given by Ai A;Bðx∓;x ⊥Þ¼θðx∓Þαi A;Bð x⊥Þ;ð9Þ in which θrepresents the Heaviside function and αi A;Bð x⊥Þ¼i gVA;Bð x⊥Þ∂iV† A;Bð x⊥Þ;ð10Þ involves a Wilson line depending on the transverse coordinate, obtainable as VA;Bð x⊥Þ¼limx∓→∞VA;Bðx∓; x⊥Þ. We now consider the classical collision problem DμFμν ¼Jν AþJν B;ð11Þ 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 ultrarelativistic limit, the total color current generated by the two nuclei Jμ¼Jμ AþJμ Bpossesses 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 τ¼ffiffiffiffiffiffiffiffiffiffiffiffiffi 2xþx− p;η¼1 2lnxþ x−;ð12Þ with proper time τand space-time rapidity η. By fixing the residual gauge freedom by imposing the temporal gauge condition Aτ¼0and 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] Aiðτ; x⊥Þτ¼0¼αi Að x⊥Þþαi Bð x⊥Þ; Aηðτ; x⊥Þτ¼0¼ig 2½αi Að x⊥Þ;αi Bð x⊥Þ;ð13Þ accompanied by ∂τAiðτ; x⊥Þτ¼0¼∂τAηðτ; x⊥Þτ¼0¼0:ð14Þ The conjugate momenta associated with the gauge fields are Pi¼τ∂τAi;P η¼1 τ∂τAη:ð15Þ The Yang-Mills action expressed in Milne coordinates, together with boost-invariance, yields the field equations SIMULATING JETS AND HEAVY QUARKS IN THE GLASMA …PHYS. REV. D 107, 114021 (2023) 114021-3 ∂τPi¼τDjFji −ig τhAη;DiAηi; ∂τPη¼1 τDiðDiAηÞ;ð16Þ along with the Gauss constraint DiPiþ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 Eq. (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 Land accompanied by periodic boundary conditions for the fields, is discretized in N2points 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 ˆ i, where ˆ iis the unit vector along xi,isgivenby Uˆ iðτ; x⊥Þ≈expigaAiτ; x⊥þa 2 ˆ i:ð17Þ Links in opposite directions can be expressed through the Hermitian operation U− ˆ iðτ; x⊥Þ≡U† ˆ iðτ; x⊥− ˆ iÞ.These gauge links are then used to construct a plaquette variable as Uˆ i ˆ jðτ; x⊥Þ≡Uˆ iðτ; x⊥ÞUˆ jðτ; x⊥þˆ iÞU− ˆ iðτ; x⊥þˆ iþˆ jÞ U− ˆ jðτ; x⊥þˆ jÞ. The Yang-Mills action can be approximated using link and plaquette variables along with the conjugate momenta Pηðτ; x⊥Þ¼1 τ∂τAηðτ; x⊥Þ; Piðτ; x⊥Þ¼−iτ ga h∂τUˆ iðτ; x⊥ÞiU† ˆ iðτ; x⊥Þ:ð18Þ Varying the discretized action yields discretized equations of motion ∂τPηðτ; x⊥Þ¼1 τ D2 iAηðτ; x⊥Þ; ∂τPiðτ; x⊥Þ¼−X j τ ga3hUˆ i ˆ jðτ; x⊥ÞþUˆ i− ˆ jðτ; x⊥Þiah −ig τhAtransp ηðτ; x⊥Þ;DF iAηðτ; x⊥Þi;ð19Þ where D2 i≡DF iDB icontains the forward DF iand backward DB igauge-covariant finite differences on the lattice, ð…Þah denotes the anti-Hermitian traceless part of a matrix, and Atransp ηðτ; x⊥Þ≡Uˆ iðτ; x⊥ÞAηðτ; x⊥þˆ iÞU† ˆ iðτ; x⊥Þ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 nontrivial path-ordering of the involved Wilson lines. Nevertheless, by sticking together infinitesimally thin sheets of color charge and regularizing the correlator as [46] hρa mð x⊥Þρb nð y⊥ÞiA;B ¼1 Nsa2g2μ2 A;Bδmnδabδð x⊥− y⊥Þ;ð20Þ where m; n ∈f1;2;…Nsgdenotes the index of the sheet, and with Nsthe 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 cutoffs are λand Λ. Furthermore, the Wilson lines are constructed as products computed for each sheet as V†ð x⊥Þ¼QNs n¼1expð−igαnð x⊥ÞÞ with αn≡αa nTa. 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 1The simulation code for the glasma fields is publicly available at https://gitlab.com/openpixi/curraun. DANA AVRAMESCU et al. PHYS. REV. D 107, 114021 (2023) 114021-4 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 dxi dt¼pi E; dpi dt¼gQaFiμ;a pμ E; dQa dt¼−gfabcAb μQcpμ E;ð21Þ where i∈fx; y; zgand a¼f1;2;…;D Agwith the dimension of the adjoint representation DA¼N2 c−1. The energy is given by E¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  p2þm2 pwith mbeing the mass of the particle,  p≡ðpx;p y;p zÞ, and fabc the structure constants for the SUðNcÞgroup. The dynamic equation for the energy p0¼Eis given by dE dt¼gQaF0i;a pi E¼gQa Ea· v; ð22Þ with Ei≡F0idenoting 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  vthe 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; Dpμ dτ¼gQaFμν;a pν m; dQa dτ¼−gfabcAb μQcpμ m;ð23Þ where ðdτÞ2¼gμνdxμdxνdenotes the relativistic proper time and in which Ed=dt¼md=dτis employed and D=dτ is the covariant derivative taken along the particle worldline. Further, we make use of the Lie-algebra-valued color charges Q¼QaTato write QaFμν;a ¼Tr½QFμν=TR, where TRis the representation-dependent Dynkin index defined through Tr½TaTb¼TRδab, with Ta∈SUðNcÞ. The above equations simplify to Dpμ dτ¼g TR Tr½QFμνpν m; dQ dτ¼−ig½Aμ;Qpμ m:ð24Þ The equation governing the evolution of the color charge may formally be solved by QðτÞ¼Uðτ;τ0ÞQðτ0ÞUðτ0;τÞ;ð25Þ where Qis rotated with Wilson lines. These Wilson lines involve the path-ordered exponential computed along the trajectory of the particle and are given by Uðτ;τ0Þ¼Pexp0 @−igZτ τ0 dτ0dxμ dτ0AμðxμÞ1 A:ð26Þ The last relation may be derived by making use of the parallel transport equation for a Wilson line d dτ Uðτ;τ0Þ¼−igdxμ dτ AμðxμðτÞÞUðτ;τ0Þ:ð27Þ The use of Wilson lines in the evolution of the color charge automatically conserves the quadratic QaQa≡q2ðRÞð28Þ and cubic classical Casimirs2 dabcQaQbQc≡q3ðRÞ;ð29Þ where dabc are the symmetric structure constants and the values q2ðRÞ,q3ðRÞdepend on the chosen representation R for the color charge Q. We go into detail about how these invariants are fixed in Sec. 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 A1. The coordinate Wong equations are given by dxμ dτ¼pμ pτð30Þ where xμ∈fx; y; ηgand ðpτÞ2¼ðpxÞ2þðpyÞ2þ τ2ðpηÞ2þm2. The Wong equations for momenta read as 2One may also define Casimir invariants for classical Lie algebras. For SU(2) and SU(3), we construct them in analogy with the group-theoretical Casimir invariants. More details are offered in Appendix B. SIMULATING JETS AND HEAVY QUARKS IN THE GLASMA …PHYS. REV. D 107, 114021 (2023) 114021-5 τdpη dτþ2pη¼g TRTr½QEη−Tr½QBxpy pτþTr½QBypx pτ; dpx dτ¼g TRTr½QExþTr½QBηpy pτ−Tr½QByτpη pτ; dpy dτ¼g TRTr½QEy−Tr½QBηpx pτþTr½QBxτpη pτ; ð31Þ and are accompanied by the dynamic equation for the temporal component dpτ dτþτpη pτpη¼g TRTr½QEητpη pτþTr½QExpx pτ þTr½QEypy pτ;ð32Þ where the color-electric and -magnetic fields are determined from the field strength tensor via Ei≡Fτi;B i≡ϵij 1 τFηj; Eη≡1 τFτη;B η≡−Fxy:ð33Þ 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 Uðτ;τ0Þ¼Pexp−igZxμðτÞ xμðτ0Þ dxμAμðxμðτÞÞ:ð34Þ As used in the glasma framework, we employ the temporal gauge Aτ¼0and 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 δp2 μðτÞ≡p2 μðτÞ−p2 μðτformÞ:ð35Þ 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 pT¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðpxÞ2þðpyÞ2 pat formation time, with fixed pTðτformÞand pηðτformÞ¼0for heavy quarks, or an initial pxðτformÞand py;ηðτformÞ¼0for 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 Sec. 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 Eiðτn;xnÞ¼1 τn Piðτn;xnÞ≈ 1 4τnPi xnτnþΔτ 2þPi xnτn− Δτ 2 þUxn;− ˆ iðτnÞPi xn− ˆ iτnþΔτ 2þPi xn− ˆ iτn− Δτ 2U† xn;− ˆ iðτnÞ; Eηðτn;xnÞ¼Pηðτn;xnÞ≈ 1 2Pη xnτnþΔτ 2þPη xnτn− Δτ 2;ð36Þ where i¼x,y, and similarly DANA AVRAMESCU et al. PHYS. REV. D 107, 114021 (2023) 114021-6 Biðτn;xnÞ¼− 1 τn DiAηðτn;xnÞ≈− 1 2τnaThUxn;ˆ iðτnÞAxnþˆ i;ηðτnÞU† xn; ˆ iðτnÞ−Uxn;− ˆ iðτnÞAx− ˆ i;ηðτnÞU† xn;− ˆ iðτnÞi; Bηðτn;xnÞ¼−Fxyðτn;xnÞ≈− 1 4ga2 ThUxn;ˆxˆyðτnÞþUxn;ˆy−ˆxðτnÞþUxn;−ˆx−ˆyðτnÞþUxn;−ˆyˆxðτnÞiah:ð37Þ 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 xn≡½ x⊥ðτnÞNGP. 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-in-cell (CPIC) methods used in the context of particles in classical Yang-Mills (CYM) plasmas [49–52].Inthe 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 Yang-Mills (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 τfas being comprised of subsequent products of “short”Wilson lines Uðτn−1;τnÞas Uðτi;τfÞ≈Uðτi;τiþ1ÞUðτiþ1;τiþ2Þ…Uðτf−1;τfÞ.These short Wilson lines may be reduced to Uðτn−1;τnÞ≃exp0 @igZxn xn−1 dx0iAiðx0Þ1 A ×expðigδηnAηðxnÞÞ ¼Uxn−1;ˆ iðτnÞUxn;ˆηðτnÞ:ð38Þ Here, Uxn;ˆ iðτnÞis a transverse gauge link along the direction ˆ iwith i¼x,yevaluated at position xn,while Uxn;ˆ ηðτ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 δηnis numerically small and it follows that ½RdxiAi;δηnAη≃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 dQ dτ¼ig½Q; Axpx pτþ½Q; Aypy pτþ½Q; Aηpη pτ;ð39Þ FIG. 1. Trajectories of charm quarks in glasma, simulated with our particle solver, where the proper time evolution of (a) positions and (b) momenta are given by Eqs. (30) and (31) for Ntp ¼100 test particles, all initialized with zero px;y;zðτformÞ. The color represents the value of the proper time difference δτ ≡τ−τform at which the coordinates or momenta are evaluated. SIMULATING JETS AND HEAVY QUARKS IN THE GLASMA …PHYS. REV. D 107, 114021 (2023) 114021-7 where the transverse gauge fields are numerically extracted from the gauge links using matrix logarithms igaAxxþa 2;y ¼lnðUˆ xðx; yÞÞ; igaAyx; y þa 2¼lnðUˆyðx; yÞÞ:ð40Þ 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ðNcÞ, 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 Qin the ensemble of partons. Here, we largely follow the seminal works on classical non-Abelian transport theory [20–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 Qa as stochastic variables with fixed values of q2¼QaQaand q3¼dabcQaQbQc. We emphasize that this is a choice and in our framework, where the Casimirs remain constant throughout the evolution, see Sec. A2, the obvious choice is to fix the Casimirs. In analogy with the trace relations for operator-valued elements of the suðNcÞcolor algebra Tr½ˆ Qa¼0; Tr½ˆ Qaˆ Qb¼TRδab; Tr½ˆ Qaˆ Qbˆ Qc¼AR 4ðdabc þifabcÞ;ð41Þ we choose to have color charges randomly distributed according to one-, two- and three-point functions hQai¼0;ð42aÞ hQaQbi¼TRδab;ð42bÞ hQaQbQci¼AR 4dabc;ð42cÞ where the representation-dependent coefficients TRand AR are given by TR¼1 2;R¼F Nc;R¼A;A R¼1;R¼F 0;R¼A:ð43Þ The Casimir invariants from Eqs. (28) and (29) constrain what values the color charges Qacan 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 q2ðRÞ¼N2 c−1 2;R¼F NcðN2 c−1Þ;R¼A ;ð44aÞ q3ðRÞ¼ðN2 c−4ÞðN2 c−1Þ 4Nc;R¼F 0;R¼A :ð44bÞ 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 FIG. 2. Diagram with the color rotation performed during a numerical time step from τn−1to τ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 Uˆx. 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Þ¼UˆxðτnÞUˆηðτnÞ. DANA AVRAMESCU et al. PHYS. REV. D 107, 114021 (2023) 114021-8 without expansion. In these systems, the accumulated momenta of heavy quarks exhibit damped oscillations with the plasmon frequency. It is likely that the longitudinal component hδp2 ziin the glasma oscillates for a similar reason (plasmon excitations), although it is not clear why only the longitudinal component is affected. 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 κinst L;TðτÞ≡d dτhδp2 L;TðτÞi:ð84Þ 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 hδp2 Liincreases more rapidly than the transverse component hδp2 Tiat early times. Even though not shown here, we checked that the longitudinal and transverse momentum broadenings have the same behaviors at larger proper times τ≫2fm=c, as already noticed in Fig. 3for static quarks. The first peak of the oscillations in hδp2 Lihappens at around δτ ¼0.8fm=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 pT. Figure 5shows the numerical results for beauty and charm quarks for various values of the initial pT∈f0;2;5;10g GeV.3For 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 as heavy quark anisotropy ≡hδp2 Li hδp2 Ti:ð85Þ 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]. FIG. 4. (a) Longitudinal and transverse momentum broadening components of beauty quarks formed at τform ≈0.02 fm=c, initialized with pTðτformÞ¼0GeV, as a function of the time difference δτ ≡τ−τform. (b) Derivatives of the accumulated momenta which give the transport coefficients according to Eq. (84). 3It should be noted that at the highest initial transverse momenta, these heavy quarks essentially behave like jets. SIMULATING JETS AND HEAVY QUARKS IN THE GLASMA …PHYS. REV. D 107, 114021 (2023) 114021-15 Focusing on the heavy quark anisotropy, we find that as the initial pTincreases, hδp2 Lidecreases and hδp2 Ti increases. Consequently, the corresponding anisotropy hδp2 Li=hδp2 Tibecomes smaller. Compared to the static quark accumulated momentum (dashed lines), beauty quarks with zero initial pThave an increase in hδp2 Liof 50% and charm quarks 50–80% throughout the proper time evolution. For the maximum initial pTtaken in our simulations, hδp2 Lifor dynamic quarks differs from that for static quarks by 50% for beauty and 30–70% for charm quarks, whereas hδp2 Tiincreases only by 20–30% compared to static quarks. These will have a complementary effect on the anisotropy. Namely, hδp2 Li=hδp2 Tifor beauty or charm quarks is 20–40% larger or smaller, depending on the initial pT, than that of infinitely massive heavy quarks formed at the same formation time. The anisotropy is higher for small initial pTheavy quarks and lower for quite large initial pT. The anisotropy is more pronounced for the zero initial pTheavy 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 τ¼0fm=c, the glasma consists of correlated domains of longitudinal color-electric and - magnetic flux tubes with a typical size of ≈1=Q2 s. This shortly lived initial phase is probed by heavy quarks with very high mass mHQ due to their early formation time τform ¼1=ð2mHQÞ. 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. 4and 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. 5for 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=Qs≈0.1fm=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 FIG. 5. (Top) Longitudinal and transverse momentum broadening components, along with their ratio (bottom). The simulations are performed for (left) beauty and (right) charm quarks for various values of initial transverse momentum (colored full lines). We compare to the static case (gray dashed line), when the quarks are considered infinitely massive and the accumulated momentum is extracted solely from color-electric correlator, see Eq. (64). DANA AVRAMESCU et al. PHYS. REV. D 107, 114021 (2023) 114021-16 κ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 pT∈f0;2;5gGeV 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. 6describes an oversimplified 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 lightlike 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 pxvalues. Similarly to the heavy quark transport coefficient κ,we distinguish between various components of the jet transport coefficient ˆ q. We define the instantaneous jet broadening coefficient ˆ qiðτÞ≡d dτhδp2 iðτÞi ð86Þ with i∈fx; y; zg. This is different from the collisional energy loss dE=dx. Since the jet propagates along the xaxis, we introduce the transverse ˆ qT≡ˆ qyand longitudinal ˆ qL≡ˆ qzjet transport coefficients. In addition, we are interested in deviations from the lightlike jet scenario by considering a finite jet mass mand an initial px. We find that jet momentum broadening essentially only depends on the ratio px=m. In the limit FIG. 6. (Colored lines) trajectories of heavy quarks propagating in a single glasma flux tube evolved up to τ¼0.2fm=c. All partons are produced at the center of a flux tube, where the energy density was locally maximal at the creation time of the glasma. We consider three cases: (left) very massive quarks with m¼200 GeV (approaching the static quark limit) with τform ¼0fm=c,(middle) beauty quarks with τform ¼0.02 fm=c and (right) charm quarks with τform ¼0.06 fm=c. The initial transverse momentum pTis varied between 0 GeV and 5 GeV. The background shows the energy density at formation time of the respective particle type, namely ϵform ≡ϵðτformÞ. SIMULATING JETS AND HEAVY QUARKS IN THE GLASMA …PHYS. REV. D 107, 114021 (2023) 114021-17 of px=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 jet anisotropy ≡hδp2 zi hδp2 yi;ð87Þ Our numerical results for jets are shown in Figs. 7and 8. Figure 7(a) shows the accumulated momentum broadening for a quark jet with m¼1GeV and initial px¼10 GeV as a function of Milne proper time τ. The longitudinal component hδp2 zi(along the beam axis) shows similar behavior as in the case of heavy quarks. After reaching a maximum at roughly τ≈0.8fm=c, the longitudinal component hδp2 zistarts to decrease at late times τ≳1fm=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. 3and we have checked that it is still present in realistic heavy quark and jets simulations. The other components, hδp2 yiand hδp2 xi, show a steady monotonic increase at late times. We also plot the jet broadening coefficient ˆ qias 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.1fm=c and a quick decay afterwards. Due to the decrease of hδp2 ziat later times (τ≳0.6fm=c), the longitudinal component ˆ qzbecomes negative. Results for jets with various values of px=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 px>5GeV. We also include the results for lightlike jets. As expected, one recovers the highly energetic jet limit by choosing a sufficiently large value for px=m in the particle solver. Compared to heavy quarks, there is little difference in the results when accounting for finite masses and momenta. Increasing px=m leads to a slight decrease in the longitudinal component hδp2 zi(at most 15%). The transverse component is affected in the FIG. 7. (a) Momentum broadening components of jets with m¼1GeV and initial px¼10 GeV, along the x,y,z-axes, as a function of proper time and (b) the derivative of the accumulated momenta that produces components of the jet transport coefficients according to Eq. (86).(Insets) Zoom-in on the very early stage. FIG. 8. Momentum broadening along (top) z-axis and y-axis, together with (bottom) their ratio, which is a measure of the momentum-broadening anisotropy. The simulations are performed for various values of px=m ∈f1;2;5;10g(colored full lines), compared to lightlike jets moving along the x-axis (gray dashed line). For large px=m, the jet becomes lightlike and our particle simulations approach the limiting case. DANA AVRAMESCU et al. PHYS. REV. D 107, 114021 (2023) 114021-18 opposite way; hδp2 yiincreases with px=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 px=m ≈1as can be seen from the lower panel in Fig. 8. Figure 9depicts jet trajectories overlaid on top of the initial energy density of the glasma for various initial momenta pT∈f10;20;50gGeV. Here, instead of fixing pxas 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 pTonly 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 boost-invariant glasma. This work can be understood as an extension of earlier studies on highly energetic jets [23,24] and heavy quarks [27–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 energies. On the other hand, studies of heavy quarks in the glasma relied on using SU(2) [29–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 nonperturbative 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 hδp2 iðτÞ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 coefficient ˆ q. Additionally, we studied anisotropy ratios of different components of hδp2 iðτÞi. As a consistency check for our simulations, we have performed nontrivial 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 κ,witha characteristic peak at early times τ≲0.1fm=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],itisof 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 hδp2 Li>hδp2 Ti, i.e. more accumulation along the beam axis compared to FIG. 9. (Colored lines) Trajectories of jets propagating out of a single glasma flux tube evolved up to τsim ¼0.2fm=c. The colors of the lines indicate the initial momentum pT∈f10;20;50gGeV. All jets are initialized with m¼1GeV. The jet trajectories are essentially straight and are barely affected by the color fields of the glasma. SIMULATING JETS AND HEAVY QUARKS IN THE GLASMA …PHYS. REV. D 107, 114021 (2023) 114021-19 the transverse plane at early times. Curiously, this effect is reversed at late times for charm quarks, where hδp2 Li<hδp2 Ti. Most remarkably, we observed that the longitudinal component hδp2 Lioscillates 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–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],wefound that the jet momentum broadening coefficient ˆ qis highly peaked at early times τ≲0.1fm=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 hδp2 Li>hδp2 Ti. For less relativistic jets with jpj∼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, hδp2 Liexhibits a single peak around τ≈0.8fm=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 hδp2 Li, 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–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–76] or due to the JIMWLK evolution [77–79]. Although generalizing our numerical setup to 3þ1dimensions 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 nonperturbative 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 hδp2i, 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 nonperturbative 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. ACKNOWLEDGMENTS D. A. acknowledges funding from the Academy of Finland, Center of Excellence in Quark Matter Project No. 346324. V. G. acknowledges funding from UniCT under “Linea di intervento 2”(HQCDyn Grant). D. M. acknowledges funding from the Austrian Science Fund (FWF) Projects No. P 34455 and No. P 34764. All simulations were performed using the GPU nodes of the Center of Theoretical Physics, University of Bucharest. D. A., D. M. and M. R. acknowledge S. Mrówczyński and C. Manuel for discussions regarding classical color charges, and K. Boguslavksi, H. Mäntysaari, and T. Lappi for many insightful discussions regarding the early stages of heavy-ion collisions. M. R. acknowledges John Petrucci for inspiration. We are grateful to T. Lappi for reviewing the manuscript. DANA AVRAMESCU et al. PHYS. REV. D 107, 114021 (2023) 114021-20 APPENDIX A: SOME DETAILS REGARDING WONG’S EQUATIONS In this part of the appendix we collect some derivations and technical details regarding Wong’s equations. 1. 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 τ¼ffiffiffiffiffiffiffiffiffiffiffiffiffi t2−z2 p;η¼1 2lntþz t−z:ðA1Þ 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 Γτ ηη ¼τ;Γη τη ¼Γη ητ ¼1 τ:ðA2Þ The Christoffel symbols of the second kind are related to the first-kind Christoffel symbols through ½ab; c¼gcdΓd ab; which in Milne coordinates read ½ηη;τ¼τ;½ητ;η¼½τη;η¼−τ:ðA3Þ 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 Duμ dτ¼gμν duν dτþ½νλ;μuνuλ:ðA4Þ 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 ηut−sinh ηuzalong with uη¼−ðsinh ηutþcosh ηuzÞ=τ. The inverse transformations are ut¼cosh ηuτþsinh ητuηand uz¼sinh ηuτþ cosh ητuη. The next step is to express derivatives with respect to τin terms of τ-derivatives. In particular, we use md=dτ¼ pτd=dτ. This allows us to write the τ-evolution of the particle coordinates from Eq. (23) as dx dτ¼px pτ;dy dτ¼py pτ;dη dτ¼pη pτ;ðA5Þ where the temporal component pτis given by pτ¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p2 Tþτ2ðpηÞ2þm2 q;ðA6Þ with p2 T¼ðpxÞ2þðpyÞ2. The Milne proper time evolution of the momenta is given by Dpν dτ¼g TR Tr½QFνμpμ pτ;ðA7Þ where the covariant derivatives can be written as Dpτ dτ¼dpτ dτþτðpηÞ2 pτ;ðA8Þ Dpi dτ¼−dpi dτ;ðA9Þ Dpη dτ¼−τ2dpη dτ−2τpη:ðA10Þ 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 Ei≡Fτi;B i≡ϵij 1 τFηj; Eη≡1 τFτη;B η≡−Fxy;ðA11Þ the spatial components of the momentum equations become τdpη dτþ2pη¼g TRTr½QEη−Tr½QBxpy pτþTr½QBypx pτ; dpx dτ¼g TRTr½QExþTr½QBηpy pτ−Tr½QByτpη pτ; dpy dτ¼g TRTr½QEy−Tr½QBηpx pτþTr½QBxτpη pτ: ðA12Þ The temporal component is given by dpτ dτþτpη pτpη¼g TRTr½QEητpη pτþTr½QExpx pτ þTr½QEypy pτ:ðA13Þ SIMULATING JETS AND HEAVY QUARKS IN THE GLASMA …PHYS. REV. D 107, 114021 (2023) 114021-21 2. Conservation of Casimirs by Wong’s equations Wong’s equations from Eq. (23) preserve the classical Casimir values of the color charge Qa. The quadratic and cubic Casimirs are given by QaQa¼q2;ðA14Þ dabcQaQbQc¼q3:ðA15Þ Taking the time derivative of the quadratic Casimir and using the antisymmetry of the structure constants immediately yields dq2 dτ¼2QadQa dτ¼−2gpμ mAb μfabcQaQc¼0:ðA16Þ Using the identity (see e.g. [84]) valid for the fundamental generators Tr½TaTbTc¼1 4ðdabc þifabcÞ;ðA17Þ the cubic Casimir can be written as q3¼4Tr½Q3;ðA18Þ where Q¼QaTa. The conservation of q3then follows from d dτ Tr½Q3¼3TrQ2dQ dτ ¼−3igTr½Q2½Aμ;Qpμ m ¼0;ðA19Þ where we have used Eq. (24) in the second line. The last line vanishes for any Qand Aμdue to the cyclic property of the trace. APPENDIX B: PROPERTIES OF CLASSICAL COLOR CHARGES Here we provide some additional mathematical details regarding classical color charges such as the choice of Casimir invariants and the integration measure. 1. Casimir invariants for generators Taand classical color charges Qa Let fTagwith a∈f1;2;…DAgbe the generators of SUðNcÞwith DA¼N2 c−1. Common choices for generators in the fundamental (quark) representation are Ta¼ σa=2for SU(2) and Ta¼λa=2for SU(3), with Pauli matrices σaand Gell-Mann matrices λa. Regardless of a particular representation, the generators satisfy ½Ta;Tb¼ ifabcTc, where fabc are the totally antisymmetric structure constants of the group. For SU(2) these constants are given by fabc ¼ϵabc. For SU(3) the nonvanishing structure constants are listed in Table I. The other representation that we are interested in is the adjoint (gluon) representation, whose generators are given by ðTaÞbc ¼−ifabc:ðB1Þ The dimensions of the fundamental (F) and adjoint (A) representations are DR¼Nc;R¼F N2 c−1;R¼A:ðB2Þ The generators are orthonormal in the sense that Tr½TaTbR¼TRδab;ðB3Þ where the Dynkin index TRdepends on the chosen representation R. For the fundamental and adjoint representations it is given by TR¼1 2;R¼F Nc;R¼A:ðB4Þ The representations of the algebra suð2Þare uniquely labeled by the value of the quadratic Casimir X a TaTa¼1DRC2ðRÞ;ðB5Þ with DRthe dimension of the representation. For the ðsuÞ3 algebra, representations are additionally labeled by the cubic Casimir X abc dabcTaTbTc¼1DRC3ðRÞ;ðB6Þ where dabc denote the symmetric structure constants (see Table II). For SU(2), dabc can be taken as zero. The values of the quadratic and cubic Casimirs are (see [84]) C2ðRÞ¼(N2 c−1 2Nc;R¼F Nc;R¼A ;ðB7aÞ TABLE I. Antisymmetric structure constants for SU(3). fabc f123 f147 f156 f246 f257 f345 f367 f458 f678 11 2−1 2 1 2 1 2 1 2−1 2ffiffi3 p 2ffiffi3 p 2 DANA AVRAMESCU et al. PHYS. REV. D 107, 114021 (2023) 114021-22 C3ðRÞ¼(ðN2 c−4ÞðN2 c−1Þ 4N2 c;R¼F 0;R¼A :ðB7bÞ Classical color charges Qacan be understood as the limit of high-dimensional representations of suðNcÞ. 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 Qaare given by X a QaQa≡q2ðRÞ;ðB8Þ and similarly X abc dabcQaQbQc≡q3ðRÞ:ðB9Þ The color charges Qaare real-valued numbers and there are DA¼N2 c−1components for a given SUðNcÞ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 dabc ¼0. Thus, due to DA¼3, the classical color charges are three-dimensional vectors constrained to a 2-sphere with radius q2. The manifold of SU(2) color charges is therefore two-dimensional. For SU(3), the manifold becomes more complicated: Firstly, due to DA¼8for SU(3), the quadratic invariant constrains the possible choices of color charges to a 7-sphere with radius q2. Secondly, the cubic Casimir invariant further constrains the color charge manifold to a six-dimensional submanifold of R8. It turns out that not all choices of q2and q3are admissible. In contrast to SU(2), where the color charge manifold exists for any value q2>0, the SU(3) color charge manifold only exists for certain values of q2and q3. We address this at the end of Appendix B2. 2. Integration measure and n-point functions The generators of SUðNcÞsatisfy the following trace relations Tr½Ta¼0;ðB10aÞ Tr½TaTb¼TRδab;ðB10bÞ Tr½TaTbTc¼AR 4ðdabc þifabcÞ;ðB10cÞ where the anomaly coefficient ARis given by AR¼1;R¼F 0;R¼A:ðB11Þ Similarly to Eq. (B10), one may impose that the averages performed over classical color charge configurations must satisfy hQai≡ZdQQa¼0;ðB12aÞ hQaQbi≡ZdQQaQb¼TRδab;ðB12bÞ hQaQbQci≡ZdQQaQbQc¼AR 4dabc;ðB12cÞ 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 fabc is antisymmetric. The integration measure dQ used in the definition of the n-point functions is constrained by the Casimir invariants. For SU(2) it reads dQ¼cRd3QδðQaQa−q2Þ;ðB13Þ and similarly for SU(3) dQ¼cRd8QδðQaQa−q2ÞδðdabcQaQbQc−q3Þ:ðB14Þ The normalization constant cRis chosen such that the color charge distributions are normalized to unity ZdQ¼1:ðB15Þ Once a normalization for the integration measure is chosen, the values for the classical Casimirs q2and q3are fixed according to Eqs. (B12). Contracting the two-point function in Eq. (B12b) with δab immediately yields q2¼TRDA, which, by virtue of TRDA¼DRC2ðRÞ, yields q2ðRÞ¼DRC2ðRÞ:ðB16Þ TABLE II. Symmetric structure constants for SU(3). dabc d118 d146 d157 d228 d247 d256 d338 d344 d355 d366 d377 d448 d558 d668 d778 d888 1ffiffi3 p1 2 1 2 1ffiffi3 p−1 2 1 2 1ffiffi3 p1 2 1 2−1 2−1 2−1 2ffiffi3 p−1 2ffiffi3 p−1 2ffiffi3 p−1 2ffiffi3 p−1ffiffi3 p SIMULATING JETS AND HEAVY QUARKS IN THE GLASMA …PHYS. REV. D 107, 114021 (2023) 114021-23 In an analogous manner, contracting the three-point function of Eq. (B12c) with dabc yields q3ðRÞ¼DRC3ðRÞ;ðB17Þ which is analogous to Eq. (B16). Thus, according to Eqs. (B2) and (B7), the classical Casimirs are given by q2ðRÞ¼(N2 c−1 2;R¼F NcðN2 c−1Þ;R¼A ;ðB18aÞ q3ðRÞ¼(ðN2 c−4ÞðN2 c−1Þ 4Nc;R¼F 0;R¼A :ðB18bÞ Similar relationships for the Casimirs of the classical color charges were established in [20–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 q2and q3 are the product of the dimension of the representation DR and the group-theoretical Casimirs C2and C3. 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 q2¼C2 and q3¼C3, i.e. match the classical and the grouptheoretical Casimirs directly. This would allows us to avoid the matching procedure of Sec. VI, where we showed that observables such as momentum broadening hδp2imust be divided by a factor of DRin 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 Qafor SU(3) quarks that satisfies both q2¼ C2and q3¼C3in 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 R8and 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 Qafor both quarks and gluons in SU(2) and SU(3). We require these solution vectors because our color charge sampling method (see Sec. 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 q2ðAÞ¼24 and q3ðAÞ¼0while SU(3) quarks in fundamental representation are labeled by q2ðFÞ¼4and q3ðFÞ¼10=3. In our numerical simulations, we use the initial color vectors  QðAÞ¼ð4.89898;0;0;0;0;0;0;0Þ;  QðFÞ¼ð0;0;0;0;−1.69469;0;0;−1.06209Þ;ðB19Þ where we use the notation  Q≡ðQ1;…;Q 8Þ. The case of SU(2) is much simpler, because only the quadratic Casimir invariant applies. Possible color charges Qaare vectors on 2-sphere with radius q2. APPENDIX C: SAMPLING CLASSICAL COLOR CHARGES VIA THE HAAR MEASURE This appendix contains an analytical proof for the matching from Eq. (56) of the one-, two- and three-point functions of the classical color charges computed via the Haar measure, as given in Eq. (55a), to the desired values provided in Eq. (42). For example, the aim is to derive that the one-point function extracted by integrating over the manifold of SUðNcÞ hQaiU¼Qa0 0ZdUUaa0 where Qa0 0is the initial color vector fixed by the quadratic and cubic Casimirs from Eqs. (44a) and (44b), does indeed give the expected value hQaiU↦hQai¼0: To this end, we need to perform the integral over matrix elements of U∈SUðNcÞmatrices, as can be seen by using Eq. (54) to further write ZdUUaa0¼ZdU1 TR Tr½TaUTa0U† ¼1 TR Ta liTa0 jk ZdUUijU† kl and similarly for the two- and three-point functions. This matching requires the evaluation of certain integrals over SUðNcÞ. In what follows, we calculate them symbolically,4 in the fundamental representation, and quote the main steps in the derivation. 1. Revising some relevant SUðNcÞintegrals According to [88], integrals over U∈SUðNcÞof the type ZdUUi1j1…UipjpU† k1l1 …U† knln;ðC1Þ 4A Python notebook using SymPy [87] where these calculations are carried out explicitly is publicly available at https:// github.com/avramescudana/sun_integrals. DANA AVRAMESCU et al. PHYS. REV. D 107, 114021 (2023) 114021-24