Progress on 3+1D Glasma simulations

We review our progress on 3+1D Glasma simulations to describe the earliest stages of heavy-ion collisions. In our simulations we include nuclei with finite longitudinal extent and describe the collision process as well as the evolution of the strongly interacting gluonic fields in the laboratory frame in 3+1 dimensions using the colored particle-in-cell method. This allows us to compute the 3+1 dimensional Glasma energy-momentum tensor, whose rapidity dependence can be compared to experimental pion multiplicity data from RHIC. An improved scheme cures the numerical Cherenkov instability and paves the way for simulations at higher energies used at LHC.


Introduction
QCD matter under extreme temperatures and densities in the form of the quark-gluon plasma is experimentally accessible in relativistic heavy-ion collisions. The possibility to conduct these experiments for a wide range of collision energies and baryon chemical potentials allows for the successive exploration of the QCD phase diagram. The highest collision energies have been achieved at LHC and RHIC, and lower collision energies are being explored in the Beam Energy Scan programs of RHIC [1] and upcoming programs at GSI FAIR [2] and JINR NICA [3]. The matter created in such collisions is initially very far from an ideal thermodynamic equilibrium. With the experimental progress also an improved theoretical understanding of the collision process from first principles is desirable.
The Color Glass Condensate (CGC) framework [4][5][6] provides such a theoretical basis for describing nuclear matter at ultrarelativistic energies. It is a classical effective field theory where hard partons within the nuclei act as sources for a e-mail: ipp@hep.itp.tuwien.ac.at b e-mail: dmueller@hep.itp.tuwien.ac.at (corresponding author) soft gluonic fields. In the simplest version, describing very large nuclei, the distribution of the color charges is given by the McLerran-Venugopalan (MV) model [7,8]. The preequilibrium stage that is created right after the collision is characterized by longitudinal color flux tubes and has been termed the Glasma [9]. In more sophisticated models like the IP-Glasma, the color charge distribution is based on fits to deep-inelastic-scattering data [10,11]. In combination with a subsequent hydrodynamical evolution, the IP-Glasma model is able to correctly reproduce many observables, including azimuthal anisotropies or event-by-event multiplicity distributions [12,13]. Nevertheless, the underlying Glasma evolution is commonly based on a boost-invariant formulation [14][15][16][17] where incoming nuclei are assumed to be Lorentzcontracted to infinitesimally thin discs. This assumption is justified for observables close to mid-rapidity at very high energies, but is a severe conceptual limitation when studying rapidity-dependent quantities or collisions at lower energies.
It is possible to break boost invariance by introducing fluctuations on top of boost invariant background fields [18][19][20]. Boost invariance is also broken by the JIMWLK evolution [21][22][23][24]. Apart from recent attempts [25,26], Glasma simulations using JIMWLK-based initial conditions still have to be performed in an effectively boost invariant manner [27]. The generated rapidity dependence of the Glasma can reproduce observables like gluon [27] and charged hadron multiplicities [26]. However, deviations from boost invariance may already arise at the classical level if one considers nuclei with finite extent in the beam direction [28]. A three-dimensional formulation has been introduced by using an extended source that is not quite aligned with the light cone [29][30][31], however in a way that violates the covariant conservation condition at order g 2 . An alternative approach to include such finite extent corrections has been developed for proton-nucleus collisions [32][33][34] which however is difficult to generalize to nucleusnucleus collisions.
In the following we review our progress towards an alternative approach of a 3+1D simulation for heavy-ion collisions including finite longitudinal extent of the nuclei [35][36][37][38][39]. The loss of boost invariance requires us to keep track of the hard color sources throughout the subsequent evolution after the collision. This is achieved using the colored particlein-cell method (CPIC), which has been originally developed to study aspects of the evolution of the quark-gluon plasma [40][41][42][43][44]. We apply this technique to study the collision process itself within the CGC framework. The simulation is performed in the laboratory frame and follows the nuclei throughout the collision process. Using this approach, we demonstrate that already a classical leading-order CGC simulation can give rise to a rapidity dependency consistent with experimental findings. Recent algorithmic developments will allow us to scrutinize our findings at even higher energies.

