Energy-momentum tensor of the dilute (3+1)D Glasma

,


I. INTRODUCTION
Relativistic heavy-ion collisions (HIC) are studied experimentally at the Large Hadron Collider (LHC) and Relativistic Heavy Ion Collider (RHIC) as a method to probe quantum chromodynamics (QCD) at extreme energies [1][2][3], where quarks and gluons are liberated from color confinement and form the Quark Gluon Plasma (QGP).The intricate spacetime evolution of the system, coupled with changing degrees of freedom, presents a difficult challenge to first principles calculations and phenomenological model building.The current state-ofthe-art understanding of the creation and evolution of the medium relies on sophisticated multi-stage modeling, typically split into a pre-equilibrium stage, a hydrodynamic stage, and a hadron gas stage [4].Bayesian analysis in the context of multi-state examination reveals correlations between initial conditions and medium properties [5][6][7], but it is not clear to what extent the properties of the initial state of a heavy-ion collision can be deduced from experimentally measured final state observables.
Despite these efforts, our current knowledge of the initial conditions in three spatial dimensions for high-energy heavy-ion collisions remains notably limited.The Color Glass Condensate (CGC) [8,9] emerges as a promising theoretical framework for providing first principles insight into the initial state, wherein the initial interaction of nuclei is characterized in terms of large and small Bjorken-x partons, leading to a dense state of gluonic matter known as the Glasma [10].Due to high gluon occupancy, the leading-order dynamics are governed by the Classical Yang-Mills (CYM) equations.The Glasma is thus characterized by strong classical color fields.
In contrast to phenomenological models, such as e.g.MC Glauber-type models [11], which aim to directly model the energy density of the collision medium, the CGC/Glasma approach provides a description for nuclei before the collision, as well as the collision itself, and the subsequent pre-equilibrium evolution of the Glasma.Since the properties of the Glasma are fully determined by the CYM dynamics and the cold nuclear matter properties of atomic nuclei before the collisions, the modeling aspect of the pre-equilibrium stage is relegated to the distribution of color charge of high-energy nuclei.The IP-Glasma [12][13][14] is one such nuclear model with notable phenomenological success.
Most CGC-based initial state models neglect the longitudinal dynamics of the medium due to the assumption of boost invariance at high energies, effectively rendering the Glasma a (2+1)D system.While boost invariance provides a good approximation for central collisions at mid-rapidity, recent experimental measurements carried out at the LHC, such as the observation of flow decorrelations [15,16] present a challenge to the assumption of boost invariance in the initial state.Consequently, it becomes essential to delve into the longitudinal dynamics of high-energy collisions.Over the years, various implementations of (3+1)D initial state models have been developed, employing either static [17][18][19] or dynamical sources [20][21][22][23][24][25][26] in the CGC framework, or phenomenologically by considering constituent quarks, color flux tubes or strings of varying lengths [27][28][29][30][31].However, these approaches are either phenomenological in nature or severely constrained by the computational resources required.
In this paper, we present the extension of the first analytical calculation of energy deposition in high-energy heavy-ion collisions obtained by solving the (3+1)D CYM equations within the weak field approximation [32].The results exhibit excellent agreement with fully nonperturbative lattice simulations that share a similar origin of rapidity dependence [32,33].We now employ the derived analytic expressions for the perturbative gauge fields on the future light cone to formulate the energymomentum tensor of the Glasma.A key focus of this study lies in presenting a concise expression for the field strength tensor of the Glasma, wherein we rigorously reduce the number of integrals from six to three, optimizing efficiency for Monte Carlo integration.Subsequently, we employ these expressions in conjunction with a threedimensional generalization of the MV model at two different collision energies corresponding to RHIC and LHC ranges to examine the rapidity profiles of the energymomentum tensor, local rest frame (LRF) energy density and flow.
In Sec.II we discuss the weak field approximation in the CGC and show how to obtain the Glasma field strength tensor and Glasma energy-momentum tensor up to leading order in the sources.The resulting integrals need to be solved numerically, for which we explain our Monte Carlo (MC) integration method.Next, we introduce a shifted Milne frame for evaluating the energymomentum tensor of the dilute Glasma and derived observables in the future light cone.In Sec.III we explain the details of our three-dimensional nuclear model, which includes a correlation length parameter that allows us to set the scale of longitudinal fluctuations in our nuclei.In Sec.IV we discuss our numerical results for several observables of the dilute Glasma.We summarize our findings and give a brief outlook in Sec.V.

II. WEAK FIELD APPROXIMATION
The CGC is an effective field theory for high energy QCD applicable to relativistic heavy-ion collisions [8,9].The hard degrees of freedom within the two nuclei are modeled as classical currents of color charge moving at the speed of light.In contrast, the soft gluon fields are sourced by the nuclei and determined by the Classical Yang-Mills equations.While the classical approximation is warranted due to the high gluon occupation numbers and weak coupling at high energies, it should also be noted that it agrees with the quantum treatment to leading order in perturbation theory.
We use the CGC to describe the initial states and collision of two relativistic nuclei.The geometry of our setup is shown in Fig. 1.First, consider a single nucleus A moving in negative z-direction through the (pre-collision) Region I at the speed of light.The corresponding color current in light cone coordinates Region II

