Broken boost invariance in the Glasma via finite nuclei thickness

TL;DR

This paper simulates the early non-boost-invariant Glasma in heavy ion collisions, revealing deviations from boost invariance and matching experimental data by extending the McLerran-Venugopalan model to finite nucleus thickness.

Contribution

It introduces a finite nucleus width parameter into the McLerran-Venugopalan model to study non-boost-invariant Glasma evolution.

Findings

Rapidity profile deviations from boost invariance

Ability to produce broad and narrow energy density profiles

Surprising agreement with RHIC experimental data

Abstract

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.

Equations

$$⟨{\rho }^a(x_T){\rho }^b(y_T)⟩=g^2{\mu }^2{\delta }^{ab}{\delta }^{(2)}(x_T-y_T),$$

$$f(x^{-})=\frac{1}{\sqrt{2\pi }L}\exp \left(-(x^{-}{)}^2/L^2\right),$$

$$J_{\text{cov}}^{+}(x_T,x^{-})$$

$$D_{\mu }{F}^{\mu \nu }=J^{\nu },$$

$$-{\Delta }_T{A}^{+}(x_T,x^{-})=J_{\text{cov}}^{+}(x_T,x^{-}),$$

$${\varphi }^a(x_T)$$

$$A^{+}(x_T,x^{-})$$

$$A_i(x_T,x^{-})=\frac{1}{ig}V(x_T,x^{-}){\partial }_i{V}^{†}(x_T,x^{-}),$$

$$J^{+}(x_T,x^{-})=V(x_T,x^{-})J_{\text{cov}}^{+}(x_T,x^{-})V^{†}(x_T,x^{-}),$$

$$V(x_T,x^{-})$$

$$\gamma =\frac{R}{2L}.$$

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

$$⟨{\epsilon }_{\text{loc}}⟩$$

$$⟨\epsilon ⟩$$

$$⟨p_T⟩$$

$$⟨p_L⟩$$

$$⟨S_L⟩$$

$${\epsilon }_{\text{loc}}({\tau }_0,{\eta }_s)\approx {\epsilon }_{\text{loc}}({\tau }_0,0)\exp \left(-\frac{{\eta }_s^2}{2{\sigma }_{\eta }^2}\right),$$

Full text

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 (Heinz and Snellings, 2013; Gale et al., 2013a). A first principles description of the initial state of heavy-ion collisions is provided by the Color Glass Condensate (CGC) framework (Gelis et al., 2010; Gelis, 2013, 2015). The CGC is a classical effective field theory for nuclear matter at ultrarelativistic energies. Models such as the IP-Glasma (Schenke et al., 2012a, b) in combination with hydrodynamics are able to correctly reproduce azimuthal anisotropies and event-by-event multiplicity distributions (Gale et al., 2013b; Snellings, 2011). Furthermore, the CGC can explain long-range rapidity correlations like the ridge (Dumitru et al., 2008; Iancu, 2014).

A Gaussian shaped rapidity profile of particle multiplicity can be found in experiments covering various energy ranges, from LHC (Abbas et al., 2013) to RHIC Beam Energy Scan (Bearden et al., 2005; Flores, 2016). This shape is well explained by the Landau model (Landau, 1953) up to RHIC energies (Bearden et al., 2005), which assumes full stopping of the colliding nuclei. The Landau model is in contrast to the Bjorken model (Bjorken, 1983) which relies on approximate boost invariance. A Gaussian profile has also been found in holographic calculations of colliding shock waves (Casalderrey-Solana et al., 2013; van der Schee, 2014; van der Schee and Schenke, 2015).