The 2+1D Glasma
From the viewpoint of the laboratory frame, a high energy nucleus is highly Lorentz-contracted along the beam axis and its dynamics are slowed down by time dilation. At collision energies available to RHIC and LHC, nuclei therefore appear to be almost infinitesimally thin, static discs moving at highly relativistic velocities. In CGC effective theory [4][5][6] the partons of high nuclei are split into hard and soft degrees of freedom. At leading order, hard partons are described in terms of classical color currents J μ and soft partons in terms of classical color fields A μ , whose dynamics are governed by the Yang-Mills (YM) field equations. For example, the classical color current associated with a nucleus moving along the negative x 3 = z axis (shown as nucleus "A" in Fig. 1) is given by where we have used light cone coordinates x ± = (x 0 ± x 3 )/ √ 2 and transverse coordinates x T = (x, y). The matrices t a are the traceless, hermitian generators of the color gauge group SU(N c ). The color charge density is given by ρ A (x + , x T ) and is assumed to be highly peaked around x + = 0 to account for the high Lorentz contraction of the nucleus. Additionally, ρ A (x + , x T ) is static in the sense that it does not explicitly depend on x − , which is a consequence of time dilation.
The color current J μ induces a classical color field A μ , which is to be determined from the non-linear YM equations.
Nucleus A Nucleus B Glasma z y x Fig. 1 Three-dimensional simulation of the collision of two nuclei. The distribution of the energy density between the two nuclei "A" and "B" right after the collision reveals the flux tube structure of the Glasma that develops between them. The simulation only covers a small part of the full collision in the transverse plane spanned by x and y. Adapted from [36] where D μ denotes the gauge covariant derivative. The non-Abelian field strength tensor is defined as where g is the YM coupling constant. If J μ is given by Eqs.
(2) can be solved analytically by employing covariant gauge ∂ μ A μ = 0. In this gauge choice, the field Eqs. simplify to where Δ T is the two-dimensional Laplace operator in the transverse plane. The other components of A μ A (x) vanish. This can be formally solved by inverting the Laplace operator: Infrared divergences are avoided by including a small regulator m, which dampens the long-range behavior of the color field. The length scale m −1 is usually identified with the confinement radius. It is also possible to regulate ultraviolet modes with a cutoff . The color field of the nucleus consists of purely transverse color-electric and -magnetic fields located near x + = 0. The color field of a high energy nucleus is therefore analogous to the Lorentz-boosted electromagnetic field of an ultrarelativistic charge. In CGC effective theory the color currents of nuclei are stochastic fields whose distribution is described by a probability functional W [ρ]. Expectation values of observables (expressed as functionals of the color fields O[A μ ]) are computed by averaging over all possible realizations of ρ: A simple and popular model for W [ρ] is the MV model [7,8], which assumes the functional to be Gaussian: where Z is a normalization constant and the function μ 2 (x + , x T ) describes the variance of the randomly distributed color charges. The description of nucleus "B" (see Fig. 1), which moves along the positive z axis is completely analogous.
A collision of two CGCs can be described by solving where the source is given by the combined color current of the colliding nuclei. We assume nucleus "A" to be centered around x + = 0 and nucleus "B" to be centered around x − = 0 such that the collision occurs at x + = x − = 0. A Minkowski diagram of the collision scenario is shown in Fig. 2. In the boost invariant limit, the color charge densities become effectively two-dimensional due to Lorentz contraction: This approximation is only correct in the limit of infinite collision energy. In this case, the space-time shown in Fig. 2 separates into four distinct regions. Invoking causality, the solution A μ (x) to Eq. (8) in regions I -III is given by the superposition of the single nucleus solutions The solution in the future light cone, which describes the Glasma, generally does not exist in closed form and has to be determined perturbatively or numerically. In the boost invariant limit however, the gauge field can be determined along the boundary of the light cone. Using proper time τ = √ 2x + x − and (space-time) rapidity η s = ln x − /x + /2 and employing temporal gauge A τ = 0 for τ ≥ 0, the color field at τ = 0 + is given by [14] A Here, the color fields where the lightlike Wilson lines are given by The initial conditions Eqs. (10) and (11) describe a highly anisotropic initial state consisting of purely longitudinal color-electric and -magnetic flux tubes [9,45]. Since these initial conditions do not depend on rapidity η s , the Glasma and any observables remain boost invariant for τ > 0. For charge densities which are not δ-shaped as in Eq. (9) and instead have finite longitudinal length along x ± , there are no rigorous derivations of generalized initial conditions. In order to allow for charge densities which explicitly break boost invariance, we have to move to a fully 3+1 dimensional description of the Glasma. In particular, it is necessary to solve the YM Eqs. (8) in a different way.