Region I
FIG. 1. Spacetime diagram of a relativistic heavy-ion collision showing time t and beam axis z.The tracks of nuclei A (B) are marked in orange (blue).The diamond-shaped region where they overlap is the interaction region, and together with its causal future, it is called Region II.The complement region, which represents the system before the collision, is called Region I.
The boldface symbol x refers to the transverse coordinates x = (x 1 , x 2 ) = (x, y) and ρ A (x + , x) is the color charge density, which takes values in the Lie algebra su(N c ). Due to Lorentz contraction, the charge density is localized in a thin sheet around x + = 0. We allow for a finite longitudinal thickness of the nucleus as opposed to the ultrarelativistic (boost-invariant) case, where the nucleus is considered to be infinitesimally thin.Within our treatment, the nucleus is still considered a static source, i.e., the current does not depend on x − and propagates at the fixed speed of light.Corrections to this behavior provide another source of sub-eikonal corrections [34][35][36][37][38].
In Eq. ( 1) and in the following, we always assume covariant gauge, which for the gauge field with the field strength tensor Using the boundary conditions where i = x, y labels the two transverse components, the YM equations are solved by with ∇ ⊥ acting in the transverse plane.This gauge field and the current satisfy the continuity equation Analogously, introducing a nucleus B moving in the opposite direction with current yields A. Dilute limit of the Glasma field We now review our treatment of the collision of nuclei A and B in terms of the dilute approximation (see [32,33] for details).In the interaction region, where the two nuclei overlap, and its causal future, denoted together in Fig. 1 as Region II, the solutions in Eqs. ( 5) and ( 8) are no longer valid and neither is a superposition.Instead, the solution to the full collision problem has to be obtained from the Classical Yang-Mills equations Here, non-calligraphic symbols refer to the interacting solutions of the field equations.The corresponding gauge field and current (11) are defined in terms of the non-interacting solutions before the collision plus additional correction terms a µ and j µ , which are zero prior to the collision and depend on x = (x + , x − , x).Therefore, in the asymptotic past To make further progress, we consider an expansion of Eqs.(10) and (11) in powers of ρ A and ρ B .Since we work in the covariant gauge, ∂ µ a µ = 0. Note that the single nuclei gauge fields A µ A/B are linear functionals of ρ A/B as expressed in Eqs. ( 5) and ( 8), i.e. they are of order ρ A/B .Similarly, the currents J µ A/B are also of order ρ A/B .If ρ B = 0, then A µ = A µ A and J µ = J µ A solve Eq. ( 9) to all orders in ρ A .The analog is true for A µ = A µ B and J µ = J µ B if ρ A = 0. Therefore, a µ and j µ capture corrections of order ρ n A ρ m B with n, m ≥ 1 in the full collision problem.Expanding Eq. ( 9) in orders of ρ A and ρ B the contributions of order ρ A ρ B are where we take a ν to be of order ρ A ρ B with all higher order terms removed.Dropping terms of higher order in ρ A/B , we implicitly assume weak sources.Therefore, we refer to this approximation as the dilute limit or weak field approximation.Considering only the contributions of order ρ A ρ B to D µ J µ = 0 yields As the nuclei are in recoilless, lightlike motion, their trajectory does not change.Therefore, j µ only represents a rotation in color space of the nucleus' currents and is limited to j − and j + components, which are localized in the same region as ρ A and ρ B , respectively.Equation ( 15) then splits into two contributions Given the initial conditions in Eq. ( 13) they are solved by where ϕ A/B correspond to the gauge fields of nuclei A/B in Eqs. ( 5) and (8).These expressions can be used in Eq. ( 14) to solve for a µ (x), which we identify as the Glasma field.We recite the solutions from [32], Here, f abc are the structure constants of the su(N c ) Lie algebra and t c are its generators.The symbol |p| denotes the modulus of the transverse vector p, and p•x = p i x i = p i x i .Quantities with a tilde are understood as Fourier transforms in the transverse plane, e.g.
where p = , and J m are the Bessel functions of the first kind.We note that these solutions are only correct inside the future light cone, i.e. outside the tracks of the nuclei, as there are additional corrections inside the tracks.However, these are co-moving with the color fields of the single nuclei, such that these corrections do not contribute to the color field of the Glasma in the future light cone and thus can be ignored.

