3-D structure of the Glasma initial state -- Breaking boost-invariance by collisions of extended shock waves in classical Yang-Mills theory

We simulate the 3+1 D classical Yang-Mills dynamics of the collisions of longitudinally extended nuclei, described by eikonal color charges in the Color Glass Condensate framework. By varying the longitudinal thickness of the colliding nuclei, we discuss the violations of boost invariance and explore how the boost invariant high-energy limit is approached. Subsequently, we develop a more realistic model of the three dimensional color charge distributions that connects longitudinal and transverse fluctuations to the $x$ and $k_\bot$ dependence of transverse momentum dependent parton distributions, and explore the resulting rapidity profiles and their longitudinal fluctuations.


I. INTRODUCTION
High-energy heavy-ion collisions at the Large Hadron Collider (LHC) and the Relativistic Heavy Ion Collider (RHIC) produce a deconfined Quark Gluon Plasma (QGP), whose space-time dynamics is well described by relativistic viscous hydrodynamics [1,2]. Over the course of the last few years sophisticated multi-stage models, including the early time pre-equilibrium dynamics, the intermediate hydrodynamic QGP phase, as well as the late stage hadronic re-scattering have been developed to describe the space-time dynamics of heavy-ion collisions. Nevertheless, a comprehensive understanding hinges to a considerable extent on a proper theoretical understanding of the initial conditions, characterizing the non-equilibrium QGP created immediately after the collision of heavy nuclei.
The Color Glass Condensate (CGC) effective theory of high-energy QCD [3,4] provides a framework to study the earliest stages of the collision which encompasses the description of dense colliding nuclei prior to the collision and their initial energy deposition [5,6]. Within this framework, the high-energy boost-invariant limit has been studied in great detail [6][7][8][9] and lead to the successful development of microscopic initial state models, such as IP-Glasma [10,11] and MC-KLN [12,13], which describe the two dimensional transverse energy density profiles.
Based on new experimental measurements at lower energies and interesting results of experimental measurements of longitudinal fluctuations at high energies, there has recently been a renewed interest in understanding the three dimensional structure of the initial state, that has lead to the development of various three dimensional initial state models [14][15][16][17]. However, so far only a limited number of first principles insights have been obtained beyond the high-energy boost invariant limit, with existing analytic studies focusing on the fragmentation region [18][19][20][21][22] and sub-eikonal corrections to particle production in momentum space [23,24]. Beyond that, there is a number of proposals that aim to generalize the boostinvariant IP-Glasma model to 3+1D, by taking into account the small x JIMWLK evolution of the color charge * pragya@physik.uni-bielefeld.de distributions [25,26], or by including sub-eikonal corrections due to the finite longitudinal thickness of the colliding nuclei based on colored particle-in-cell methods (CPIC) [27][28][29].
In this work we develop a framework to perform 3+1D classical Yang-Mills simulations of the initial energy deposition in heavy-ion collisions, which as in [27][28][29] take into account the finite longitudinal extent of the colliding nuclei. While the colliding nuclei are still propagating on the light-cone, the collision region no longer remains a point but stretches to a diamond like region in spacetime, such that the formation of the Glasma begins as soon as the two nuclei begins to overlaps. Similar to previous studies [27][28][29] this also entails the explicit propagation of the color currents of the colliding nuclei, which is in contrast to the boost-invariant case where they are considered to be static. Within a simple model of the color charge distribution of each nucleus, we perform a detailed investigation of the dynamics during and shortly after the collision as a function of the longitudinal thickness of the colliding nuclei, and contrast our results with the high-energy limit of infinitely thin shocks. Subsequently, we develop a more physical model that connects the color charge distributions in the colliding nuclei to parton distributions inside the nuclei, and we discuss the rapidity profiles and fluctuations that emerge within this model. This work is organized as follows: We first introduce the general formalism to discretize 3+1D Yang-Mills equation in II, followed by the initial condition for gauge fields and currents. In Section III, we present our results for 3+1D collision of individual color charges based on MV-model, followed by exploring rapidity dependence of the observable. We present our results for collision with realistic charge distributions in Section IV, where we investigate the fluctuation at lower energies.

II. GENERAL FORMALISM FOR 3+1D COLLISION IN YANG-MILLS THEORY
Based on the Color Glass Condensate (CGC) effective description of high-energy QCD, the initial state energy deposition and early time dynamics in high-energy heavy-ion collisions can be described semi-classically by solving classical Yang-Mills equations of motion for the gluon fields A µ in the presence of fluctuating color charges ρ, which characterize the nuclear parton content. Even though a complete analytical treatment of the Yang-Mills equation is not possible, it is remarkable that the initial state immediately after the collision (τ = 0 + ) can be determined analytically for boost-invariant collisions in the high-energy limit [9,30]. Beyond the time of the collision τ = 0, where a far-from-equilibrium Glasma is produced [6,31], the classical Yang-Mills dynamics becomes highly non-linear, and additional analytic insights can only be obtained in the limit where one or both of the two sources are considered to be weak such that the equations of motion linearize.
Nevertheless, important insights into the early time non-equilibrium dynamics have been established based on numerical simulations [8,10,32,33], where in order to describe the non-linear dynamics of the boost invariant Glasma, the effectively 2+1 dimensional classical Yang-Mills equations are discretized on a lattice and solved numerically in the forward light cone i.e. for τ > 0. Below we explain how to generalize this setup to simulate the collision of nuclei with a finite longitudinal thickness in 3+1D collisions, where in contrast to the boost invariant high-energy limit, the entire space-time dynamics of the collision has to be simulated numerically, including the explicit evolution of color charges before, during and after the collision.