In its original formulation collisions in the CGC picture are assumed to be boost-invariant (Kovner et al., 1995a, b; Krasnitz and Venugopalan, 1999; Lappi, 2011) 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 (Jalilian-Marian et al., 1997; Iancu et al., 2001; Mueller, 2001; Ferreiro et al., 2002), 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 (Schenke et al., 2014; Schenke and Schlichting, 2016). 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 at the classical level (Özönder and Fries, 2014). In the case of proton-nucleus collisions methods have been developed to systematically include finite width corrections of the nucleus (Altinoluk et al., 2014, 2016; Altinoluk and Dumitru, 2016). 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 boost-invariant 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 (Gelfand et al., 2016) based on the colored particle-in-cell method (CPIC) (Dumitru and Nara, 2006; Moore et al., 1998; Hu and Müller, 1997; Dumitru et al., 2007), which is a non-Abelian extension of the particle-in-cell method for the simulation of Abelian plasmas (Esirkepov, 2001; Verboncoeur, 2005). 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 $\tau$ and space-time rapidity $\eta_{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. $\tau=0$. To solve this problem numerically we simulate the continuous color charge densities of the nuclei with a large number of color-charged point-like 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 boost-invariant approximation. For a more detailed description we refer the reader to (Gelfand et al., 2016).

The initial conditions in our simulation differ from the traditional approach as well. Instead of starting at $\tau=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 (McLerran and Venugopalan, 1994, 1994), extended by a thickness parameter in longitudinal direction. The transverse charge density $\rho^{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[\rho]$ with the two-point correlation function

[TABLE]

where $\mu$ 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 $\rho^{a}(x_{T},x^{-})=f(x^{-})\rho^{a}(x_{T})$ with a longitudinal profile function $f(x^{-})$, where $x^{\pm}\equiv\left(t\pm z\right)/\sqrt{2}$ are the usual light-cone coordinates. For the longitudinal profile we choose a Gaussian

[TABLE]

where we introduce the thickness parameter $L$. In the limit of $L\rightarrow 0$ we have $f(x^{-})\propto\delta(x^{-})$ and restore the boost-invariant limit of the original MV model. Note that this model explicitly neglects non-trivial longitudinal color structure (Fukushima, 2008). The only non-vanishing component of the light-like color current of the nucleus is then given by

[TABLE]

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 $\partial_{\mu}A^{\mu,a}=0$. Using this ansatz we can solve the Yang-Mills equations

[TABLE]

in the covariant gauge by finding a solution to the two-dimensional Poisson equation

[TABLE]

which is solved by

[TABLE]

where $\tilde{\rho}^{a}(k_{T})$ is the Fourier transform of $\rho^{a}(x_{T})$. We introduced an infrared regulator $m$ and an ultraviolet cutoff $\Lambda$, since the MV model is both infrared and UV divergent. The regularization in (6) should be read as a modification of the charge densities $\rho^{a}(x_{T})$ while the field equations remain unchanged.

Our numerical method requires the gauge fields to satisfy the temporal gauge condition $A_{0}^{a}=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

[TABLE]

[TABLE]

with the temporal Wilson line

[TABLE]

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 $\sqrt{s_{NN}}$ 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 $\gamma$. It seems natural that $L$ should be proportional to the Lorentz contracted length of the nucleus $\gamma^{-1}R$, where $R$ is the nuclear radius. One possibility is to define (Gelfand et al., 2016)

[TABLE]

This way $L$ is fixed by the geometry of the Lorentz-contracted nucleus. We determine the saturation momentum $Q_{s}$ using the estimation $Q_{s}^{2}\approx\left(\sqrt{s_{NN}}\right)^{0.25}\,\mbox{GeV}^{2}$ with $\sqrt{s_{NN}}$ given in GeV (Lappi, 2006, 2008; Kharzeev and Levin, 2001). The coupling constant $g$ is set by the one-loop beta function at the saturation scale $Q_{s}$, which gives values close to $g\approx 2$. The relation between the MV model parameter $\mu$ and $Q_{s}$ is non-trivial (Lappi, 2008) and for simplicity we choose $0.75\,g^{2}\mu\simeq Q_{s}$ as suggested in (Schenke et al., 2012b). The Lorentz gamma factor is given by $\gamma=\sqrt{s_{NN}}/(2m_{N})$ with the nucleon mass $m_{N}\approx 1\>\mbox{GeV}$ and consequently we can find $L$ via Eq. (11), setting $R$ to the radius of a Gold nucleus. Since our results depend on the IR regulator $m$ we try out different values: $m=0.2\,\mbox{GeV}$, $m=0.4\,\mbox{GeV}$ and $m=0.8\,\mbox{GeV}$. The UV modes are regulated by $\Lambda=10\>\mbox{GeV}$. The simulation is performed with the gauge group $\mbox{SU}(2)$ instead of SU$(3)$, which should give reasonable qualitative results (Ipp et al., 2011).

During the simulation we record the components of the energy-momentum tensor $T^{\mu\nu}$ in the laboratory frame and average over a number of collision events to obtain the expectation value $\left\langle T^{\mu\nu}\right\rangle$. In the MV model most of the $T^{\mu\nu}$ components vanish after averaging over all initial conditions due to homogeneity and isotropy in the transverse plane. Therefore the energy-momentum tensor reduces to

[TABLE]

where $\left\langle\varepsilon\right\rangle$ is the energy density, $\left\langle p_{L}\right\rangle$ and $\left\langle p_{T}\right\rangle$ are the longitudinal and transverse pressure components and $\left\langle S_{L}\right\rangle$ is the longitudinal component of the Poynting vector. The local rest frame energy density is obtained by diagonalizing the energy-momentum tensor $\left\langle T_{\,\,\nu}^{\mu}\right\rangle$:

[TABLE]

Given the electric and magnetic fields $E_{i}^{a}$ and $B_{i}^{a}$ we have

[TABLE]

Using proper time $\tau=\sqrt{t^{2}-z^{2}}$ and space-time rapidity $\eta_{s}=\frac{1}{2}\ln\left[\left(t-z\right)/\left(t+z\right)\right]$ we can plot the profile of the energy density as a function of $\eta_{s}$ for various values of $\tau$. 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 space-time coordinate of the maximum of $\left\langle p_{T}(t,z)\right\rangle$, 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 $\sqrt{s_{NN}}=200\,\mbox{GeV}$ and $m=0.2\,\text{GeV}$, we obtain an approximate Gaussian profile of the local rest frame energy density,

[TABLE]

in the range $\eta_{s}\in(-1,1)$. Using a Gaussian fit we compute $\sigma_{\eta}$ and extrapolate to higher $\eta_{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 $\sim 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 (Lappi, 2008; Fukushima, 2008; Fujii et al., 2009; Fukushima and Gelis, 2012). 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\,\text{GeV}$ agrees with the measured rapidity profile of pion multiplicities for the most central collisions at RHIC. Of course, a direct comparison of $\varepsilon_{\text{loc}}$ and $dN/dy$ profiles is not strictly valid: the gluon number distribution can be somewhat broader than the energy density (Hirano and Nara, 2004). 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 (Schenke and Schlichting, 2016), 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 (Landau, 1953) with $\sigma_{\text{Landau}}=\sqrt{\ln\gamma}$.

The profiles have been computed at $\tau_{0}=1\,\mbox{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 $\tau\gtrsim 0.3\,\mbox{fm/c}$ and afterwards become independent of $\tau$. We also observe free-streaming behavior signaled by $\varepsilon_{\text{loc}}\propto 1/\tau$ for a wide range of $\eta_{s}$ and longitudinal velocities $v_{z}\sim\frac{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 $\sqrt{s_{NN}}=130\,\mbox{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-boost-invariant creation of the Glasma we look at the space-time distribution of the transverse pressure $\left\langle p_{T}(z,t)\right\rangle$ 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. $\left\langle p_{T}(z,t)\right\rangle$ is equivalent to the energy density generated by longitudinal fields $\left\langle\varepsilon_{L}(z,t)\right\rangle$. 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 $\tau=0$. In our simulations we see very different behavior: the longitudinal fields are mostly centered around the maximum in the extended collision region $t\sim z\sim 0$ and decrease rather quickly along the $t=\pm 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 $\tau$ there is less Glasma at higher values of rapidity $\eta_{s}$, leading to the observed Gaussian profiles of the energy density. Furthermore, we see a mostly constant spatial distribution of $\left\langle p_{T}(z,t)\right\rangle$ in Figure 4 for later times $t$. A similar distribution has been found in holographic models of heavy-ion collisions (Casalderrey-Solana et al., 2013; van der Schee and Schenke, 2015).

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\,\mbox{fm}$ and $192^{2}$ cells for the transverse plane with an area of $\left(6\,\mbox{fm}\right)^{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 $\eta_{s}\in(-1,1)$ due to numerical issues: at high space-time rapidity $\eta_{s}$ near the boundary of the light cone, the computation of $\varepsilon_{\text{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\rightarrow 0$, which converge to the finite, boost-invariant result at mid-rapidity for high $\sqrt{s_{NN}}$. This leads to an unwanted modification of $\varepsilon_{\text{loc}}$ by large transverse fields of the nuclei at larger $\eta_{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 (Schenke and Schlichting, 2016). 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 (Casalderrey-Solana et al., 2013). 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.