B. Field strength tensor
The next step is to compute the leading order contribution to the field strength tensor f µν from the leading order gauge fields a µ .In general, the non-Abelian field strength tensor (cf.Eq. ( 3)) contains an Abelian part ∂ µ a ν − ∂ ν a µ and a commutator term [a µ , a ν ].In the dilute limit, the perturbative Glasma fields a µ are of order ρ A ρ B , and the commutator term only contributes to higher-order corrections.Thus, to leading order, the Glasma field strength tensor can be defined as1 After a lengthy derivation, which can be found in Appendix A, we arrive at a remarkably simple result.The components of the leading order field strength tensor are given by where the rapidity-like integration variable η ′ ∈ (−∞, ∞), and the v-integral spans the entire transverse plane, v = d 2 v.The symbols V and V ij are defined as where β i are the gradients of the color potentials ϕ, In covariant gauge, these gradients may be identified with the only non-trivial components of the field strength tensors of the single nuclei, The vector components w j are given by These concise expressions for the Glasma field strength tensor are our main analytical result and we stress that the boost-invariant limit β i,a A/B (x ± , x) = δ(x ± )α i,a A/B (x) of our expressions is consistent with previous boostinvariant results [39,40].
Our solutions for the field strength tensor have an intuitive geometric interpretation illustrated in Fig. 2. Given the spacetime point x, at which we evaluate f µν , the integration is restricted to the causal past of x along lightlike paths emitted from the collision region.To see this, we note that the modulus of v in Eqs. ( 25) -( 28) can be interpreted as a time coordinate τ ′ = |v| analogous to Milne proper time τ .Together with the rapidity-like coordinate η ′ , we find that the pair (τ ′ , η ′ ) can be viewed as Milne coordinates for a slice with fixed |v| of the past light cone attached to the spacetime point x.We emphasize that the displacement four-vector measures the spacetime distance between the source of emission and the point x.Since this vector is lightlike, i.e. v µ v µ = 0, it is apparent that the three-dimensional integrals in Eqs. ( 25) -( 28) only sum over lightlike paths ending at x. Furthermore, we note that the integrands are non-zero only in the four-dimensional spacetime volume where the two colliding nuclei overlap.Intuitively, this region corresponds to points in spacetime from which gluons are emitted due to interactions among the colliding color fields.These gluons then propagate along lightlike paths into the future light cone until they arrive at x, forming the Glasma field strength at that particular point.
Compared to the gauge fields in Eqs. ( 20)-( 22), the components of the field strength tensor f µν are much easier to evaluate numerically.Firstly, the gauge fields are six-dimensional integrals, whereas the field strength tensor only involves three-dimensional integrals.This simplification is possible by exploiting the closure relation of the Bessel functions J m , see Eq. (A19) in the Appendix for details.Secondly, having eliminated the Bessel functions, the integrands of f µν are now free of oscillating terms.This way, the numerical evaluation and convergence of Eqs. ( 25) -( 28) is greatly simplified.

C. Monte Carlo integration
Despite the relatively simple structure of the expressions, the integrals for the components of the field strength tensor have to be solved numerically in practical computations of the (3+1)D Glasma.In this work, we focus on a Monte Carlo integration approach.We consider the prototypical integral where we assume that the function h(u + , u − , u) has compact support along the two light cone directions, u ± min < u ± < u ± max with longitudinal widths 2. Geometric illustration of the integration paths for the field strength tensor f µν .The integration domain is the past light cone of the spacetime point (x + , x − , x), parametrized by the Milne coordinates (τ ′ = |v|, η ′ , v).For a slice at fixed transverse displacement τ ′ = const, we highlight the conic section inside the collision region that contributes to the integral (orange solid line).
a slight error because the color fields ϕ A/B (x ± , x), which enter the function h, generally only exhibit exponential decay in both longitudinal and transverse directions.
The constants u ± max/min can nonetheless be chosen large enough such that our numerical results are barely affected by such a cut-off.Similarly, we restrict the transverse integration domain over u to a rectangular region u i min < u i < u i max with u i max > u i min , where the color fields overlap.Furthermore, we evaluate the integral I only in the region where Eqs. ( 25)-( 28) are valid, that is for x ± > u ± max .With these strong assumptions, the integrals are guaranteed to converge. 2 Our integration strategy is based on rewriting the integral in Eq. ( 36) using a coordinate transformation and making use of the restricted integration domain.Since Monte Carlo integration relies on random sampling, suboptimal sampling strategies can lead to unnecessary eval- 2 The conditions for convergence can be weakened if we assume that the color potentials ϕ A/B (u ± , u) are not compact along the light cone directions, but merely decay to zero for |u ± | → ∞.
In this case, the overlap of the two nuclei shown in Fig. 2 is unbounded and the evaluation point x is always inside the overlap region.The convergence of f +− and f ij is guaranteed if the color fields do not exhibit any singularities and decay faster than On the other hand, the components f ±i in Eqs. ( 26) and ( 27) impose a stronger decay due to the presence of the exponential factors e ±η ′ .Rewriting these integrals using light cone coordinates v ± yields e ±η ′ ∝ v ± /v ∓ (see Appendix A).To ensure convergence for v ∓ → ∞, the color potentials ϕ A/B have to asymptotically decrease as |u ± | −s with s > 3/2.Since nuclear models impose at least exponential decay, these conditions are always fulfilled.Nonetheless, we stress that the derivation of Eqs. ( 25)-( 28) in principle assumes compact support.Evaluating the integrals inside the overlap yields an incomplete result for the Glasma field strength.
uations of the integrand in regions where it is zero.We can increase efficiency by determining strict bounds on the coordinate ranges.First, we use polar coordinates (v, θ) for the transverse integration over v, where v = v e(θ) with the unit vector in the v-direction e(θ) = (cos θ, sin θ).Secondly, the aforementioned restrictions imposed on the integration domain allow us to put bounds on the integrals over η ′ , v, and θ.Considering only the restrictions on u + , we find that η ′ is restricted to (η + min (v), η + max (v)) with for a given value of v > 0. Analogously, we find a restriction from Combining both inequalities for η ′ , the integration domain must be restricted to η ′ ∈ (η min , η max ) with η min = max(η + min , η − min ) and η max = min(η + max , η − max ).We also deduce bounds on the integration variable v. From Eq. (37) it follows that Saturating the restrictions on u ± , we find that v ∈ (v min , v max ) with The angle θ is also restricted.Given valid choices of η ′ and v, the angular integration corresponds to a circle of radius v with center x in the transverse plane, intersecting the rectangular integration domain only at certain points.Only those segments of the circle that lie inside the restricted integration domain contribute to the integral, as illustrated in Fig. 3, where the intersection of the circle with the integration domain is a single arc.Determining these arcs then allows for more efficient sampling of θ.
We perform random sampling as follows: we first choose a radius v with probability density p(v) ∝ v inside the bounds (v min , v max ).Given v, we evaluate the bounds on η ′ and sample uniformly from the interval given by η min (v) and η max (v).Finally, we determine where the circle defined by v and x intersects the rectangular domain and obtain the segments from which we sample θ u v x b FIG. 3. Illustration of two partially overlapping nuclei (orange blobs) separated by the impact parameter b in the transverse plane spanned by u.To evaluate the field strength tensor at x, we restrict the angular integration with radius v to the arc (thick blue) within the rectangular region where the two nuclei overlap.
uniformly.Drawing M samples for all three coordinates, we approximate the integral by where ) and C is a Jacobian factor, which is proportional to the integration volume.