A. Discretization of gauge fields and currents
We follow standard procedure in the context of classical-statistical lattice gauge theory simulations [5,7] and start by discretizing the classical Yang-Mills Hamiltonian in temporal axial A t = 0 gauge on a three dimensional N x × N y × N z lattice with lattice spacing a µ in thê µ direction [6,8] where x denotes the lattice site, i, j = x, y, z are spatial Lorentz indices, a = 1 · · · N 2 c − 1 is the color index and g µν = (+1, −1, −1, −1) denotes the Minkowski metric. We note that within the Hamiltonian formulation in Eq. 1, the electric field strength is represented in terms of the lattice electric field variables while the magnetic field strength is given in terms of the lattice plaquette variables which are formed of SU(N c ) group valued lattice gauge links and the pre-factors in the definitions in Eqns. 2, 3 have been arranged such that all lattice variables are explicitly dimensionless. Within the CGC framework the incoming nuclei are described by eikonal color currents J µ R/L (x) which propagate along the light cones and provide a source for the classical gluon fields, where the subscripts R and L represent the nuclei coming in from the right and left respectively. Before the collision the initial conditions for the currents J µ R/L (x) in covariant (∂ µ A µ = 0) gauge are given in terms of the color charge densities ρ a where t a are generators in the fundamental representation. However, during the collision both currents J µ R/L (x) will receive a color-rotation, which in the 3+1 D setup has to be calculated by solving dynamical equations of motions for the currents J µ R/L (x). We therefore also discretize the color currents on the lattice, where by keeping track of the relevant light-cone (±) components, the currents are defined as with "dynamical" (dyn) and "static" (stat) currents on alternating half-integer time slices, as usual in a leap-frog scheme.

B. Evolution equations & Gauss Law
By performing the variation of the lattice Hamiltonian w.r.t to electric fields and gauge fields, one obtains the Hamiltonian equations of motion for the lattice gauge link and electric field variables. We employ a leap frog algorithm with time step a t = 0.8 × min(a z , a ⊥ ), where gauge links are defined at every full time step whereas the electric fields are calculated for every half-integer time step, such that the update rule for the lattice gauge links takes the form whereas for the evolution of the lattice electric fields one also has to take into account the coupling to the eikonal currents, such that the update rule for the lattice electric fields is given by Due to the explicit appearance of the currents on the rhs of Eq. 10, the color currents J ± x (t) also have to be treated as dynamical degree of freedom as -in contrast to the boost-invariant high-energy limit -they are present not only on the infinitesimal boundary of the light-cone but throughout the entire simulation volume and thus affect the evolution of the classical Yang-Mills fields. Since in the eikonal limit the different components of the current are related by J 0 R/L (x) = ∓J z R/L (x), it is straightforward to construct the dynamical equation of motion of the currents from the (covariant) charge conservation equation D µ J µ = 0 as ∂ 0 J ± (x) = ∓D z J ± (x), and we employ the following update rules where D F/B µ denote the forward and backward covariant derivatives We note that due to the leap-frog discretization in Eqns. 11, the dynamical currents do not propagate exactly at the speed of light, but instead satisfy the same lattice dispersion relation as the lattice gauge fields, which converges to a light-like dispersion in the continuum limit. This is different from CPIC method where the color charges have light-like dispersion and the dispersion of the gauge field depends upon the details of the numerical scheme [34]. However, as in CPIC, the lattice version of the Gauss law constraint is automatically satisfied at each time, as long as it is satisfied by the initial conditions, as can be checked straightforwardly by evaluating the time derivative of Eq. 13, based on the equations of motion for the lattice gauge links, electric fields and currents as in Eqns. 8-11.

