Broken boost invariance in the Glasma via finite nuclei thickness

We simulate the creation and evolution of non-boost-invariant Glasma in the early stages of heavy ion collisions within the color glass condensate framework. This is accomplished by extending the McLerran-Venugopalan model to include a parameter for the Lorentz-contracted but finite width of the nucleus in the beam direction. We determine the rapidity profile of the Glasma energy density, which shows deviations from the boost-invariant result. Varying the parameters both broad and narrow profiles can be produced. We compare our results to experimental data from RHIC and find surprising agreement.

Heavy ion collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) provide insight into the properties of nuclear matter under extreme conditions. The evolution of the Quark-Gluon Plasma (QGP) that is created in such collisions is well described by relativistic viscous hydrodynamics [1,2]. A first principles description of the initial state of heavyion collisions is provided by the Color Glass Condensate (CGC) framework [3][4][5]. The CGC is a classical effective field theory for nuclear matter at ultrarelativistic energies. Models such as the IP-Glasma [6,7] in combination with hydrodynamics are able to correctly reproduce azimuthal anisotropies and event-by-event multiplicity distributions [8,9]. Furthermore, the CGC can explain long-range rapidity correlations like the ridge [10,11].
A Gaussian shaped rapidity profile of particle multiplicity can be found in experiments covering various energy ranges, from LHC [12] to RHIC Beam Energy Scan [13,14]. This shape is well explained by the Landau model [15] up to RHIC energies [13], which assumes full stopping of the colliding nuclei. The Landau model is in contrast to the Bjorken model [16] which relies on approximate boost invariance. A Gaussian profile has also been found in holographic calculations of colliding shock waves [17][18][19].
In its original formulation collisions in the CGC picture are assumed to be boost-invariant [20][21][22][23] and were thus only understood as an approximation valid close to midrapidity. This approach implicitly assumes infinitely thin Lorentz-contracted nuclei and entails a classical, boost-invariant evolution of the Glasma at leading order.
Only at the next-to-leading order the boost invariance is broken by a change of the initial conditions through JIMWLK evolution [24][25][26][27], and non-boost-invariant rapidity profiles can be obtained. Recently such rapidity dependencies have been found to agree reasonably well with experimental data where observables like charged particle multiplicities show a Gaussian rapidity profile [28,29]. On the other hand it has been suggested that if one considers nuclei with finite extent in the beam direction, deviations from boost invariance may arise already Nucleus B Nucleus A Glasma x y z Figure 1. A 3+1 dimensional colored particle-in-cell simulation of the collision of two thick sheets of relativistic nuclear matter. The simulation box covers a small part of the full transverse extent of the nuclei in the x, y plane. This figure shows a density plot of the energy density of both nuclei A and B, and the three-dimensional Glasma that is created in the collision.
at the classical level [30]. In the case of proton-nucleus collisions methods have been developed to systematically include finite width corrections of the nucleus [31][32][33]. However, so far there has been no consistent simulation of the subsequent three-dimensional evolution for heavy-ion collisions even at the classical level.
In this letter, we show that Gaussian rapidity profiles of energy density can arise already from 3+1 dimensional purely classical CGC simulations, if incoming nuclei have a finite extent in the beam direction. As one would expect, the Gaussian profiles become broader at higher collision energy. In principle, we can cover the wide range from very thin nuclei with almost boostinvariant behavior to thick nuclei at low collision energies and narrow Gaussian profiles. We simulate the collision in the laboratory frame, see Fig. 1, which makes it necessary to include the propagating nuclei already before and during the collision.
The nuclei in the CGC picture consist of hard partons that are surrounded by soft gluons. The hard partons can be described as classical color charges moving at the speed of light, while the soft gluons form a highly occupied coherent non-Abelian gauge field. The collision of two such infinitely thin condensates produces the Glasma whose evolution can be described classically by solving the Yang-Mills equations for early proper times. At finite nuclei width, the collision region is not pointlike anymore, and the nuclei, the collision, and the evolution of the Glasma can not be described separately and require one consistent simulation that covers all these steps.
A suitable numerical method was developed in our previous publication [34] based on the colored particlein-cell method (CPIC) [35][36][37][38], which is a non-Abelian extension of the particle-in-cell method for the simulation of Abelian plasmas [39,40]. In contrast to the traditional approach of simulating the Glasma, where the field equations are solved in the forward light-cone parametrized by proper time τ and space-time rapidity η s , we describe the collision in the laboratory frame using the lab-frame time t and the longitudinal beam direction z. The most striking difference of this approach is the explicit inclusion of the nuclei in the simulation, whereas in a boost-invariant simulation the information about the nuclei and their color currents is completely encoded in the initial conditions at the boundary of the light-cone, i.e. τ = 0. To solve this problem numerically we simulate the continuous color charge densities of the nuclei with a large number of color-charged pointlike particles, mimicking the dynamics of the continuous cloud of color charges on a lattice. This enables us to describe the full 3+1 dimensional collision and the subsequent evolution of the Glasma beyond the boostinvariant approximation. For a more detailed description we refer the reader to [34].
The initial conditions in our simulation differ from the traditional approach as well. Instead of starting at τ = 0, our simulation begins before the collision with the nuclei well-separated in the longitudinal direction. Here we quickly review how to solve the Yang-Mills equations in the covariant gauge and the transformation to the temporal gauge in the laboratory frame for a single nucleus. We base our model of the initial state on the McLerran-Venugopalan (MV) model [41,42], extended by a thickness parameter in longitudinal direction. The transverse charge density ρ a (x T ) as a function of the transverse coordinate x T is a random variable following the usual gauge-invariant Gaussian probability functional W [ρ] with the two-point correlation function where µ is the MV model parameter controlling average color charge density and g is the Yang-Mills coupling constant. For a nucleus moving in the positive z direction, we embed this two-dimensional charge density into the three-dimensional laboratory frame via where we introduce the thickness parameter L. In the limit of L → 0 we have f (x − ) ∝ δ(x − ) and restore the boost-invariant limit of the original MV model. Note that this model explicitly neglects non-trivial longitudinal color structure [43]. The only non-vanishing component of the light-like color current of the nucleus is then given by where t a are the generators of the gauge group SU(N ). The subscript "cov" denotes that this defines the color current in the covariant gauge ∂ µ A µ,a = 0. Using this ansatz we can solve the Yang-Mills equations in the covariant gauge by finding a solution to the twodimensional Poisson equation which is solved by whereρ a (k T ) is the Fourier transform of ρ a (x T ). We introduced an infrared regulator m and an ultraviolet cutoff Λ, since the MV model is both infrared and UV divergent. The regularization in (6) should be read as a modification of the charge densities ρ a (x T ) while the field equations remain unchanged. Our numerical method requires the gauge fields to satisfy the temporal gauge condition A a 0 = 0. Switching to this gauge from the covariant gauge renders the fields purely transverse. The transverse field components and the color current are given by Having set up the initial conditions the system can be evolved forward in time via the equations of motion on the lattice.
In order to fix the other parameters of the initial conditions for a given collision energy √ s N N we take the following approach: the longitudinal thickness L introduced in our implementation of the MV model somehow needs to be related to the Lorentz factor γ. It seems natural that L should be proportional to the Lorentz contracted length of the nucleus γ −1 R, where R is the nuclear radius. One possibility is to define [34] γ = R 2L .
This way L is fixed by the geometry of the Lorentz-contracted nucleus. We determine the saturation momentum Q s using the estimation Q 2 s ≈ √ s N N 0.25 GeV 2 with √ s N N given in GeV [44][45][46].
The coupling constant g is set by the one-loop beta function at the saturation scale Q s , which gives values close to g ≈ 2.
The relation between the MV model parameter µ and Q s is non-trivial [45] and for simplicity we choose 0.75 g 2 µ Q s as suggested in [7]. The Lorentz gamma factor is given by γ =  (3), which should give reasonable qualitative results [47]. During the simulation we record the components of the energy-momentum tensor T µν in the laboratory frame and average over a number of collision events to obtain the expectation value T µν . In the MV model most of the T µν components vanish after averaging over all initial conditions due to homogeneity and isotropy in the transverse plane. Therefore the energy-momentum tensor reduces to where ε is the energy density, p L and p T are the longitudinal and transverse pressure components and S L is the longitudinal component of the Poynting vector. The local rest frame energy density is obtained by diagonalizing the energy-momentum tensor T µ ν : Given the electric and magnetic fields E a i and B a i we have Using proper time τ = √ t 2 − z 2 and space-time rapidity η s = 1 2 ln [(t − z) / (t + z)] we can plot the profile of the energy density as a function of η s for various values of τ . Due to the extended collision region there is some ambiguity in setting the coordinate origin for the transformation to the comoving frame. We choose the spacetime coordinate of the maximum of p T (t, z) , which is generated by the initially purely longitudinal Glasma fields. This coordinate origin is slightly later than the space-time coordinates defined by the maximum overlap of the nuclei.
In Figure 2 we see a calculation of the rapidity profile of the local rest frame energy density for a RHIC-like collision: choosing parameters √ s N N = 200 GeV and m = 0.2 GeV, we obtain an approximate Gaussian profile of the local rest frame energy density, in the range η s ∈ (−1, 1). Using a Gaussian fit we compute σ η and extrapolate to higher η s . At higher values of the IR regulator m the rapidity profiles become more narrow, which shows that the transverse size of the correlated color structures ∼ m −1 has an effect on broken boost invariance. In the CGC literature it is well-known that quantities such as the initial energy density or the gluon multiplicity can be sensitive to the choice of the infrared regulator [43,45,48,49]. Here we add another example of a strong dependency on the infrared regulation. Even though the results of our model should be regarded as more qualitative than quantitative, we compare our findings to experimental results: it is an interesting observation that the rapidity profile of the energy density for m = 0.2 GeV agrees with the measured rapidity profile of pion multiplicities for the most central collisions at RHIC. Of course, a direct comparison of ε loc and dN/dy profiles is not strictly valid: the gluon number distribution can be somewhat broader than the energy density [50]. The correct approach would be to use our results as initial conditions for hydrodynamic simulations in order to make a more direct connection with measured observables. The subsequent hydrodynamic expansion of the system likely increases the width of the profiles further as mentioned in [29], which would favor the more narrow curves (b) and (c) in Figure 2. Under these assumptions, the width of the rapidity profile of measured charged particle multiplicities can be seen as an upper limit for realistic rapidity profiles computed from our simulation. We also compare to the Gaussian rapidity profile of the hydrodynamic Landau model [15] with σ Landau = √ ln γ. The profiles have been computed at τ 0 = 1 fm/c, which roughly corresponds to the transition from the Glasma to the QGP. In general we observe that the rapidity profiles quickly converge for τ 0.3 fm/c and afterwards become independent of τ . We also observe free-streaming behavior signaled by ε loc ∝ 1/τ for a wide range of η s and longitudinal velocities v z ∼ z t . This implies that there is only negligible flow of energy between different rapidity directions. Our results suggest that the rapidity dependence of the Glasma is fixed early on in the collision and remains unchanged thereafter. To see the effects of increased longitudinal thickness, we repeat the same calculation for √ s N N = 130 GeV in Figure 3. We fit to a Gaussian shape and find that, as expected, the profiles become more narrow as compared to the RHIC-like case. This time, our results are also more narrow than Landau's prediction.
Investigating the reason behind the non-boostinvariant creation of the Glasma we look at the spacetime distribution of the transverse pressure p T (z, t) in and along the forward light-cone in Figure 4. As can be seen from Eq. (15), the transverse pressure is solely due to the presence of longitudinal fields, i.e. p T (z, t) is equivalent to the energy density generated by longitudinal fields ε L (z, t) . In the boost-invariant case the longitudinal fields would be constant along the boundary of the light-cone as determined by the boost-invariant Glasma initial conditions at τ = 0. In our simulations we see very different behavior: the longitudinal fields are mostly centered around the maximum in the extended collision region t ∼ z ∼ 0 and decrease rather quickly along the t = ±z boundaries. Since the initially longitudinal fields are the starting point of the evolution of the Glasma, this sharp decrease means that for some fixed proper time τ there is less Glasma at higher values of rapidity η s , leading to the observed Gaussian profiles of the energy density. Furthermore, we see a mostly constant spatial distribution of p T (z, t) in Figure 4 for later times t. A similar distribution has been found in holographic models of heavy-ion collisions [17,19].
The results in Figures 2, 3 and 4 have been computed on a lattice with 2048 cells in the longitudinal direction with a length of 6 fm and 192 2 cells for the transverse plane with an area of (6 fm) 2 and a statistical average over 15 events. The tetragonal lattice (with much smaller longitudinal than transverse lattice spacing) is more suited for describing the collision of Lorentz-contracted nuclei. By varying the lattice spacings, the results were checked for discretization errors. Although accessing wider rapidity ranges is possible in principle, we restricted the results to η s ∈ (−1, 1) due to numerical issues: at high space-time rapidity η s near the boundary of the light cone, the computation of ε loc via Eq. (13) becomes very sensitive to cancellations in the square root term. Furthermore, the fields of the nuclei, being proportional to the longitudinal profile f (x − ), become increasingly larger compared to the Glasma fields as L → 0, which converge to the finite, boost-invariant result at mid-rapidity for high √ s N N . This leads to an unwanted modification of ε loc by large transverse fields of the nuclei at larger η s . We cannot strictly separate the Glasma from the nucleus fields, which is a clear disadvantage to simulations performed strictly in the forward light cone [29]. For these reasons we currently can not obtain reliable results for collision energies present at the LHC and further improvements in the numerical scheme are needed. Concluding, we found that Gaussian rapidity profiles in energy density can arise from CGC collisions of nuclei with finite longitudinal thickness in classical 3+1 dimensional Yang-Mills simulations. The width of the profiles is controlled by the energy of the incoming nuclei, but depends also crucially on the infrared modes of the fields. Presumably the infrared dependencies of our results could be fixed by more realistic initial conditions obtained from a JIMWLK evolution. It would be interesting to understand the mechanism behind the creation of Glasma fields in a non-boost-invariant setting with finite nucleus thickness, and in particular also the connection to infrared regulation. In any case, our result shows that there is a mechanism that breaks boost invariance at the classical level in the Glasma related to the finite thickness of the colliding nuclei, and further investigations are needed to see how this effect compares to other previously studied approaches like the rapidity dependence of JIMWLK. It is also exciting to see that by lifting the assumption of boost invariance, we have shown that our weak coupling results are in qualitative agreement with the strong coupling results exhibiting a Gaussian rapidity profile of the rest frame energy density and similar transverse pressure distributions in the laboratory frame [17]. That is to say, our result could be interpreted to indicate that the difference between the previous weak and strong coupling simulations is due to the initial conditions, not the weak or strong coupling dynamics.
This work has been supported by the Austrian Science Fund FWF, Project No. P 26582-N27 and Doctoral program No. W1252-N27. The computational results presented have been achieved using the Vienna Scientific Cluster. We would like to thank Aleksi Kurkela for useful discussions.