D. Energy-momentum tensor
From the Glasma field strength tensor given in Eqs. ( 25)- (28), it is straightforward to obtain the Glasma energy-momentum tensor in the future light cone.In addition to the energymomentum tensor, we compute the local rest frame (LRF) energy density ϵ LRF and fluid velocity u µ by solving the Landau condition with u µ the only timelike eigenvector of T µ ν and ϵ LRF its eigenvalue.We parametrize the future light cone in Milne coordinates For the origin S (orange), the hyperbola enters the trajectories of the nuclei (gray).This is never the case for the Milne frame centered at S (blue).
where τ and η s are proper time and spacetime rapidity, respectively.Contrary to the boost-invariant setup, it is a priori not clear where the origin of the Milne frame should be placed.As depicted in Fig. 4, putting the origin at the center of the interaction region (S) will lead to some of the τ = const hyperbolas cutting into the tracks of the nuclei.In order to avoid this, we shift the origin of our Milne frame forward in time by some offset δt that is slightly larger than the nuclear radius in the longitudinal direction.None of the τ = const hyperbolas with respect to the new origin S cross the trajectories of the nuclei.We evaluate all observables in terms of this shifted Milne frame, but we omit the bars in τ and ηs from now on to reduce clutter.

III. NUCLEAR MODEL
We study a three-dimensional generalization of the McLerran-Venugopalan (MV) model [41,42].We assume a Gaussian probability functional for the color charges of our nuclei, propagating either in x + or in x − direction, which is completely fixed by the 1-and 2-point functions.The 1-point function in our model vanishes, The 2-point function is given by with the MV scale µ, which, in the non-perturbative case, is related to the saturation momentum Q s ≈ g 2 µ.Here, is a boosted Woods-Saxon (WS) profile with WS radius R and skin depth parameter d.The Lorentz contraction factor γ beam is determined by the beam velocity and the normalization constant c ensures The function T determines the overall shape of the nuclei.
On the other hand, fluctuations within the nuclei are controlled by the function where ξ is the longitudinal correlation length with 0 ≤ ξ ≤ 2R l and R l = R/( √ 2γ beam ) is the boosted WS radius.A similar model was used for single nucleons in [25].The R l dependent factor in Eq. ( 53) guarantees consistency with two limiting cases w.r.t. the longitudinal correlation length ξ.
In the limit ξ → 0, Eq. 53 becomes and thus the correlator in Eq. ( 50) can be written as In this limit, our model is similar to the traditional MV model with longitudinal extent (see e.g.[43] or [44]).Introducing a transverse charge density we can integrate out the longitudinal dependence in Eq. ( 55) and obtain The normalization in Eq. ( 52) now ensures that we recover the two-dimensional MV model [41,42] at the center of our nuclei, x = 0. We refer to this limit as the MV model limit.
On the other hand, we may take ξ → 2R l and obtain lim In this case, the charge densities can be written as where the two-dimensional Gaussian random field χ ⊥ is given by The charge density does not fluctuate along the longitudinal coordinate and is merely modulated by the enveloping profile T (x ± , x).This is the coherent limit of our model.Initial conditions with these particular longitudinal correlations have been previously studied using (3+1)D lattice simulations [20,21], albeit without the non-trivial transverse structure imposed by T (x ± , x).
To sample charge distributions from our nuclear model for arbitrary ξ, we first generate random Gaussian noise χ a (x ± , x), which fulfills We then move to Fourier space with respect to x ± and introduce a new field After transforming back to position space, we obtain the single nucleus color charge densities as Note that the individual color components ρ a are statistically independent.Once we have sampled charge densities ρ A/B for both nuclei, we obtain the corresponding single nucleus color fields by solving the transverse Poisson equations ( 5) and ( 8) in momentum space.We obtain where we have introduced an infrared (IR) cutoff m and an ultraviolet (UV) cutoff Λ UV .In practice, we sample the nuclear model on a discrete lattice.The sampling procedure can be carried out numerically using the Fast Fourier Transform (FFT) algorithm, see e.g.[24,25].
In Fig. 5 we show one component of the color field ϕ(x ± , x) of a single nucleus sampled from our nuclear model.We cut along the center of one transverse direction and compare different values of the IR cutoff m and longitudinal correlation length ξ, showing positive values in orange and negative values in blue.While the magnitude of m governs the size of transverse fluctuations, ξ determines the longitudinal structure of our nuclei.The right panels show the coherent limit ξ = 2R l .