C. Initial Conditions for 3+1D collisions
Since for 3+1 D collisions of extended nuclei, one has to simulate the dynamics of the color charges before, during and after the collision, the initial conditions for the above evolution equations have be formulated at Minkowski time t 0 < 0 before the collision, where the colliding nuclei are well separated from each other. Since the color charges inside the two nuclei do not interact with each other before the wave-packets overlap, the initial conditions are then determined by the superposition of the analytic solutions for the gauge fields in the presence of the individual color charges of the two nuclei, as illustrated in Fig. 1.
Specifically, in the covariant ∂ µ A µ = 0 gauge, the solution to classical Yang-Mills equations before the collision takes the form where we have explicitly assumed that -at the ini- We note that in ∂ µ A µ = 0 gauge the gauge potentials A ± , only have support in the vicinity of the two lightcones, where color charges are present, as seen in top panel of Fig. 1. However, for real time lattice simulation it is convenient to employ the temporal axial (A 0 = 0) gauge condition, and the corresponding initial conditions can be obtained by performing a gauge transformation, which eliminates the gauge potentials A ± prior to the collision. By following previous works [35][36][37], the corresponding gauge transformation can be expressed as where the color rotations V L/R associated with the left and right moving nuclei are determined by such that V L/R are given by the light-like Wilson lines By performing the corresponding gauge transformation, the initial conditions in temporal axial gauge then take the form which is illustrated on the right panel of Fig. 1 [8,35,36], as indicated in Fig. 1.
So far we have discussed the structure of the initial conditions in the continuum, and we will now address the corresponding lattice implementation. Starting from a given distribution of color charges ρ L/R (t 0 , z, x ⊥ ) at initial time t 0 < 0 discretized on a spatial x, y, z grid, we first compute the covariant gauge A ± cov. (t 0 , z, x ⊥ ) fields according to Eq. 14, and subsequently construct the discretized version V L/R (t 0 , z, x ⊥ ) of the light-like Wilson lines. Since before the collision, the Wilson lines V R (x + , x ⊥ ) are independent of x − , the light-like Wilson lines V R (t 0 , z, x ⊥ ) can be defined to end on the lattice points x, y, z at initial time t 0 as illustrated in the left panel of Fig. 2.
If we consider the left moving nucleus, which is initially located on the right hand side of the lattice, the corresponding Wilson line V R (t 0 , z, x ⊥ ) is equal to the identity for all points z which at the time t 0 < 0 are located to the left of the incoming color charges. Starting from V R (t 0 , z = 0, x ⊥ ) = 1 at the left boundary of the lattice (z i = 0), the Wilson lines V R (t 0 , z, x ⊥ ) for z > 0 can be constructed successively based on the relation [38] Exploiting again the invariance of the A − cov. (x + , x ⊥ ) gauge fields under shifts along the x − direction, the ad-ditional color rotation can then be approximated as as indicated in the right panel of Fig. 2. Based on this procedure, the lattice discretized version of the Wilson lines for the right and left sitting nucleus is then given by By applying the gauge transformation in Eq. 4 to the lattice gauge links, the initial conditions for the lattice gauge links are then determined as Next, in order to initialize the lattice electric fields E x,i t + dt 2 , we make use of the update rule for the gauge links in Eq. 8, to express which relates the electric fields to the gauge links at times t 0 and t 0 + a t . By following the same procedure as outlined above, the gauge links U x,i (t + a t ) are constructed from the color charges propagated by a single time step according to such that the evolution of the color charges in the initialization step satisfies the lattice dispersion. While the above procedure provide initial conditions for the lattice gauge links and electric fields in temporal axial gauge, the color charges ρ L/R (t 0 , z, x ⊥ ) are still given in covariant gauge. Instead of performing a gauge transformation of the charges, we exploit Gauss Law to determine the static currents as where the factor of √ 2 comes from the transformation between Minkowski and light-cone coordinates. Subsequently, the initial value of the dynamical currents J dyn is set by performing half a time step of evolution as Starting from the initial conditions outlined above, we simulate the dynamics of the collision in Minkowksi coordinates x µ = (t, x, y, z) by solving the classical Yang-Mills equations for the lattice gauge links and electric fields, along with the evolution equations for the eikonal currents. While in principle the evolution can be performed up to arbitrary late times, in practice the incoming nuclei will approach to boundary of the lattice at some finite time after the collisions. Since we implement fixed boundary conditions along the longitudinal direction 1 the color charges can not pass through the boundary and we have to stop our simulations before this occurs.
When investigating the initial energy deposition during the collision and early-time dynamics of the Glasma, we will primarily focus on the evolution of energy momentum tensor T µν (x), which we compute as , with spatial Lorentz indices i, j and E loc (x) and B loc (x) are local electric and magnetic fields calculated using smeared operator definitions, such that they are defined at same position i.e (x + a i /2 + a t /2) as Besides the space-time evolution of the energy momentum tensor, we will also consider the evolution of the field intensity of longitudinal ( ) and transverse (⊥) components of the (chromo-) electromagnetic fields Since the focus of our investigation is to understand the time evolution and the longitudinal structure of the fields, we will frequently perform averages over the transverse, which we denote as Due to the computational complexity of 3 + 1 D simulations, we will perform all of our simulations for the SU (2) gauge group rather than the physical gauge group SU (3) of QCD. Even though we do not expect to see any qualitative changes between the SU (2) and SU (3) dynamics, quantitative values of observables will change, and should therefore not be compared directly to experimental results.

