Progress on 3+1D Glasma simulations

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 Tlusty:2018rif and upcoming programs at GSI FAIR Senger:2020pzs and JINR NICA Geraksiev:2019fon . 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 Gelis:2010nm ; Gelis:2012ri ; Gelis:2015gza 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 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 McLerran:1993ka ; McLerran:1993ni . The pre-equilibrium stage that is created right after the collision is characterized by longitudinal color flux tubes and has been termed the Glasma Lappi:2006fp . In more sophisticated models like the IP-Glasma, the color charge distribution is based on fits to deep-inelastic-scattering data Schenke:2012hg ; Schenke:2012wb . 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 Gale:2012rq ; Snellings:2011sz . Nevertheless, the underlying Glasma evolution is commonly based on a boost-invariant formulation Kovner:1995ja ; Kovner:1995ts ; Krasnitz:1998ns ; Lappi:2011ju where incoming nuclei are assumed to be Lorentz-contracted 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 Gelis:2013rba ; Fukushima:2011nq ; Berges:2012cj . Boost invariance is also broken by the JIMWLK evolution JalilianMarian:1996xn ; Iancu:2000hn ; Mueller:2001uk ; Ferreiro:2001qy . Apart from recent attempts McDonald:2018wql ; McDonald:2020oyf , Glasma simulations using JIMWLK-based initial conditions still have to be performed in an effectively boost invariant manner Schenke:2016ksl . The generated rapidity dependence of the Glasma can reproduce observables like gluon Schenke:2016ksl and charged hadron multiplicities McDonald:2020oyf . However, deviations from boost invariance may already arise at the classical level if one considers nuclei with finite extent in the beam direction Ozonder:2013moa . A three-dimensional formulation has been introduced by using an extended source that is not quite aligned with the light cone Lam:1999wu ; Lam:2000nz ; Ozonder:2012vw , 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 Altinoluk:2014oxa ; Altinoluk:2015gia ; Altinoluk:2015xuy which however is difficult to generalize to nucleus-nucleus 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 Gelfand:2016yho ; Ipp:2017lho ; Ipp:2017uxo ; Ipp:2018hai ; Muller:2019bwd . 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 particle-in-cell method (CPIC), which has been originally developed to study aspects of the evolution of the quark-gluon plasma Hu:1996sf ; Moore:1997sn ; Dumitru:2005hj ; Dumitru:2006pz ; Schenke:2008gg . 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.

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 Gelis:2010nm ; Gelis:2012ri ; Gelis:2015gza 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^{\mu}$ and soft partons in terms of classical color fields $A^{\mu}$, whose dynamics are governed by the Yang-Mills (YM) field eqs. 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

$$J_{\tiny\mathrm{A}}^{\mu}(x)=\delta^{\mu}_{-}\rho_{\tiny\mathrm{A}}^{a}(x^{+},\mathbf{x}_{T})\mathbf{t}^{a},$$

(1)