IV. NUMERICAL RESULTS FOR THE (3+1)D STRUCTURE OF THE GLASMA
We now discuss the (3+1)D spacetime structure of the Glasma emerging from our model.We compute the integrals in Eqs. ( 25)-( 28) using the Monte Carlo approach outlined in Sec.II C. From the Monte Carlo integration, we obtain estimates for the Glasma field strength tensor f µν and energy-momentum tensor T µν given in Eq. (45).The raw data files for T µν are published in [45,46].We solve the Landau condition Eq. ( 46) to find the local rest frame energy density ϵ LRF and flow velocity u µ .We obtain estimates for the magnitude of the Monte Carlo errors by performing standard jackknife analysis and we choose the number of Monte Carlo samples large enough to make them negligible.In practice, we found that M = 10 5 samples are sufficient.
Table I shows the physical parameters of our calculations and the values they are allowed to take.We consider two different collision energies, comparable to Pb+Pb collisions at LHC and Au+Au collisions at RHIC.The choice of collision energy enters our computations through the beam Lorentz contraction factor γ beam .Furthermore, we adapt the parameters of our Woods-Saxon model to the different nuclei in the same way as [13].The values for the correlation length ξ correspond to longitudinal fluctuations at the scale of the nucleon size, ξ = 0.1 R l , the coherent limit, ξ = 2 R l , and an intermediate value, ξ = 0.5 R l .We use N ev = 10 independent collision events to compute event averages of our observables.As shown in Fig. 4, we use a collision origin shifted by δt = (R + d)/γ beam to avoid contributions from the background fields at larger rapidities.We note that the dilute limit differs from the nonperturbative Glasma in one key aspect: due to the perturbative expansion, we find that the MV scale parameter µ, or more generally the energy scale g 2 µ, appears in our analytic results only as a prefactor in f µν or T µν .The scale g 2 µ does not appear as a transverse momentum scale in the dilute limit.Instead, the transverse structure of the nuclei and the resulting Glasma is determined solely by the infrared regulator m.We thus interpret m as the analog of the saturation momentum Q s in the nonperturbative Glasma, and g 2 µ as an energy scale that requires calibration using either experimental results or a fit to a non-perturbative lattice simulation.For example, using the dilute Glasma as an initial stage for a full simulation of a heavy-ion collision, the parameter g 2 µ may be fixed such that the charged particle multiplicity at midrapidity is correctly reproduced at either RHIC or LHC.In the present work, we simply set g 2 µ = 1 GeV, keeping in mind that phenomenological applications require a properly chosen value.
Figure 6 shows a perspective plot of the threedimensional local rest frame energy density ϵ LRF for a single Au+Au collision event with an impact parameter b = R at RHIC energy at τ = 0.4 fm/c.We show transverse and longitudinal slices of the energy distribution in Fig. 7.In the left panel, we observe the typical almond shape induced by the non-zero impact parameter.The right panel depicts elongated, approximately boostinvariant structures with varying transverse extents reminiscent of Glasma color flux tubes [10,47].These structures are analogous to the flux tube structure employed in several initial state models for describing longitudi- nal correlations in the initial stage of heavy-ion collisions [48,49].