III. 3+1D COLLISIONS OF INDIVIDUAL COLOR CHARGES
Based on the above simulation framework, we will now study the initial energy deposition and early time dynamics of the Glasma. Before we turn to simulations involving realistic models of the color charge distributions of the colliding nuclei, it proves insightful to first consider the collision of individual ensembles of color charges, to test the framework and develop an intuitive picture of the underlying dynamics.
We follow the previous works [6,39], and sample the transverse distribution of color charges based on the McLerran-Venugopalan (MV) model as where m ∼ Λ QCD is an infrared regulator enforcing color neutrality on large distance scales, Λ is an ultra-violet cut-off suppressing color charge fluctuations on short distance scales, such that the model has a well defined continuum limit and if not stated otherwise we employ m/Q s = 1 and Λ/Q s = 5 in our simulations.
Subsequently, the three dimensional color charge distribution ρ a L/R (x, y, z) = ρ a (2D) L/R (x, y)a z T (z) is obtained by multiplying the transverse color charge distribution with the same Gaussian profile at each point where the dimensionful parameter R controls the longitudinal extend of the nucleus.
Since the initial color charge distributions are characterized in terms of the dimensionful scales Q s and R, the latter can be used to set the scale of the lattice calculation by specifying the value of Q s a ⊥ and R/a z . Generally, the (transverse) lattice spacing has to be chosen sufficiently small to avoid discretization errors Q s a ⊥ 1, while at the same time the transverse simulation volume N ⊥ a ⊥ should be large compared to the color charge correlation length ∼ 1/m. Similarly, the longitudinal color charge distribution has to be smooth on the scale of a single lattice spacing a z R, while at the same time the longitudinal extend of the lattice N z a z has to be sufficiently large to allow for a long enough time evolution after the collision. Within the above setup, we have varied both the lattice spacings Q s a ⊥ , R/a z as well as the lattice length Q s N ⊥ a ⊥ and R/N z a z to check that discretization errors do not play a significant role, and if not stated otherwise we will present results for N ⊥ = 128, N z = 2048 with Q s a ⊥ = 0.125 and R/a z = 16 in the following.
A. Stable propagation of color charges before and after the collision Before we address the dynamics of the collision, we briefly verify that -within our numerical setupthe color charges of the individual nuclei propagate in a stable fashion. We illustrate this behavior in Fig. 3, where the top panel shows the evolution of the longitudinal profiles of the color charge distribution ρ(t, x, y, z) = ρ a (t, x, y, z)ρ a (t, x, y, z) at a randomly chosen point x, y in the transverse plane. By comparing the initial color charge distribution ρ a (t 0 , x, y, z) determined according to Eq. 25, to the charge density ρ a (t 0 , x, y, z) = V ab (t 0 , x, y, z)ρ cov. b (t 0 , x, y, z) reconstructed from the color charge distribution in covariant gauge, we observe an excellent agreement demonstrating that the re-construction of the charge density based on Gauss's law works as expected. By comparing the dynamically evolved charge distribution ρ a (t, x, y, z) to the translated initial conditions ρ a (t 0 , x, y, z − c(t − t 0 )), we can further confirm that for sufficiently small lattice spacing a z R, the numerical dispersion of the currents is small, such that over the relevant time scales the nuclei propagate in a stable fashion at almost the speed the light. Beside the stable propagation of the color charges, it is also important that the gluon fields induced by the color charges propagate in a stable way along side the charges, as can be seen from the bottom panel of Fig. 3, where we show the evolution of the longitudinal profiles of the average electric and magnetic field strengths E 2 B 2 (t, z) ⊥ . By comparing the longitudinal field strength profiles in Fig. 3 at different times, one again concludes the nuclei propagate in a stable fashion over the relevant time scales; moreover the analytic result for the evolution of the field strength is given as: 33) which follows directly from Eqns. 14 and 31, is well reproduced by our real-time lattice simulations, indicating the residual discretization errors are indeed small.