where we have used light cone coordinates $x^{\pm}=(x^{0}\pm x^{3})/\sqrt{2}$ and transverse coordinates $\mathbf{x}_{T}=(x,y)$. The matrices $\mathbf{t}^{a}$ are the traceless, hermitian generators of the color gauge group SU($N_{c}$). The color charge density is given by $\rho_{\tiny\mathrm{A}}(x^{+},\mathbf{x}_{T})$ and is assumed to be highly peaked around $x^{+}=0$ to account for the high Lorentz contraction of the nucleus. Additionally, $\rho_{\tiny\mathrm{A}}(x^{+},\mathbf{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^{\mu}$ induces a classical color field $A^{\mu}$, which is to be determined from the non-linear YM eqs.

$$D_{\mu}F^{\mu\nu}=\partial_{\mu}F^{\mu\nu}+ig\left[A_{\mu},F^{\mu\nu}\right]=J^{\nu},$$

(2)

where $D_{\mu}$ denotes the gauge covariant derivative. The non-Abelian field strength tensor is defined as

$$F^{\mu\nu}=\partial^{\mu}A^{\nu}-\partial^{\nu}A^{\mu}+ig\left[A^{\mu},A^{\nu}\right],$$

(3)

where $g$ is the YM coupling constant. If $J^{\mu}$ is given by eq. (1), the YM eqs. (2) can be solved analytically by employing covariant gauge ${\partial_{\mu}A^{\mu}=0}$. In this gauge choice, the field eqs. simplify to

$$-\Delta_{T}A_{\tiny\mathrm{A}}^{-}(x)=\rho_{\tiny\mathrm{A}}^{a}(x^{+},\mathbf{x}_{T})\mathbf{t}^{a},$$

(4)

where $\Delta_{T}$ is the two-dimensional Laplace operator in the transverse plane. The other components of $A^{\mu}_{\tiny\mathrm{A}}(x)$ vanish. This can be formally solved by inverting the Laplace operator:

$$A^{\mu}_{\tiny\mathrm{A}}(x^{+},\mathbf{x}_{T})=\delta^{\mu}_{-}\intop^{\Lambda}\frac{d^{2}\mathbf{k}_{T}}{(2\pi)^{2}}\frac{\rho_{\tiny\mathrm{A}}^{a}(x^{+},\mathbf{k}_{T})}{\mathbf{k}^{2}_{T}+m^{2}}e^{i\mathbf{k_{T}}\cdot\mathbf{x_{T}}}\mathbf{t}^{a}.$$

(5)

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 $\Lambda$. 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[\rho]$. Expectation values of observables (expressed as functionals of the color fields $\mathcal{O}[A_{\mu}]$) are computed by averaging over all possible realizations of $\rho$:

$$\langle\mathcal{O}\rangle=\int\mathcal{D}\rho\,\mathcal{O}[A_{\mu}[\rho]]\,W[\rho].$$

(6)

A simple and popular model for $W[\rho]$ is the MV model McLerran:1993ka ; McLerran:1993ni , which assumes the functional to be Gaussian:

$$W[\rho]=\frac{1}{Z}\exp{\left(\!-\!\int\!d^{2}\mathbf{x}_{T}dx^{+}\frac{\rho^{a}(x^{+},\mathbf{x}_{T})\rho^{a}(x^{+},\mathbf{x}_{T})}{2g^{2}\mu^{2}(x^{+},\mathbf{x}_{T})}\right)},$$

(7)

where $Z$ is a normalization constant and the function $\mu^{2}(x^{+},\mathbf{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

$$D_{\mu}F^{\mu\nu}(x)=J^{\nu}_{\tiny\mathrm{A}}(x)+J^{\nu}_{\tiny\mathrm{B}}(x),$$

(8)

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:

$$\rho_{({\tiny\mathrm{A}},{\tiny\mathrm{B}})}(x)\approx\delta(x^{\pm})\rho_{({\tiny\mathrm{A}},{\tiny\mathrm{B}})}(\mathbf{x}_{T}).$$

(9)

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^{\mu}(x)$ to eq. (8) in regions I - III is given by the superposition of the single nucleus solutions $A^{\mu}_{\tiny\mathrm{A}}(x)$ and $A^{\mu}_{\tiny\mathrm{B}}(x)$. 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 $\tau=\sqrt{2x^{+}x^{-}}$ and (space-time) rapidity $\eta_{s}=\ln\left(x^{-}/x^{+}\right)/2$ and employing temporal gauge $A^{\tau}=0$ for $\tau\geq 0$, the color field at $\tau=0^{+}$ is given by Kovner:1995ja

	$\displaystyle A^{i}(\tau=0^{+},\mathbf{x}_{T})$	$\displaystyle=\alpha^{i}_{\tiny\mathrm{A}}(\mathbf{x})+\alpha^{i}_{\tiny\mathrm{B}}(\mathbf{x}),$		(10)
	$\displaystyle A^{\eta}(\tau=0^{+},\mathbf{x}_{T})$	$\displaystyle=\frac{ig}{2}\left[\alpha^{i}_{\tiny\mathrm{A}}(\mathbf{x}),\alpha^{i}_{\tiny\mathrm{B}}(\mathbf{x})\right].$		(11)

Here, the color fields $\alpha^{i}_{({\tiny\mathrm{A}},{\tiny\mathrm{B}})}$ are the light cone (LC) gauge $A^{\mp}=0$ solutions

	$\displaystyle A^{i}_{({\tiny\mathrm{A}},{\tiny\mathrm{B}})}(x^{\pm},\mathbf{x}_{T})$	$\displaystyle=\theta(x^{\pm})\alpha^{i}_{({\tiny\mathrm{A}},{\tiny\mathrm{B}})}(\mathbf{x}_{T}),$		(12)
	$\displaystyle\alpha^{i}_{({\tiny\mathrm{A}},{\tiny\mathrm{B}})}(\mathbf{x}_{T})$	$\displaystyle=\frac{1}{ig}V_{({\tiny\mathrm{A}},{\tiny\mathrm{B}})}(\mathbf{x})\partial^{i}V^{\dagger}_{({\tiny\mathrm{A}},{\tiny\mathrm{B}})}(\mathbf{x}),$		(13)

where the lightlike Wilson lines are given by

	$\displaystyle V_{({\tiny\mathrm{A}},{\tiny\mathrm{B}})}^{\dagger}(\mathbf{x}_{T})$	$\displaystyle=\lim_{x^{\pm}\rightarrow\infty}V_{({\tiny\mathrm{A}},{\tiny\mathrm{B}})}^{\dagger}(x^{\pm},\mathbf{x}_{T}),$		(14)
	$\displaystyle V_{({\tiny\mathrm{A}},{\tiny\mathrm{B}})}^{\dagger}(x^{\pm},\mathbf{x}_{T})$	$\displaystyle=\mathcal{P}\exp\bigg{(}ig\!\!\intop_{-\infty}^{x^{\pm}}\!\!dx^{\prime\pm}A_{({\tiny\mathrm{A}},{\tiny\mathrm{B}})}^{\mp}(x^{\prime\pm},\mathbf{x}_{T})\bigg{)}.$		(15)

The initial conditions eqs. (10) and (11) describe a highly anisotropic initial state consisting of purely longitudinal color-electric and -magnetic flux tubes Fries:2006pv ; Lappi:2006fp . Since these initial conditions do not depend on rapidity $\eta_{s}$, the Glasma and any observables remain boost invariant for $\tau>0$. For charge densities which are not $\delta$-shaped as in eq. (9) and instead have finite longitudinal length along $x^{\pm}$, 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.

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

$$\rho_{({\tiny\mathrm{A}},{\tiny\mathrm{B}})}(x)\approx\lambda(x^{\pm})\rho_{({\tiny\mathrm{A}},{\tiny\mathrm{B}})}(\mathbf{x}_{T}),$$

(16)

where $\lambda$ is a normalized function, which determines the longitudinal shape of the charge density. The width of $\lambda$ 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^{\pm}$ and more general color charge densities $\rho_{({\tiny\mathrm{A}},{\tiny\mathrm{B}})}(x^{\pm},\mathbf{x}_{T})$ are possible as well.

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 space-time region in which the nuclei are able to interact. The non-perturbative 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 Gelfand:2016yho . 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 $\tau$. 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

$$A^{i=x,y}_{({\tiny\mathrm{A}},{\tiny\mathrm{B}})}(t_{0},\mathbf{x})=\frac{1}{ig}V_{({\tiny\mathrm{A}}{\tiny\mathrm{B}})}(t_{0},\mathbf{x})\partial^{i}V^{\dagger}_{({\tiny\mathrm{A}}{\tiny\mathrm{B}})}(t_{0},\mathbf{x}),$$

(17)

where $\mathbf{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) Gelfand:2016yho . Since the fields vanish exponentially fast between the two nuclei, the superposition of both color fields

$$A^{\mu}(t_{0},\mathbf{x})=\delta^{\mu}_{i=x,y}\left(A^{i}_{\tiny\mathrm{A}}(t_{0},\mathbf{x})+A^{i}_{\tiny\mathrm{B}}(t_{0},\mathbf{x})\right),$$

(18)

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,\mathbf{x})=0$, $\forall t\in\mathbb{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 Hu:1996sf ; Moore:1997sn ; Dumitru:2005hj ; Dumitru:2006pz ; Schenke:2008gg . 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 Gelfand:2016yho and Muller:2019bwd .

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 $\Delta t$, the YM eqs. can be written in the standard leapfrog scheme

	$\displaystyle E_{i}(t+\frac{\Delta t}{2},\mathbf{x})=$	$\displaystyle-\sum_{j}\frac{\Delta t}{(a^{j})^{2}}\left[U_{i,j}(t,\mathbf{x})+U_{i,-j}(t,\mathbf{x})\right]_{\mathrm{ah}}$
		$\displaystyle+\Delta t\,j_{i}(t,\mathbf{x})+E_{i}(t-\frac{\Delta t}{2},\mathbf{x}),$		(19)
	$\displaystyle U_{i}(t+\Delta t,\mathbf{x})=$	$\displaystyle\exp\bigg{(}i\Delta t\,E_{i}(t+\frac{\Delta t}{2},\mathbf{x})\bigg{)}U_{i}(t,\mathbf{x}),$		(20)

where $U_{i}(t,\mathbf{x})\simeq\exp\big{(}iga^{i}A_{i}(t,\mathbf{x})\big{)}$ are the gauge links, $E_{i}(t,\mathbf{x})$ is the chromo-electric field on the lattice and $\left[X\right]_{\mathrm{ah}}$ denotes the anti-hermitian, traceless part of a matrix $X$. The plaquette variables $U_{i,j}(t,\mathbf{x})$ are defined as

$$U_{i,j}(t,\mathbf{x})=U_{i}(t,\mathbf{x})U_{j}(t,\mathbf{x}+\hat{a}_{i})U^{\dagger}_{i}(t,\mathbf{x}+\hat{a}_{i})U^{\dagger}_{j}(t,\mathbf{x}).$$

(21)

The color currents of the nuclei $j_{i}(t,\mathbf{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 $\rho$ 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:

$$\rho(t,\mathbf{x})\approx\sum_{k}Q_{k}(t)\delta^{(3)}(\mathbf{x}-\mathbf{x}_{k}(t)),$$

(22)

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 $\mathbf{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.

		$\displaystyle\frac{\rho(t+\frac{\Delta t}{2},\mathbf{x})-\rho(t-\frac{\Delta t}{2},\mathbf{x})}{\Delta t}=$
		$\displaystyle\quad\sum_{i}\frac{j_{i}(t,\mathbf{x})-U^{\dagger}_{i}(t,\mathbf{x}-\hat{a}_{i})j_{i}(t,\mathbf{x}-\hat{a}_{i})U_{i}(t,\mathbf{x}-\hat{a}_{i})}{a^{i}},$		(23)

which is the discrete analogue of gauge covariant continuity eq.

$$D_{\mu}J^{\mu}(t,\mathbf{x})=0.$$

(24)

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 $\mathbf{x}$ to a neighbouring lattice site $\mathbf{y}$, parallel transport is applied to the color charge accordingly:

$$Q_{k}(t+\frac{\Delta t}{2})=U^{\dagger}_{\mathbf{x}\rightarrow\mathbf{y}}(t)Q_{k}(t-\frac{\Delta t}{2})U_{\mathbf{x}\rightarrow\mathbf{y}}(t),$$

(25)

where $U_{\mathbf{x}\rightarrow\mathbf{y}}(t)$ is the appropriate gauge link connecting $\mathbf{x}$ and $\mathbf{y}$. At the same time, the movement of the particle generates a color current $j_{i}(t,\mathbf{x})$ in accordance with eq. (3). In CPIC, this treatment of particles is known as the nearest-grid-point scheme Hu:1996sf . 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

		$\displaystyle\sum_{i}\bigg{(}E_{i}(t+\frac{\Delta t}{2},\mathbf{x})-\tilde{E}_{i}(t+\frac{\Delta t}{2},\mathbf{x}-\hat{a}_{i})\bigg{)}=$
		$\displaystyle\qquad\rho(t+\frac{\Delta t}{2},\mathbf{x}),$		(26)

remains satisfied throughout the simulation and gauge covariance on the lattice is retained. In eq. (3) the parallel transported electric field is given by

	$\displaystyle\tilde{E}_{i}(t+\frac{\Delta t}{2},\mathbf{x}-\hat{a}_{i})=$
	$\displaystyle\qquad U^{\dagger}_{i}(t,\mathbf{x}-\hat{a}_{i})E_{i}(t+\frac{\Delta t}{2},\mathbf{x}-\hat{a}_{i})U_{i}(t,\mathbf{x}-\hat{a}_{i}).$		(27)

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 GODFREY1974504 , 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 Ipp:2018hai 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 spatial plaquette terms (see fig. 3 (a)) with time-averaged generalizations of these terms. Figure 3 (b) 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 eqs., which is solved in an iterative manner. Our semi-implicit method therefore trades computational performance for numerical stability and accuracy.

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 Gelfand:2016yho ; Ipp:2017lho ; Ipp:2017uxo . In these simulations we use initial conditions of the form eq. (16), where $\lambda$ is a Gaussian of width $L$ along $z$. We also assume $\mu^{2}$ to be constant in the transverse plane. The longitudinal thickness $L$ has been set to $L=m_{N}R/\sqrt{s_{NN}}$ with collision energy $\sqrt{s_{NN}}$, nuclear radius $R$ and nucleon mass $m_{N}\approx 1\>\mbox{GeV}$. The saturation momentum $Q_{s}$ grows with collision energy as $Q_{s}^{2}\approx\left(\sqrt{s_{NN}}\right)^{0.25}\,\mbox{GeV}^{2}$ Lappi:2006hq ; Lappi:2007ku ; Kharzeev:2001gp . The MV model parameter $\mu$ can be determined from $0.75\,g^{2}\mu\simeq Q_{s}$ with coupling $g\approx 2$ Schenke:2012hg . The ultraviolet modes are regulated by $\Lambda=10\>\mbox{GeV}$, and we vary the infrared regulator $m$ in the range from 0.2 to $0.8\,\mbox{GeV}$ to check for its dependency. The simulations presented use the gauge group $\mbox{SU}(2)$ instead of SU$(3)$, which should give qualitatively comparable results Ipp:2010uy . The simulations have been performed on a lattice with $2048\times 192^{2}$ cells with finer resolution along the longitudinal direction. The simulation box corresponds to a volume of $\left(6\,\mbox{fm}\right)^{3}$.

The main observables can be obtained from the components of the energy-momentum tensor $T^{\mu\nu}$ which are extracted from the gluon fields of the simulation. We average over $15$ collision events to obtain the expectation value $\langle T^{\mu\nu}\rangle$. Due to the symmetries of the MV model, certain components vanish, and the remaining contributions of the averaged energy-momentum tensor are given by

$$\langle T^{\mu\nu}\rangle=\left(\begin{array}[]{cccc}\langle\varepsilon\rangle&0&0&\langle S_{L}\rangle\\ 0&\langle p_{T}\rangle&0&0\\ 0&0&\langle p_{T}\rangle&0\\ \langle S_{L}\rangle&0&0&\langle p_{L}\rangle\end{array}\right),$$

(28)

where $\langle\varepsilon\rangle$ is the energy density in the laboratory frame, $\langle p_{L}\rangle$ and $\langle p_{T}\rangle$ are the longitudinal and transverse pressure components and $\langle S_{L}\rangle$ is the longitudinal component of the Poynting vector. The energy density $\langle\varepsilon\rangle$ 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 $\langle\varepsilon_{\text{loc}}\rangle$, which can be expressed in terms of proper time $\tau=\sqrt{t^{2}-z^{2}}$ and space-time rapidity $\eta_{s}=\ln\left[\left(t-z\right)/\left(t+z\right)\right]/2$.

Figure 4 shows the local energy density $\varepsilon_{loc}(\tau_{0},\eta_{s})$ as a function of space-time rapidity $\eta_{s}$ for $\sqrt{s_{NN}}=200\,\mbox{GeV}$. 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 $\tau_{0}=1\,\mbox{fm}/c$ where the Glasma turns into the QGP. Already at times $\tau_{0}\gtrsim 0.3\,\mbox{fm}/c$, the system enters a free-streaming evolution with longitudinal velocity $v_{z}\approx z/t$ and the shape of the profile does not change anymore Ipp:2017uxo . The width of the profiles depends on the energy $\sqrt{s_{NN}}$ and becomes flatter with increasing energy as expected from the recovery of boost invariance at higher energies Ipp:2017lho . 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 $\lambda_{D}\propto 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/\lambda_{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 Bearden:2004yx which are indicated by the blue data points in figure 4. At these energies, the experimental results also agree with the Gaussian rapidity profile as obtained by the hydrodynamic Landau model Landau:1953gs where the width is given by $\sigma_{\text{Landau}}=\sqrt{\ln\gamma}$ with the Lorentz gamma factor $\gamma$. 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 Abbas:2013bpa . 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 Abbas:2013bpa . 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 figure 4 shows the rapidity profiles of two different quantities: the local energy density $\varepsilon_{\mathrm{loc}}(\tau,\eta_{s})$ as a function of space-time rapidity $\eta_{s}$ and the distribution of particle multiplicities $dN/dy$ as a function of momentum rapidity $y$. This comparison is justified because our simulations show Ipp:2017lho ; Ipp:2017uxo that the 3+1D Glasma settles into a state of free-streaming flow $v_{z}\approx z/t$ and vanishing longitudinal pressure $p_{L}\approx 0,$ similar to the 2+1D Glasma. Free-streaming flow implies that – ignoring a subsequent hydrodynamical phase – we can identify momentum rapidity $y$ with space-time rapidity $\eta_{s}$, similar to the case of the Bjorken model Bjorken:1982qr . It also implies that $T^{\mu\nu}$ becomes diagonal in the $(\tau,\eta_{s})$ frame. Due to vanishing longitudinal pressure, the energy-momentum tensor in the rest frame is anisotropic and reads

$$\langle T^{\mu\nu}\rangle=\mathrm{diag}(\langle\varepsilon_{\mathrm{loc}}\rangle,\,\langle\varepsilon_{\mathrm{loc}}\rangle/2,\,\langle\varepsilon_{\mathrm{loc}}\rangle/2,\,0)^{\mu\nu},$$

(29)

where we used $2\,p_{T}=\varepsilon_{\mathrm{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 $\varepsilon_{\mathrm{loc}}(\tau,\eta_{s})$. This can be checked using energy-momentum conservation

$$\nabla_{\mu}T^{\mu\nu}=0,$$

(30)

which simply reduces to

$$\partial_{\tau}\varepsilon_{\mathrm{loc}}(\tau,\eta_{s})=-\varepsilon_{\mathrm{loc}}(\tau,\eta_{s})/\tau.$$

(31)

Equation (31) leads to $\varepsilon_{\mathrm{loc}}\propto 1/\tau$, 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^{\mu\nu}$ 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 $\varepsilon_{\mathrm{loc}}$ to the momentum rapidity profile of charged particles as shown in figure 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 Schenke:2016ksl .

To better understand how rapidity dependence of observables develops in our simulations, we can look at the transverse pressure $\langle p_{T}(z,t)\rangle$ in figure 5. From the energy-momentum 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 Gelfand:2016yho

	$\displaystyle\langle p_{T}\rangle$	$\displaystyle=$	$\displaystyle\frac{1}{2}\langle E_{L}^{2}+B_{L}^{2}\rangle,$		(32)
	$\displaystyle\langle p_{L}\rangle$	$\displaystyle=$	$\displaystyle\frac{1}{2}\langle E_{T}^{2}+B_{T}^{2}-E_{L}^{2}-B_{L}^{2}\rangle,$		(33)

where the square implies a summation over color indices. In the boost invariant case, the initial conditions of the Glasma are specified at $\tau=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\sim z\sim 0$ and decrease rather quickly along the light cone boundaries. Accordingly, there is less Glasma being produced at larger values of rapidity $\eta_{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 Casalderrey-Solana:2013aba ; vanderSchee:2015rta .

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 laboratory 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 Ipp:2017uxo . 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 Casalderrey-Solana:2013aba ; vanderSchee:2015rta .

Algorithmic improvements in the form of a new semi-implicit solver Ipp:2018hai 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 Fukushima:2007ki . 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 Schenke:2012wb ; Schenke:2012hg . 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 Schenke:2016ksl . It would be a highly desirable but presumably very non-trivial task to generalize the JIMWLK evolution eqs. to be applicable to three-dimensional 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 Ipp:2020mjc .

The authors thank for a Short Term Scientific Mission (STSM) to CEA Saclay for scientific discussions with Jean-Paul Blaizot, François Gelis and Edmond Iancu within the framework of the COST action CA15213 “Theory of hot matter and relativistic heavy-ion collisions (THOR)”. This work has been supported by the Austrian Science Fund FWF No. P28352.