A. Longitudinal structure and flow
The three-dimensional Glasma generated in collisions of longitudinally extended nuclei has a particular longitudinal structure and flow properties that are qualitatively different from the boost-invariant scenario.For practical reasons, we compute the energy-momentum tensor in a (shifted) Milne frame, parameterized by proper time τ and spacetime rapidity η s , as discussed in Sec.II D. We find that the longitudinal flow u η of the collision medium, i.e. the rapidity component of the four-velocity u µ , deviates from the idealized Bjorken case, particularly at large rapidities.This is in contrast to the boost-invariant Glasma, where u η exhibits only local fluctuations [51] and thus T τ τ can be considered as a reasonable measure for the energy density of the Glasma.The Milne frame centered at the collision point is naturally adapted to the symmetries of the boost-invariant Glasma, which (assuming negligible dynamics in the transverse direction) expands with u τ = 1 and u η = 0, known as free streaming.As argued in Sec.II D, there is no unique Milne frame for a collision of longitudinally extended nuclei.The shifted origin of the Milne frame compared to the collision center, along with the extended collision region leads to a violation of the free streaming behavior.The three-dimensional Glasma picks up a considerable longitudinal velocity u η in regions of large rapidity.To avoid ambiguity due to the choice of frame, it is therefore sensible to study frame-independent quantities such as the LRF energy density ϵ LRF , or the transverse pressure T xx + T yy , which is unaffected by longitudinal boosts.
In Fig. 8 we investigate the rapidity dependence of various notions of energy density in the dilute Glasma.We see that the LRF energy density ϵ LRF is in almost perfect agreement with the sum of transverse pressures T xx + T yy .However, note that the T τ τ profiles differ fundamentally from the ϵ LRF profiles at large η s .In addition to the aforementioned energy densities and pressures, we also show γϵ LRF , which is the LRF energy density, weighted by the local flow component γ(τ, η s , x) ≡ u τ (τ, η s , x).This quantity is interesting because γϵ LRF is used to determine the transverse geometric structure of the collision medium in terms of eccentricity coefficients (see Sec. IV D).We find that its rapidity profile is typically wider, but similarly well-behaved at large rapidities as ϵ LRF and transverse pressure.
Figure 8 reveals a considerable difference between T τ τ , which is the energy density reported by an observer moving with u τ = 1 and u η = 0, and ϵ LRF at large rapidities.Intuitively, this means that the flow of the threedimensional Glasma differs from longitudinal Bjorken flow, which makes the Milne frame an ill-adapted coordinate system.Consequently, the Milne frame energymomentum tensor T µν picks up sizeable off-diagonal components at large rapidities.Upon diagonalization and consequential transformation to the LRF, ϵ LRF becomes much smaller than T τ τ and we find significant contributions to the local u η .In addition to longitudinal flow, the mismatch T τ τ ≫ ϵ LRF also arises from longitudinal pressure, as shown in Fig. 9. Particularly for small ξ and m, we find that the longitudinal pressure τ 2 T ηη dominates over the transverse pressure at large rapidi- 8. Different formulations of energy density integrated over the transverse plane as a function of spacetime rapidity ηs at fixed proper time τ = 0.4 fm/c.The results are normalized at ηs = 0 and are averages over 10 central collision events at RHIC energy.The local rest frame energy density ϵLRF (emerald dots) agrees with the transverse pressure T xx + T yy (blue triangles).For small ξ and m (upper left) the T τ τ component increases for large ηs.
ties.This is consistent with the fact that the energymomentum tensor is traceless.Similar to T τ τ , it is apparent that longitudinal pressure is strongly affected by the choice of coordinate system.
Comparison of the components of the energymomentum tensor T µν and the local rest frame energy density ϵLRF integrated over the transverse plane as a function of spacetime rapidity ηs at fixed proper time τ = 0.4 fm/c.The results are normalized to T τ τ (ηs = 0) and are averages over 10 central collision events at RHIC energy.Tracelessness of T µν is manifest with the longitudinal pressure τ 2 T ηη (red triangles) rising and compensating for the large T τ τ (yellow triangles) at extreme ηs.For small m (top row) τ 2 T ηη is negative at central ηs.
We study longitudinal flow by computing an ϵ LRFweighted average of u η in the transverse plane, i.e.
In Fig. 10 we show this weighted average over η s for constant τ = 0.4 fm/c.Regions of positive η s correspond to negative u η and vice versa.This indicates a longitudinal expansion that is slower than Bjorken flow, to which the Milne frame is naturally adapted, and leads to the strong deviations of T τ τ from the rest frame energy density ϵ LRF .We find that the shape of the longitudinal flow can be explained in a simple model, depicted by the solid curves, where we assume that the resulting four-velocity at a given point in the forward light cone can be obtained by a superposition of Bjorken flows starting from each point in the collision region taking into account our forward-shifted coordinate system.This model is derived in Appendix B. Since the data fit the analytic results extraordinarily well, we conclude that the longitudinal flow can mostly be attributed to the extended collision region and the shifted origin of our coordinate system.This reinforces the role of ϵ LRF and T xx + T yy as the more fundamental and physical notions of energy density in the three-dimensional Glasma, as compared to the Milne frame energy density T τ τ .

B. Time and energy dependence
In addition to the rapidity dependence of the energy density, we also study its time and collision energy dependence.In Fig. 11 we show the rapidity profiles of the transverse energy per unit spacetime rapidity dE ⊥ /dη s for different values of proper time τ .These data are obtained by integrating the sum of transverse pressures T xx +T yy over the transverse plane and multiplying with τ to correct for the expected leading order time dependence: We see that the transverse energy per unit rapidity stabilizes on a time scale of ∼ 1.0 fm/c.This is a wellknown result for the boost-invariant case [52,53], where the energy density exhibits a characteristic 1/τ behavior associated with free streaming for τ ≳ 0.2 fm/c.We find that this holds for the three-dimensional Glasma as The rapidity profiles stabilize at late times τ ≳ 0.6 fm/c.For collisions at LHC energy, we observe a plateau around midrapidity, which is absent at RHIC energy.
well, even at extreme rapidities.For LHC energy we observe a stable plateau of roughly ±3.5 units in rapidity around η s = 0 which can be identified with the boostinvariant regime.This plateau is absent at RHIC energy, where the profiles appear to be approximately Gaussian.
A similar plateau for LHC energy appears in the curves of the longitudinal flow in Fig. 10.

C. Limiting fragmentation
In Fig. 12 we investigate rapidity profiles of the differential transverse energy dE ⊥ /dη s for different collision energies √ s N N .After shifting the profiles by the respective beam rapidity Y beam = arcosh(γ beam ), we see remarkable agreement between LHC and RHIC setups regarding the shape of the flanks at large rapidities.In other words, for extreme rapidities (the fragmentation region) the rapidity distribution of the transverse energy is independent of the collision energy.This is an indication of limiting fragmentation, first observed experimentally for charged particle distributions in the context of proton-antiproton collisions [54], after a theoretical description was given in [55].It was subsequently studied experimentally e.g. for p+A collisions [56], d+Au collisions [57] and Au+Au collisions at RHIC [58,59].From our results, we expect limiting fragmentation to hold for LHC energies as well.We can confirm limiting fragmen- The results are averaged over 10 central collision events with error bands for one standard deviation.We find that the flanks of the rapidity profiles coincide irrespective of the collision energy, known as limiting fragmentation.
tation also for the local rest frame energy density ϵ LRF as well as T τ τ .Therefore, we conclude that it is truly a universal effect that governs the overall scaling behavior of the energy-momentum tensor in the dilute approximation.