B. Evolution of the fields during and after the collision
Now that we have established the validity of our setup, we will analyze the energy deposition and early time dynamics of the collisions. Before we present our numerical results, we briefly recall the structure of the Glasma fields in the high-energy boost invariant limit [5,9,30] which will serve as a basis for comparison. Before the collision, the incoming nuclei feature the transverse electric and magnetic fields known as the Weiszaecker-Williams fields (WW) localized in narrow strips along the light-cones. Even though the structure of the fields during the collisions is not analytically accessible, it is well established that the initial state immediately after the collision (τ = 0 + ) features boost invariant longitudinal electric and magnetic fields in the forward light-cone  where V (x ⊥ ) is defined in Eq. 15 and we adapted the usual τ, η coordinates, with τ = √ t 2 − z 2 and η = atanh(z/t). Subsequently, for τ > 0 the initial Glasma flux tubes in Eq. 35 begin to expand in transverse space and loose their coherence, such that after a time scale ∼ 1/Q s the Glasma features longitudinal and transverse color fields of comparable magnitude, resulting in a state of approximately vanishing longitudinal pressure p L = T η η 0 [9]. While the evolution of the boostinvariant Glasma produced in the forward light-cone has been explored to detailed extent within numerical simulations [9,10,33] and (semi-)analytic calculations [5], we further note that also in the boost-invariant high-energy limit the eikonal charges receive a color rotation during the collision, and the transverse electric and magnetic fields continue to exist in narrow strips along the lightcones.
Beyond the high-energy boost invariant limit, the formation of the Glasma begins as soon as the color charge distributions of the incoming nuclei start to overlap and persists over an extended period of time until the colliding nuclei have passed through each other. Now in order to analyze the formation of the 3+1D Glasma, we first consider the evolution of the fields at the center of the collision (z = 0), where one has τ = t and η = 0 such that the descriptions in Minkowski (t, z) coordinates and Milne (τ, η) coordinates coincide. We present our results in Fig. 4, where we show the time evolution of the longitudinal and transverse electric and magnetic fields strength during and after the collision for two different values of Q s R = 1/4 and Q s R = 1/16 corresponding to different longitudinal thickness of the colliding nuclei. Since we want to focus on the creation of the Glasma, we have subtracted the field strength associated with Weiszaecker-Williams fields of the colliding charges, i.e. we consider E 2 Glasma (t, z) = E 2 (t, z) − E 2 W W (t, z), and we have defined the origin of the coordinate system such that at t = 0 the charge distributions of the colliding nuclei maximally overlap with each other. During the collision longitudinal electric and magnetic fields build up monotonically as the two nuclei pass through each other, while the transverse electric and magnetic components experience rapid changes. By the time that the incoming nuclei have passed each other, which corresponds to Q s t 0.25 in the top panel and Q s t 1 in the bottom panel, the transverse electric and magnetic fields become small again and longitudinal electric and magnetic fields dominate the energy density. While close to the boost-invariant limit for Q s R = 1/16 the longitudinal electric and magnetic fields have approximately the same magnitude, the longitudinal magnetic field strength is suppressed compared to the longitudinal electric one away from the boost-invariant limit for Q s R = 1/4. Eventually, as the colliding nuclei have passed through each other, transverse electric and magnetic fields are regenerated from their longitudinal counterparts, until at late times the different components become of comparable magnitude. Similar to the boost-invariant case, the different field intensities at the center of the collision (z = η = 0) then start to decay as approximately ∝ 1/τ as indicated by the black dashed lines in Fig. 4 Next we will investigate the dependence of the initial energy deposition and early time dynamics on the longitudinal thickness. In Fig. 5 we present the evolution of transverse magnetic and longitudinal electric fields for different values of Q s R characterizing the longitudinal thickness of the colliding nuclei. 2 Alongside the results from 3+1D numerical simulations, we also show results for the boost invariant limit, obtained by performing 2 + 1D boost-invariant classical Yang-Mills simulations for the same color charge distributions, using the infinite Wilson lines V L/R (x ⊥ ) as an input. Starting from the collision of thick nuclei with Q s R = 1/2, where the collision takes a significant amount of time and the evolution of B 2 ⊥ (t, z = 0) and E 2 (t, z = 0) shows a smooth transition between the different stages, the decrease/increase of the transverse magnetic/longitudinal electric field strength during the collision sharpens significantly as the collision becomes shorter and shorter for smaller values of Q s R. Conversely, the evolution of the fields after the collision for Q s t Q s R is rather insensitive to the longitudinal thickness, and the results for 3+1D collisions smoothly approach the boost-invariant result as Q s R → 0.
So far we have focused on the time evolution in the center of the collision (z = 0), and we will now analyze space-time evolution of longitudinal profiles of the collision in more detail. Instead of showing results for the individual field strength components, we will focus on the evolution of the dominant components T 00 , T 0z , T zz and T ii of the energy momentum tensor, and similar to our previous result subtract the contributions T 00 W W = T zz W W = ±T 0z W W of the Weiszaecker-Williams fields to the energy momentum tensor. By performing three independent simulations, where in the first case we simulate the full collision, while the other two simulations simply propagate of the left/right moving charges, the subtraction takes into account the non-zero dispersion of the color charges due to residual discretization errors, and the energy-momentum tensor of the Glasma T µν Glasma = T µν f ull − T µν W W,L − T µν W W,R vanishes identically before the two nuclei collide.
Our results for the space-time evolution evolution of the energy-momentum tensor are compactly summarized in Fig. 6, where the different panels show the t, z dependence of the various components of T µν Glasma ⊥ in the labframe averaged over the transverse plane. By focusing on the evolution of the transverse pressure (T xx + T yy )/2, one clearly observes the energy deposition in the central region where (T xx + T yy )/2 increases during the collision, exhibits a pronounced peak and subsequently decreases due to the rapid longitudinal expansion of the Glasma. However, in addition to the energy deposition in the central region, we also observe rather large contributions to T 00 , T 0z and T zz in the vicinity of the two light-cones. While it is intuitively clear that the nonequilibrium plasma produced away from central region should feature sizeable velocities in the longitudinal direction and therefore contribute significantly to T 00 , T 0z and T zz in the lab-frame, the magnitude of contributions is surprisingly large compared to the transverse pressure. Even though we can not clearly rule out that these contributions may arise due to artifacts of the lattice discretization, we have checked explicitly that the observed behavior remains unchanged when we decrease the lattice spacing, as discussed in more detail in Appendix A. Since to the best of our knowledge such behavior has not been reported previously in the context of 3+1D Glasma simulations, clarifying the detailed structure of the fields in the vicinity of the light-cone will require further numerical and analytical investigations in future.
While the results in Fig. 6 were obtained for the collision of rather thick nuclei (Q s R = 1/2), it is also interesting to investigate how the space-time profiles change when varying the longitudinal thickness Q s R of the colliding nuclei. We investigate this behavior in Fig. 7, where we present heat-map plots of the space-time evolution of the transverse pressure (T xx +T yy )/2 for Q s R = 1/2, overlayed with the τ, η coordinate system in the for- transverse pressure is significantly reduced towards the edges of the forward light-cone, indicating a non-trivial space-time rapidity profile around mid-rapidity.