Simulating the Glasma in 3+1D
The motivation for relaxing Eq. (9) is to go beyond the approximation of infinite collisional energy and be able to describe observables such as the energy momentum tensor of the Glasma in a rapidity dependent manner. A simple generalization is given by where λ is a normalized function, which determines the longitudinal shape of the charge density. The width of λ should be directly related to the Lorentz-contracted diameter of the nucleus. It should be noted that Eq. (16) is a special case where the color structure does not depend on the longitudinal coordinate x ± and more general color charge densities A direct consequence of allowing for finite longitudinal extent is that the collision event is not just a single point in the Minkowski diagram (see Fig. 2), but an extended spacetime region in which the nuclei are able to interact. The nonperturbative nature of the color fields of the nuclei and the extended interaction time generally make analytical calculations in this space-time region intractable. Our approach to the 3+1D Glasma is therefore to simulate not just the evolution in the future light cone, i.e. region IV in Fig. 2, but the whole collision [35]. This setup is formulated in the laboratory frame using (t, z) coordinates and the time evolution is performed in the direction of t instead of proper time τ . The initial conditions of such a time evolution in t are specified in the following way: at some initial time t 0 < 0 sufficiently far away from the collision time t = 0, the color charge densities of the nuclei are non-overlapping in z. In LC gauge, the color field of each nucleus is given by where x = (x, y, z) is a three-dimensional coordinate vector. Equation (17) solves the classical YM Eqs. (2) with the current given by Eq. (16) [35]. Since the fields vanish exponentially fast between the two nuclei, the superposition of both color fields is a valid solution to the YM Eqs. (8) at time t 0 . Evolving this initial condition for the YM field to t > 0 yields a genuinely 3+1D description of the rapidity dependent Glasma. Equations (17) and (18) are compatible with the temporal gauge condition A 0 (t, x) = 0, ∀t ∈ R, which is a convenient gauge choice for an evolution along x 0 = t.
One of the main differences to the standard 2+1D Glasma is that in the laboratory frame the system not only consists of the color fields but also the color currents of the nuclei. In order to solve the YM Eqs. which include color currents, we make use of the CPIC method [40][41][42][43][44]. CPIC allows for consistent simulations of color charged point particles coupled to color fields on a lattice. For a comprehensive description of our numerical methods we refer to [35,39].
The numerical treatment of the fields follows the common real-time lattice gauge theory approach: by discretizing Minkowski space-time as a hypercubic lattice with spacings a i and time-step Δt, the YM Eqs. can be written in the standard leapfrog scheme where , x) is the chromo-electric field on the lattice and [X ] ah denotes the anti-hermitian, traceless part of a matrix X . The plaquette variables U i, j (t, x) are defined as The color currents of the nuclei j i (t, x), which enter on the right hand side of Eq. (19), require careful treatment. The main idea of using the CPIC method to describe collisions in the CGC framework is to replace the continuous color charge distributions ρ of the nuclei by a large number of auxiliary particles with time-dependent color charges Q k (t) such that the original color charge distribution is sufficiently well approximated on a lattice: where k is the particle index. Similar to the boost invariant case, we assume these auxiliary particles to be recoil-less. Thus, the trajectories of the particles x k (t) are fixed and not part of the dynamics of the system. The time-dependence of the color charges Q k (t) is determined from the discretized continuity Eq. Adapted from [38] which is the discrete analogue of gauge covariant continuity equation In our setup, the color charge Q k (t) of each particle is mapped to its nearest grid point on the spatial lattice in each time-step. Whenever the nearest grid point of a particular point charge changes from one lattice site x to a neighbouring lattice site y, parallel transport is applied to the color charge accordingly: where U x→y (t) is the appropriate gauge link connecting x and y. At the same time, the movement of the particle generates a color current j i (t, x) in accordance with Eq. (23). In CPIC, this treatment of particles is known as the nearestgrid-point scheme [40]. By evolving the color charges of the particles in this manner, the discretized field Eqs. (19) and (20) are solved consistently in the sense that the discrete Gauss law remains satisfied throughout the simulation and gauge covariance on the lattice is retained. In Eq. (26) the parallel transported electric field is given byẼ We find that numerically stable 3+1D Glasma simulations using the leapfrog scheme Eqs. (19) and (20) require high lattice resolution with particularly fine lattice spacing along the beam axis z. In practice, this can be computationally prohibitive. This problem is related to a subtle numerical instability inherent to the leapfrog scheme known as the numerical Cherenkov instability [46], which also affects traditional electromagnetic plasma simulations. This unphysical instability is caused by lattice artifacts which modify the propagation of wave modes on the lattice: high frequency modes propagate at a lower phase velocity compared to low frequency modes (i.e. numerical dispersion). In contrast, the auxiliary particles move at the speed of light by design. This situation is reminiscent of charged particles moving through a medium in which the in-medium speed of light is lower than the particle velocity, which leads to the well-known phenomenon of Cherenkov radiation. The particular lattice discretization used in Eqs. (19) and (20) leads to a similar, although purely numerical generation of Cherenkov radiation. Due to the numerical instability the nuclei are not able to propagate stably unless a very small lattice spacing is chosen along z, which reduces the mismatch in propagation velocities. Fortunately, these numerical problems can be solved by modifying the lattice discretization of the color fields. In [38] we show how an improved numerical scheme restores the correct propagation of wave modes along the beam axis such that artificial Cherenkov radiation is effectively avoided. This modification mainly amounts to replacing specific spa-tial plaquette terms (see Fig. 3a) with time-averaged generalizations of these terms. Figure 3b shows one example of such a time-averaged term. In comparison to the leapfrog scheme Eqs. (19) and (20), which is an explicit finite difference method, our numerical method is in the form of a system of semi-implicit equations, which is solved in an iterative manner. Our semi-implicit method therefore trades computational performance for numerical stability and accuracy.