D. Eccentricity
Up until now, we have highlighted the longitudinal structure of the three-dimensional Glasma, but our approach allows us to study the transverse structure as well, namely in the form of (transverse) eccentricity coefficients.From the local rest frame energy density ϵ LRF (τ, η s , x) we obtain the n-th eccentricity via the standard formula x γ ϵ LRF rn (68) with ϵ LRF = ϵ LRF (τ, η s , x) and γ = γ(τ, η s , x) = u τ (τ, η s , x) the local Lorentz factor.In the above, r = (x − x 0 ) 2 + (y − y 0 ) 2 and φ = arctan (y − y 0 )/(x − x 0 ) are polar coordinates in the center of mass frame in the transverse plane.We evaluate the transverse center of mass, at mid-rapidity, η s = 0. 0.0 0.2 In Fig. 13 we show the rapidity dependence of the eccentricities ⟨|ε 2 | 2 ⟩ and ⟨|ε 4 | 2 ⟩ for both RHIC and LHC energies and different combinations of ξ and m.We choose a large impact parameter of b = R and consequently see significant contributions to ε 2 and ε 4 .The eccentricities are independent of rapidity except for the most forward and backward bins, where they fall off.We note that the central plateau is narrower for RHIC than for the LHC setup.The fourth eccentricity ε 4 exhibits similar qualitative behavior to the second eccentricity ε 2 but has a lower overall value.In any case, ε 2 and ε 4 show only minor dependence on the longitudinal correlation length ξ.
Evidently, it would also be interesting to study the longitudinal decorrelation of the event geometry in the dilute Glasma by considering un-equal rapidity correlations of the eccentricities (see e.g.[17,60,61]).However, the fluctuations of the eccentricities ε n (η s ) are largely driven by fluctuations in nucleon positions, which are currently not included in our model and we therefore defer this analysis to a forthcoming study.
In Fig. 14 we show the real part of the first eccentricity ε 1 obtained from Eq. ( 68) for n = 1.Since the impact parameter lies in the x-direction we only get negligible contributions to the imaginary part.Figure 14 shows how the Glasma is skewed in the transverse plane.Specifically, for longitudinal fluctuations at the scale of the nucleon size, ξ = 0.1 R l , regions of positive rapidity feature more energy density for positive x.This is what one would expect from the collision geometry in the sense that nucleus B is shifted to positive x and moves in positive z-direction.Interestingly, in the coherent limit ξ = 2.0 R l , the opposite is the case and we see more energy in regions of positive x for negative rapidity.In our results, the energy distribution appears to depend strongly on the longitudinal structure of the nuclei-even changing signs for different ξ.
Generically, we find that the dilute (3+1)D Glasma exhibits both non-trivial longitudinal flow ⟨u η ⟩ ̸ = 0 and a skewed energy distribution in the transverse plane, ε 1 ̸ = 0.In [62] it was demonstrated, using a parameterized Glauber model, that longitudinal flow and skewness may both contribute to the directed flow of pions in off-central Au+Au collisions.We find that these features also arise naturally in the dilute (3+1)D Glasma.Thus, an interesting prospect would be to study the impact of the longitudinal structure on observables such as directed flow.
Finally, we remark that the RHIC and LHC curves' shapes are similar for large rapidities, but around midrapidity, a plateau becomes visible for LHC that is absent in the RHIC curves.

V. CONCLUSION
We presented the semi-analytical treatment of the three-dimensional, dilute Glasma in the weak field approximation.Starting from the expressions for color fields obtained in [32] within the dilute limit of the Color Glass Condensate effective field theory, we derived concise expressions for the field strength tensor and energymomentum tensor of the dilute Glasma.We found that our analytical results have a transparent geometric interpretation in terms of integrating over the causal past of the specific spacetime point at which we evaluate the energy-momentum tensor.
The main advantage of the dilute approximation is that it leads to analytic results for the Glasma field strength tensor.The corresponding integrals have to be solved numerically, but, compared to non-perturbative real-time (3+1)D lattice simulations [23][24][25][26], the numerical evaluation using Monte Carlo methods can be implemented much more efficiently than previous simulations at comparable spacetime volumes and enables system sizes that were not achievable previously.As such, we are able to describe in unprecedented detail the energy-momentum tensor and derived observables of the three-dimensional Glasma, accounting for both transverse and longitudinal structure within the colliding nuclei.Furthermore, our calculation of the field strength tensor does not require numerical time-stepping common to traditional simulation approaches.Instead, the field strengths can be evaluated at arbitrary points in the future light cone.This allows the direct evaluation of the energy-momentum tensor at arbitrary hypersurfaces and simplifies the coupling to hydrodynamics or kinetic theory.Especially for situations where many independent events have to be considered to obtain a given observable, we believe that the dilute Glasma offers a novel and economic approach to the pre-equilibrium stage of relativistic heavy-ion collisions.
By introducing a generalized McLerran-Venugopalan model with a fully three-dimensional structure for the color charge distributions of our nuclei we were able to parametrize longitudinal correlations within the nuclei.This led to a longitudinally extended overlap region for symmetric nuclear collisions and required us to shift the Milne frame origin into the future light cone outside of the overlap region.Within this setup, we explored the longitudinal and transverse structures of Pb+Pb and Au+Au collisions at LHC and RHIC energies.We carefully discussed how to interpret the Milne frame components of the energy-momentum tensor of the resulting Glasma.We found that the energy distribution of the Glasma exhibits broken boost-invariance and nontrivial longitudinal flow at large rapidities.Remarkably, we recovered limiting fragmentation as a generic feature of the dilute approximation.We believe that this can be demonstrated analytically, as will be detailed in a forthcoming manuscript.We also studied the pre-equilibrium, transverse structure in terms of eccentricity coefficients ε n .In particular, we find that off-central collisions lead to a rapidity-odd first-order eccentricity coefficient ε 1 .
There are multiple possible extensions to this work.As a direct application of our approach, the dilute Glasma can be used as fully (3+1)D initial conditions for effective kinetic theory [63][64][65][66] or hydrodynamic simulations of the QGP [2,19,67,68].Proceeding through these later stages of heavy-ion collisions would allow us to compute observables directly comparable to those measured in ex-periments and study their dependence on our fully threedimensional nuclear model.Furthermore, we are actively working on improving our nuclear model, for example by considering nucleons [25] or nucleon hot spots [69][70][71] as additional substructures within the nuclei.We are excited to explore the implications of such improved models on the various different aspects of the (3+1)D spacetime structure of the initial state of high-energy heavy-ion collisions in the future.To compute the dilute Glasma field strength tensor given in Eq. (24) we require first derivatives of the Glasma gauge fields a µ (x).For the f +− component, we start by working out which can be simplified by first rewriting the ∂ We dropped the boundary terms because φb B (x − − v − , q) is zero for finite x − if v − → ∞ and goes to zero if v − = 0 forming to Milne coordinates and normalizing the total velocity, we obtain a result u η tot (η s ) for fixed τ .We solve Eq. (B4) analytically and show the resulting curve in Fig. 10.