C. Space-time rapidity profiles
So far, we have leveraged our framework to study the space-time picture of the collision in Minkowksi space. Now we will look at the non-trivial rapidity dependence of the observables, which arises naturally by including the longitudinal thickness of the colliding nuclei. Even though the mapping of (t, z) data into (τ, η) coordinates is in principle straightforward, the availability of information on a discrete t, z grid 3 poses additional challenges, as a straightforward interpolation between data points can become problematic in the vicinity of the light cones. Due to these technical difficulties, we will restrict ourselves to an investigation of the rapidity range η ∈ (−1.25, 1.25), and show the corresponding result for different τ as a function of η. We further note that the conversion of Minkowski (t, z) space results to Milne coordinates τ, η can in fact be quite sensitive to the definition of the origin of the coordinate system, and we will always fix the origin t = 0, z = 0 at the space-time point, where the center of mass of the two nuclei coincides.
In Fig. 8, we present the evolution of transverse pressure τ (T xx +T yy )/2 as a function of η for different values of Q s R. Different color codings in Fig. 8 correspond to the results obtained at different proper times τ , and the scaling of the transverse pressure by a factor of τ has been chosen such that -beyond time scales τ ∼ 1/Q s - the quantity τ (T xx + T yy )/2 shown in Fig. 8 becomes independent of proper time τ in the boost-invariant highenergy limit. Starting from the collision of thin nuclei with Q s R = 1/16, we observe the emergence of a boostinvariant plateau for η ∈ (−0.8, 0.8), as already seen in the bottom panel of Fig. 7, where shortly after the collision the contours of constant transverse pressure follows the line of constant proper time τ . When increasing the longitudinal thickness of the colliding nuclei, the transverse pressure of the Glasma created around midrapidity decreases and we see how the profiles are no longer flat around the central rapidity even for the late times Q s t Q s R. Empirically, we find that the resulting rapidity profiles can be reasonably well described by the following functional form, which is indicated by the black dashed lines in Fig. 8, where τ P T (η = 0) is the pressure at mid-rapidity and η R describes the rapidity width. By looking at the extracted values of η R and τ P T (η = 0) in Tab. I, one observes that the rapidity width η R exhibits a strong dependence on the width Q s R of the colliding nuclei, whereas the transverse pressure τ P T (η = 0) of the Glasma at midrapidity only decreases slowly with increasing thickness of the colliding nuclei, as can also be seen directly from Fig. 8. When analyzing the rapidity dependence of the other components of the energy-momentum tensor, it is convenient to switch to the local rest frame (LRF), defined by the condition that u µ LRF is a time-like eigenvector of the energy-momentum tensor T µ ν u ν LRF = LRF u µ LRF . By diagonalizing the average stress-energy tensor T µ ν,Glasma of the Glasma, 4 one gets the energy density and longitudinal pressure P L in this frame as We show our results in Fig. 9, where we present results for the energy density τ LRF and longitudinal pressure 4 Note that, as discussed previously, the energy-momentum tensor of the Glasma T µν Glasma is obtained by subtracting the contributions of the Weiszaecker-Williams fields of the colliding nuclei, prior to the diagonalization procedure P L LRF for different values of Q s R. Starting with the evolution of the longitudinal pressure depicted in the top panel, we find that for collisions of thick nuclei, the longitudinal pressure almost vanishes as the two nuclei have passed through each other, as seen for Q s R = 1/2 at late times. With decreasing thickness, the longitudinal pressure starts out from negative values around mid-rapidity, and subsequently relaxes towards zero, in qualitive agreement with the well established behavior in the high-energy boost-invariant limit [9,30]. Conversely, the rise of the longitudinal pressure P L LRF at larger rapidities signifies clear deviations from boost invariance, and can be attributed to the spurious presence of fields on/near the light-cones in Fig. 6 In the bottom panel of Fig. 9 we present energy density for different longitudinal extent of nuclei, at late times where longitudinal pressure almost vanishes (top panel), such that LRF 2P T and τ LRF const is approximately constant. Similar to Fig. 8, we again notice the emergence of a boost-invariant plateau in the high-energy limit (Q s R = 1/16), whereas for the collision of thick nuclei at lower energies (Q s R = 1/2, 1/4) we see contrasting result which again signifies broken boost-invariance.

IV. 3+1 D COLLISIONS WITH REALISTIC COLOR CHARGE DISTRIBUTIONS
So far we have considered a simplistic model of color charge distributions inside each nucleus, to perform a detailed investigation of the longitudinal dynamics of the Glasma during and shortly after the collision. Evidently, to connect these simulations to realistic heavy-ion collisions, it is necessary to develop a more physical model of the color charge distributions, which reflect both the longitudinal and transverse structure of the colliding nuclei. Similar to the discussion in the boost invariant high energy limit [10,39,40], the basic idea of our construction will be to connect the color charge distributions ρ L/R (x ± , x ⊥ ), to measurements of hadronic structure functions from the deep-inelastic scattering experiments. Below we develop a model of the three-dimensional structure of the color charge distribution based on the smallx transverse momentum distribution (TMDs), and subsequently perform simulations within this framework to study the effect of longitudinal fluctuations of the color charge distributions.