Results
Here we review the most important results that have been obtained from simulations of collisions of nuclei with finite longitudinal extent using the standard leapfrog scheme [35][36][37]. In these simulations we use initial conditions of the form Eq. (16), where λ is a Gaussian of width L along z. We also assume μ 2 to be constant in the transverse plane. The longitudinal thickness L has been set to L = m N R/ √ s N N with collision energy √ s N N , nuclear radius R and nucleon mass m N ≈ 1GeV. The saturation momentum Q s grows with collision energy as Q 2 s ≈ √ s N N 0.25 GeV 2 [47][48][49]. The MV model parameter μ can be determined from 0.75 g 2 μ Q s with coupling g ≈ 2 [10]. The ultraviolet modes are regulated by = 10 GeV, and we vary the infrared regulator m in the range from 0.2 to 0.8 GeV to check for its dependency. The simulations presented use the gauge group SU(2) instead of SU(3), which should give qualitatively comparable results [50]. The simulations have been performed on a lattice with 2048 × 192 2 cells with finer resolution along the longitudinal direction. The simulation box corresponds to a volume of (6 fm) 3 .
The main observables can be obtained from the components of the energy-momentum tensor T μν which are extracted from the gluon fields of the simulation. We average over 15 collision events to obtain the expectation value T μν . Due to the symmetries of the MV model, certain components vanish, and the remaining contributions of the averaged energy-momentum tensor are given by where ε is the energy density in the laboratory frame, p L and p T are the longitudinal and transverse pressure components and S L is the longitudinal component of the Poynting vector. The energy density ε as given in the laboratory frame is depicted in Fig. 1. By diagonalizing the energy-momentum tensor, we obtain the local rest frame energy density ε loc , which can be expressed in terms of proper time τ = √ t 2 − z 2 and space-time rapidity η s = ln [(t − z) / (t + z)] /2. The black solid lines show the rapidity profiles as extracted from of our simulations, which can be fitted to a Gaussian shape (dashed continuing lines). The profiles have been extracted at τ 0 = 1 fm/c where the Glasma turns into the QGP. Already at times τ 0 0.3 fm/c, the system enters a free-streaming evolution with longitudinal velocity v z ≈ z/t and the shape of the profile does not change anymore [37]. The width of the profiles depends on the energy √ s N N and becomes flatter with increasing energy as expected from the recovery of boost invariance at higher energies [36]. We also find a strong dependency on the infrared regulator m, where higher values of m make the rapidity profiles narrower. While this strong dependence on the infrared regulator seems unexpected, it may indicate that the screening length λ D ∝ 1/m plays an important role in generating a deviation from boost invariance. It is not only the longitudinal thickness L, but the dimensionless ratio L/λ D , which seems to determine the shape of the profiles. Interestingly, the simulation results agree well with measured rapidity profiles of pion multiplicities at RHIC [51] which are indicated by the blue data points in Fig. 4. At these energies, the experimental results also agree with the Gaussian rapidity profile as obtained by the hydrodynamic Landau model [52] where the width is given by σ Landau = √ ln γ with the Lorentz gamma factor γ . For other particle species, the experimental momentum rapidity distribution of particle multiplicities deviates from the Landau picture but can still be fitted to a Gaussian distribution with slightly larger width [53]. However, it should be noted that especially at higher energies the Landau model predicts rapidity profiles which are too narrow as measured by the ALICE collaboration [53]. Whether our simulations of the 3+1D Glasma can describe these wider profiles at LHC energies will be the topic of future work.
It is important to note that Fig. 4 shows the rapidity profiles of two different quantities: the local energy density ε loc (τ, η s ) as a function of space-time rapidity η s and the distribution of particle multiplicities d N/dy as a function of momentum rapidity y. This comparison is justified because our simulations show [36,37] that the 3+1D Glasma settles into a state of free-streaming flow v z ≈ z/t and vanishing longitudinal pressure p L ≈ 0, similar to the 2+1D Glasma. Freestreaming flow implies that -ignoring a subsequent hydrodynamical phase -we can identify momentum rapidity y with space-time rapidity η s , similar to the case of the Bjorken model [54]. It also implies that T μν becomes diagonal in the (τ, η s ) frame. Due to vanishing longitudinal pressure, the energy-momentum tensor in the rest frame is anisotropic and reads where we used 2 p T = ε loc due to conformal symmetry. In contrast to the Bjorken model, there are a priori no restrictions on the rapidity dependence of the energy density ε loc (τ, η s ). This can be checked using energy-momentum conservation which simply reduces to Equation (31) leads to ε loc ∝ 1/τ , but does not impose any additional constraints on the rapidity dependence of the energy density. On the other hand, the Bjorken model requires that T μν is isotropic in the rest frame, which combined with Eq. (30), then leads to boost invariance. As a first approximation, it is therefore reasonable to compare the space-time rapidity profile of ε loc to the momentum rapidity profile of charged particles as shown in Fig. 4. For a more quantitative comparison, our results should be used as input for a subsequent hydrodynamic simulation, which may slightly increase the width of the profiles [27].
To better understand how rapidity dependence of observables develops in our simulations, we can look at the Space-time distribution of the normalized transverse pressure p T (t, z) / p T (0, 0) [36]. The parameters are the same as in Fig. 4 with m = 0.2 GeV. The transverse pressure corresponds to longitudinal chromo-magnetic and -electric fields and thus to the the longitudinal component of the energy density ε L (t, z) . Contrary to the boost-invariant case, this quantity falls off steeply along the boundary of the light cone transverse pressure p T (z, t) in Fig. 5. From the energymomentum tensor, one finds that the transverse pressure is linked to the longitudinal fields of the Glasma, whereas longitudinal pressure involves both, transverse and longitudinal field components [35] where the square implies a summation over color indices. In the boost invariant case, the initial conditions of the Glasma are specified at τ = 0. This corresponds to the boundary of the forward light cone, where the longitudinal fields would be constant. In contrast, in our simulation the longitudinal fields are peaked around the collision region at t ∼ z ∼ 0 and decrease rather quickly along the light cone boundaries. Accordingly, there is less Glasma being produced at larger values of rapidity η s which produces the Gaussian profiles. It is interesting to see that our weak coupling results are in qualitative agreement with strong coupling results from holographic models of heavy-ion collisions that also exhibit similar transverse pressure distributions [55,56].