FIG. 4 .
FIG.4.Constant proper time hyperbolas for two different choices of the Milne frame origin separated by the distance δt.For the origin S (orange), the hyperbola enters the trajectories of the nuclei (gray).This is never the case for the Milne frame centered at S (blue).
FIG. 5. Slices through the center of a single component of a nucleus' color field for different values of infrared cutoff m and longitudinal correlation length ξ.Positive (negative) values are shown in orange (blue).

FIG. 6 .
FIG. 6. Perspective view of the local rest frame energy density ϵLRF of a single collision event at τ = 0.4 fm/c with impact parameter b = R. RHIC nuclear parameters with ξ = 0.5 R l , m = 0.2 µ are used.Lighter colors correlate with larger ϵLRF.A movie showing different viewing angles is included in the Supplemental Material [50].

3 ]FIG. 7 .
FIG. 7. Transverse (left) and longitudinal (right) slices through the local rest frame energy density ϵLRF of a single collision event at τ = 0.4 fm/c with impact parameter b = R. RHIC nuclear parameters with ξ = 0.5 R l , m = 0.2 µ are used.The white lines represent R and R + d boundaries of both nuclei.

ξFIG. 10 .
FIG. 10.Longitudinal flow u η with local rest frame energy density and integrated over the transverse plane as a function of spacetime rapidity ηs at τ = 0.4 fm/c.Discrete points represent the average of 10 central collision events.Lines are analytic predictions from simple model considerations (see Appendix B).

FIG. 11 .
FIG. 11.Differential transverse energy dE ⊥ /dηs as a function of spacetime rapidity ηs for different values of proper time τ .The points represent the average of 10 central collision events.The rapidity profiles stabilize at late times τ ≳ 0.6 fm/c.For collisions at LHC energy, we observe a plateau around midrapidity, which is absent at RHIC energy.

0 1 ξ 1 ξ
FIG.12.transverse energy as a function of spacetime rapidity ηs at fixed proper time τ = 0.4 fm/c.The curves are normalized to the respective LHC values at ηs = 0 and are offset by the respective LHC and RHIC beam rapidities Y beam .The results are averaged over 10 central collision events with error bands for one standard deviation.We find that the flanks of the rapidity profiles coincide irrespective of the collision energy, known as limiting fragmentation.

FIG. 13 .
FIG.13.⟨|ε2| 2 ⟩ and ⟨|ε4| 2 ⟩ at τ = 0.4 fm/c as a function of spacetime rapidity ηs.The impact parameter b = R and ⟨•⟩ refers to an average over 10 events with the error bands representing one standard deviation in the eventby-event statistics.

1 ξRe ε 1 FIG. 14 .
FIG.14.Real part of the eccentricity ε1 at τ = 0.4 fm/c for an impact parameter b = R.The points represent the average of 10 events and the error bands represent the standard deviation of the event-by-event statistics.

ACKNOWLEDGMENTS
ML, DM and KS are supported by the Austrian Science Fund FWF No. P34764.ML and KS further acknowledge funding from the Doktoratskolleg Particles and Interactions (DK-PI, FWF doctoral program No. W1252).SS acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the CRC-TR 211 'Strong-interaction matter under extreme conditions' -project number 315477589 -TRR 211.PS is supported by the Academy of Finland, by the Centre of Excellence in Quark Matter (project 346324), and European Research Council, ERC-2018-ADG-835105 YoctoLHC.The computational results presented have been achieved in part using the Vienna Scientific Cluster (VSC).
Appendix A: Derivation of the dilute Glasma field strength tensor

−
. Integrating by parts then gives

TABLE I
. Physical model parameters and their values in our calculations with R (L) denoting RHIC (LHC) setups.The correlation length ξ and the impact parameter b are given as multiples of the Woods-Saxon radius and are therefore different for RHIC and LHC setups.Without loss of generality, we always put the impact parameter in transverse x-direction.