A. Connection to small-x TMDs
Generally speaking, the three-dimensional parton of nucleons and nuclei is encoded in an underlying Wigner distribution [41,42] that contains information on the position and momenta of single partons inside a nucleon or nucleus. By disregarding position or momentum information, the Wigner distribution reduces to a transverse momentum parton distribution (TMD) or respectively to a generalized parton distribution (GPD). Conversely, if both position and momentum information are discarded, one obtains the standard collinear parton distribution function (PDF).
Even though a modeling of the color charge distribution based on the Wigner function would be desirable, little is known about this fundamental object, and we will therefore consider a parametrization of the color charge densities based on TMDs, with the three dimensional spatial structure of nucleons and nuclei imposed by hand according to a Monte Carlo Glauber model. Specifically, we will assume that the position and momentum dependence of the color charge distribution inside a nucleon can be factorized as where T x+y 2 captures the spatial structure of the colliding nucleon, and thus varies on length scales ∼ R p and ∼ R p /γ, where R p is the proton radius, whereas the Fourier transform of Γ (x − y) describes the transverse and longitudinal momentum dependence of color charges inside the nucleus, such that e.g. in the transverse plane Γ (x − y) typically varies on distance scales ∼ 1/Q s . Now in order to constrain the behavior of Γ (x − y), we consider the operator definition of the dipole gluon TMD for a left moving nucleus [43,44] where x 2 is the longitudinal momentum fraction and k ⊥ is the transverse momentum of the gluons. The gauge links U [ξ,ξ ] and U [ξ ,ξ] connecting the points ξ and ξ ensures a gauge invariant definition of the TMD distributions.
Within the Color Glass Condensate effective theory, the matrix element p A |...|p A / p A |p A , is replaced by an average . over the color charge distribution [43,44] Based on eqn. (4,14) the non-abelian field strength tensor F i− (ξ) and gauge links U [ξ,ξ ] , U [ξ ,ξ] , can be calculated as functional of the color charge density ρ, such that the is in fact entirely determined by specifying the n−point correlation functions of the color charge density. Evidently, the general relation between the correlation functions of ρ and the dipole TMD is rather complicated, and we will thus simplify the problem by considering Gaussian correlations of color charges in the dilute limit, where the expression in Eq. (40), can be expanded to lowest non-trivial order in ρ's. By evaluating the twopoint correlation functions according to Eq. (38) using average and difference coordinates, we then obtain where p − = s N N /2 is the large light-component of momentum of the nucleon in the lab frame and S ⊥ = dξ + d 2 ξ ⊥ T (ξ ⊥ , ξ + ) is the transverse area of the nucleon. Based on Eq. (41), whereΓ denotes the Fourier transform of Γ(x−y), one then concludes that the k ⊥ dependence of the x 2 G (2) (x 2 , k ⊥ ) determines the transverse structure of the correlation function Γ(x−y), whereas the longitudinal structure of the correlation function Γ(x−y) is related to the x 2 dependence of the TMD. Now in order to employ Eq. (41) to determine the correlation function, we still need a parametrization of the gluon TMD x 2 G (2) (x 2 , k ⊥ ) as input to the calculation. We will employ the GBW Model [45] -a simple phenomenological model that has been fit to small-x deepinelastic scattering data -where with the saturation scale Q s (x) parameterized as [45] Q with λ = 0.144 and the value of Q 0 being set to 0.5 GeV.
Since the color charges are assumed to be x − independent, the light-cone component k + vanishes identically for each source, such that in terms of the spatial momenta k − = − √ 2k z and the momentum fraction is given by x 2 = k − p − . Based on Eq. (42), the correlation functioñ Γ(k ⊥ , k z ) of the initial color charges is then obtained as which can be used to sample individual configurations of the color charge density as discussed below.

B. Sampling of realistic color charge distributions
When describing the color charge distributions of atomic nuclei, we follow the Monte Carlo Glauber Model and sample the position x i of the i = 1, · · · , A individual nucleons according to a Wood-Saxon distribution. Each individual nucleon is assigned a three-dimensional thickness profile such that the overall thickness of the nucleus is given by which according to Eq. (41) is normalized such that d 3 x T (x) = 2πR 2 p A, where R p = 2GeV −1 is the proton radius. Since the spatial distribution T typically varies on distance scales ∼ R p 1/Q s , the color charge distribution inside the nuclei is then sampled according to ρ a L/R (t, x, y, z) = ga x a y a z T (x, y, z) (47) where ζ(x ⊥ , z) = 1/a 3/2 χ RN G (x ⊥ , z) are Gaussian random numbers andζ(k ⊥ , k z ) denotes their Fourier transform. Based on the initial conditions for the charge density profiles in Eq. (47), we then proceed as described in Sec. II to set up the initial conditions and simulate the dynamics of the collision. We note that due to the presence of longitudinal fluctuations, a finer discretization in the longitudinal (z) direction is required in this case, and unless stated otherwise, we will employ Q 0 a ⊥ = 0.33 and R p /a z = 256 in our numerical studies.