Conclusions and outlook
In this paper we reviewed our progress on 3+1D Glasma simulations. Our simulations allow to explore the creation of the Glasma in heavy-ion collisions beyond the commonly assumed boost-invariant case. We do this by introducing a finite longitudinal extent for the incoming nuclei corresponding to realistic Lorentz contractions as found for example at RHIC. Without the usual simplifications of boost-invariance, we have to keep the color currents of the hard partons in the Glasma simulation. This is achieved using CPIC in the lab-oratory frame. Using the MV model, we demonstrated that our approach can give rise to Gaussian rapidity profiles in the energy density. These profiles depend on the energy of the incoming nuclei, but also on an infrared regulator. For energies used at RHIC we obtain qualitative agreement of these profiles [37]. This is remarkable as it shows that boost invariance can be broken already at the classical level if the longitudinal structure is properly taken into account. This nicely complements findings from holographic models where similar profiles can be found [55,56]. Algorithmic improvements in the form of a new semiimplicit solver [38] will allow for further explorations of our boost-invariance breaking simulations. By modifying the standard Wilson gauge action we achieve a dispersion-free propagation along the longitudinal direction which cures the numerical Cherenkov instability which has plagued previous simulations. This sets the basis for more accurate and larger simulations valid for larger ranges of rapidity which are necessary for a comparison of rapidity profiles at LHC collision energies. One crucial aspect that shall be studied in this context is the role of longitudinal color fluctuations. These are usually approximated as an infinitely thin stack of uncorrelated sheets of color charge [57]. Dispersion-free propagation will allow for the fine-grained simulation of collisions and the study of the effect of internal longitudinal color structures on the creation of the Glasma. In principle, it should be straightforward to also include more realistic sub-nucleonic color structure in the transverse direction as is the case in the IP-Glasma model [10,11]. Here, one is essentially limited by the large computational requirements of such three-dimensional simulations.
On a more conceptual level, it would be interesting to better understand the relation between the boost-invariance breaking that we find at the leading classical order and a similar breaking that can be found at next-to-leading order from the JIMWLK evolution [27]. It would be a highly desirable but presumably very non-trivial task to generalize the JIMWLK evolution Eqs. to be applicable to threedimensional color distributions. Another extension to our work would be the deviation from the eikonal approximation and the inclusion of dynamical colored particles. This could also be a way to accommodate three-dimensional extensions of the calculation of quantities like energy loss or momentum broadening from the Glasma [58].