C. Numerical results for realistic charges
We now proceed to the study of the collision dynamics for realistic charge profiles, and consider head on (b = 0) Au-Au collisions for center of mass energies of 130 and 200 GeV. Since the basic features of the reaction dynamics remain essentially the same as for the simplistic charge profiles discussed in Sec. II, we will focus our investigation on the violations of boost invariance and study the longitudinal fluctuations which emerge naturally within our framework.
We illustrate the full 3 + 1 D structure of the events in Fig. 10, which shows the different phases of the collision for a head-on Au-Au collision at √ s N N = 200 GeV. In the left panel, we show the energy density T 00 W W associated with the color fields of the incoming nuclei, which are well separated at t = −0.37 fm/c before the collision. Since at late times the rapid longitudinal expansion of the system stretches out the longitudinal profiles, some of the features of the Glasma are not clearly visible in the above figure and it is advantageous to visualize the corresponding structure in Milne coordinates. We illustrate this in Fig. 11, where we show the three dimensional profile for a head-on Au-Au collision at 130 GeV in x, y and η coordinates at τ 0.4 fm/c, which characterises the rapidity fluctuations. We note that due to limited availability of points near the light cone, we restrict ourselves to the central rapidity range η ∈ [−0. 8  ward in various models of the longitudinal structure of the initial state [46,47]. By careful inspection of individual flux tubes in Fig. 11 one also observes longitudinal fluctuations, albeit the amplitude of the longitudinal variations is significantly smaller compared to the variations in the transverse plane.
Now in order to further analyze the longitudinal fluctuations, we consider the rapidity profiles of the transverse pressure P T (τ, η) ⊥ averaged over the transverse plane. In Fig. 12 we show the ratio of P T (τ, η) ⊥ relative to its value P T (τ, η = 0) ⊥ at mid-rapidity at τ 0.75 fm/c for head-on Au+Au collisions at two different energies √ s N N = 130, 200 GeV in the top and bottom panels.
Different curves in each panel correspond to the results for five individual events (labeled as Seed 0-4), along with the symmetrized average over all events. Generally, the fluctuations in the accessible rapidity window −0.6 < η < 0.6 are relatively small 1%, and appear to decrease with increasing center of mass energy, as the longitudinal profile gets stretched out over a larger rapidity range. It is also interesting to observe that the dominant fluctuation in individual events around mid-rapidity appears to be a forward/backward asymmetry, and it would be interesting to extend this analysis to larger rapidities. However, this would require simulations on significantly finer and larger lattices, and is therefore beyond the scope of the present work.

V. CONCLUSION & OUTLOOK
We developed a framework to perform 3 + 1D simulations of initial energy deposition in heavy-ion collision based on the effective theory of CGC, which takes the finite longitudinal extent of the colliding nuclei into account. Based on a simple model of the color charge distribution, we investigated the detailed dynamics, during and shortly after the collision. While in low energy collisions, where the longitudinal extent of the incoming nuclei Q s R is non-negligible, significant violations of boost invariance can be observed, the results smoothly approach the boost invariant limit [30,33] at high energies where the longitudinal thickness Q s R → 0 becomes sufficiently small.
Subsequently, we developed a more elaborate model of the three dimensional color charge distribution in a large nucleus, where the large scale structure of the nucleus is determined by the longitudinal and transverse positions of nucleons, while small scale fluctuations in the longitudinal and transverse directions are determined by the x and k ⊥ dependence of transverse momentum dependent parton distributions. Based on this model, we obtained first results regarding the three dimensional structure and its fluctuations at two different center of mass energies, which show encouraging trends e.g. the longitudinal rapidity profiles and fluctuations appear to become stretched with increasing center of mass energy, which was not necessarily the case in a previous attempt to generalize the IP-Glasma initial state to 3 + 1 dimensions [25].
Due to the significant computational cost of performing 3 + 1D classical Yang-Mills simulations of the space-time dynamics our numerical results have so far been limited to the central rapidity window for a few head-on collisions, and it would certainly be interesting to extend the analysis to larger rapidities and higher center of mass energies and perform a more systematic study of the various effects as a function of the centrality of the events. Evidently, to make contact with experimental observations, such phenomenological studies should be performed within the physical SU (3) gauge group of QCD, where the formalism developed within this paper can be applied in exactly the same way, albeit further increasing the computational cost of the simulations. Beyond the improvement of numerical simulations (see also [28]), it would also be important to develop further analytical insights into the 3 + 1D space-time evolution of the Glasma, which could e.g. be obtained by analyzing the perturbative dilute limit along the lines of [23].
Below we provide additional results for simulations where we vary the various discretization parameters, to illustrate the results presented in the main part of this work do not suffer from significant discretization artifacts.
We present a compact summary of the results in Fig. 13, where for a fixed value of Q s t 4, subtracted T 00 and T 0z components of the energy momentum tensor are shown for the thick nuclei (Q s R = 0.5). In the first two panels: T 00 obtained with the lattice discretization used throughout the manuscript (R/a z = 16; dt/a z = 0.08) is compared against the values obtained with finer longitudinal lattice spacing (R/a z = 64) and finer time step (dt/a z = 0.04). Lattice dimension is taken to be 128 2 ×1024 and 128 2 ×2048 for R/a z = 16 and R/a z = 64 respectively. We observe that approaching the continuum limit doesn't shrink the spurious fields in the proximity of light cone. In the other two panels, similar result is shown for T 0z , which again shows that the effect of these contributions do not change with finer lattice spacing.