Space-time structure of 3+1D color fields in high energy nuclear collisions

We perform an analytic calculation of the color fields in heavy-ion collisions by considering the collision of longitudinally extended nuclei in the dilute limit of the Color Glass Condensate effective field theory of high-energy QCD. Based on general analytic expressions for the color fields in the future light cone, we evaluate the rapidity profile of the transverse pressure within a simple specific model of the nuclear collision geometry and compare our results to 3+1D classical Yang-Mills simulations.


I. INTRODUCTION
The space-time evolution of the Quark-Gluon-Plasma (QGP), produced in high-energy heavy-ion collisions (HIC) at the Large Hadron Collider (LHC) and Relativistic Heavy Ion Collider (RHIC), can be accurately described by relativistic viscous hydrodynamics [1,2]. However, the early stage of HICs, which provides the initial conditions for the subsequent hydrodynamic evolution, still requires a comprehensive understanding. Over the course of years, sophisticated pre-equilibrium models incorporating fluctuations at nucleonic and sub-nucleonic level [3][4][5][6][7] and frameworks including intermediate kinetic theory evolution [8,9] have been developed in order to obtain a complete theoretical description in quantitative agreement with the experimental observations. Even though these models have played a significant role to improve the understanding of the initial state and properly characterize the subsequent QGP dynamics [10][11][12], the current state of the art modeling is often carried out at the level of an effectively 2+1D boost-invariant description and tends to ignore the longitudinal dynamics of heavy-ion collisions.
The development of initial state models has benefited from first principles insights into the initial stage of heavy-ion collisions provided by the Color Glass Condensate (CGC) effective theory of high-energy QCD [13,14], where small x partons in high energy nuclei are described by classical gluon fields, whose dynamics is governed by the classical Yang-Mills (CYM) equations. On the one hand, the CGC framework has lead to the development of various initial state models such as IP-Glasma [6,15] and MC-KLN [4,16], which along with relativistic viscous hydrodynamics have been successful in describing azimuthal anisotropies and charged hadron multiplicity [6,[17][18][19][20].
On the other hand, a plethora of (semi-)analytic calculations have been carried out within the CGC framework using expansions in color source densities [21][22][23][24] or near-field expansions [25][26][27][28][29] in the boost-invariant limit, which have been further exploited to study the correlation function of the initial energy momentum tensor [23,27,[30][31][32][33] and jet momentum broadening in the early stages [34,35]. Such numerical and (semi-)analytical results have also been important to guide the development of simple parametric initial state models such as IP-Jazma [36] or TrENTo [11,37], and to the development of a comprehensive understanding of the transverse dynamics of the fireball near mid-rapidity.
Beyond the boost-invariant description of heavy-ion collisions, recent experimental studies, e.g. of rapidity-dependent factorization breakdown [38][39][40], urge us to understand the dynamics of heavyion collisions beyond mid-rapidity, and have triggered an increased interest in the longitudinal structure of the initial state [41][42][43][44][45] and corrections to the eikonal limit [46,47]. While various implementations of 3+1D classical Yang-Mills equations have been developed either by varying the strengths of the classical sources [48], by generalizing the IP-Glasma model to 3+1D using JIMWLK rapidity evolution [49,50], or by taking finite thickness of colliding nuclei into account [51][52][53][54][55], we are not aware of any analytical calculation of the energy deposition in heavy-ion collision beyond the boost-invariant high-energy limit.
In this paper, we present the first analytical calculation of the initial energy deposition in heavyion collisions by solving the 3+1D classical Yang-Mills equations within the dilute limit of the Color Glass Condensate effective field theory of high-energy QCD. Similar to the previous studies [51][52][53][54][55], the longitudinal dynamics of the Glasma is analyzed by taking the finite extent of the colliding nuclei into account, while other corrections, such as the offset of the beam trajectory from the light cone, are neglected. We solve the linearized Yang-Mills equations and obtain analytic expressions for the perturbative gauge fields in the future light cone (analogous to the results of [23,24,56] for the boost-invariant case). Subsequently, we derive an analytic expression for the transverse pressure beyond the high-energy boost-invariant limit, for a specific realization of nuclear collision geometry.
We establish the effectiveness of our perturbative calculation by comparing it to non-perturbative 3+1D classical Yang-Mills simulations [53,55] for various thicknesses of the colliding nuclei.
This work is organized as follows: Starting in Section II, we set up the formalism to study 3+1D collisions in the dilute limit and develop an auxiliary field approach to obtain the analytic expressions for the color fields produced in the future light cone. We then employ a simple model of nuclear collision geometry to derive analytic expressions for the transverse pressure in Section III and compare our (semi-)analytic results to 3+1D classical Yang-Mills simulation in Section IV. We conclude with Section V.

II. GENERAL FORMALISM
The Color Glass Condensate (CGC) framework provides an effective description of nucleusnucleus collisions at high energies in terms of Yang-Mills theory. In the CGC formalism, hard partons of nuclei are modeled as color charges, which source soft partons in the form of classical color fields. Using light cone coordinates x ± = (x 0 ± x 3 )/ √ 2, the color current of a nucleus moving along the negative x 3 direction (denoted as "A") is given by where ρ a A denotes the color charge density per unit transverse area, x ⊥ = (x 1 , x 2 ) denotes the transverse coordinates, δ µ − is the Kronecker delta of the "−" light cone component and t a are the generators of the SU (N c ) gauge group. The color current depends only on one of the two light cone coordinates (in this case x + ) and is assumed to be localized around x + = 0. The color field A µ sourced by Eq. (1) is a solution to the Yang-Mills equations with the gauge covariant derivative and the non-Abelian field strength tensor given by In covariant gauge, ∂ µ A µ = 0, and using appropriate boundary conditions in the asymptotic past with all other components of A µ vanishing. The current and color field in Eqs. (1) and (5) solve the gauge covariant continuity equation Similarly, we can consider a nucleus moving along x − (denoted as "B") with the analogous current and color field In order to describe a collision of two nuclei using the Yang-Mills equations we need to solve the collision problem given by with initial conditions specified in the asymptotic past where we use calligraphic letters (A, J ) for the single nuclei solutions, while non-calligraphic letters (A, J) denote the solution for the collision problem (cf. discussion around Eq. (21)). In general, there are no closed form solutions for Eq. (9). However, in the ultrarelativistic limit where nuclei become infinitesimally thin, i.e.
where the color fields α i A/B are given by with the lightlike Wilson lines Equations (16) and (17) serve as initial conditions at τ = 0 + for the Yang-Mills equations in the future light cone for τ > 0, which are typically solved numerically on a lattice [15,[57][58][59] or through other approximations such as a Taylor expansion in proper time τ [25][26][27][28][29]. By construction, solutions in the boost-invariant high-energy limit do not depend on the space-time rapidity η.
To model collisions at finite energies, it is necessary to go beyond the boost-invariant approximation given by Eq. (15) and allow for a more general structure of the color charge densities which exhibit a non-trivial dependence on the light cone coordinates. Numerical solution methods using 3+1 dimensional real-time lattice simulations have been previously developed by the authors either based on the colored particle-in-cell method (CPIC) [51][52][53][54] or based on dynamically updated color currents [55]. Such methods rely on a direct solution of Eqs. (9) - (11) in the (t, z) coordinate frame and allow to numerically determine fully non-perturbative solutions for the Glasma in 3+1 dimensions. However, due to the large lattice sizes required for stable and accurate simulations, exploring realistic heavy ion collision scenarios using these methods is highly computationally demanding.
In this work we explore a different approach based on the weak field approximation to obtain semi-analytical approximations to Eqs. (9) -(11) beyond the boost-invariant limit. The weak field (or dilute) approximation is a perturbative expansion in the color charge densities ρ A/B (x ± , x ⊥ ) of the projectile and target. It relies on an explicit split of the gauge field A µ into background fields where the background fields A µ A/B and background currents J µ A/B are given by the single nuclei solutions Eqs. (1), (5) and Eqs. (7), (8). Evidently, the details of this perturbative expansion depend on the choice of gauge and we adapt covariant gauge throughout this paper. Since the background fields in Eqs. (5) and (8) We further note, that this gauge choice simplifies our calculation, because the covariant gauge background fields also simplify in covariant gauge. While the color field a µ describes the (dilute) Glasma itself, the currents j µ represent perturbations of the color currents J µ A and J µ B of nuclei A and B, due to non-Abelian color rotation. Expanding to the first non-trivial order O(ρ A ρ B ) in the color charge densities, the background field equations remain of the same form absorbing all terms of O(ρ n A ρ 0 B ) and respectively O(ρ 0 A ρ n B ), while the perturbative field equations that account for the interaction of the nuclei read with Since the perturbations represent the Glasma created from the collision of the two nuclei, we assume that both perturbative fields and currents vanish in the asymptotic past Assuming that the color charges of the colliding nuclei do not change their trajectories and considering the initial conditions in Eq. (29), the solution to Eq. (26) is straightforward: In covariant gauge, ∂ µ a µ = 0, Eq. (25) simplifies to where S µ (x) are a µ -independent source terms given by Here, we have introduced a slightly unusual but very useful shorthand Due to the choice of covariant gauge, we can independently solve for the four independent components of a µ in Eq. (32). Analyzing the x ± dependence of the source terms in Eqs.
where G ret (z) denotes the retarded propagator ensuring causality in compliance with the initial conditions.

A. Gauge field solutions in the future light cone
We now focus on carrying out the integration in Eq. (37) as far as possible to find simple expressions for the gauge field a µ in terms of color potentials φ a A/B of the colliding nuclei. We start by noting that the source terms in Eqs. (33) - (35) can be stated in a more unified way by performing a partial Fourier transform over transverse coordinates. We use with p ⊥ = d 2 p ⊥ (2π) 2 and rewrite the charge densities in terms of potentials viaρ a (x + , p ⊥ ) = p 2 ⊥φ a (x + , p ⊥ ).
Here, we use p ⊥ · x ⊥ = p i x i and p 2 ⊥ = p i p i . The source terms are then simply where we defined the auxiliary source term S d as Our strategy for performing integrations in Eq. (37) is most easily demonstrated using the transverse gauge field a i , where it is easy to see that the transverse integration only acts on the phase factor: By performing a change of variables z µ = x µ − y µ we can then solve the integral over y ⊥ as where cos φ = z ⊥ ·(p ⊥ +q ⊥ ) |z ⊥ ||p ⊥ +q ⊥ | is the azimuthal angle between z ⊥ and p ⊥ +q ⊥ and we denote |p + q| = |(p + q) ⊥ |. Evaluating the φ and |z ⊥ | integrals, we then obtain where τ z = √ 2z + z − and the two Heaviside functions imply that this term only contributes in the future light cone. By inserting this result into Eq. (44), we obtain which is more compactly written as with the auxiliary field a d given by Carrying out the same steps for a + and a − , we find analogous expressions and it is straightforward to check that Eqs. (48) -(51) satisfy the gauge condition We can further simplify the expressions for a + and a − by explicitly computing the derivatives and integrals with respect to x ± in Eqs. (50) and (51). Starting with the derivative term in Eq. (50), we use integration by parts to find The second line in the above expression is the boundary term for z − → 0, which generally does not vanish in contrast to the z − → ∞ boundary. However, it is proportional to the color potential φ b B (x − , q ⊥ ), which vanishes inside the future light cone. If we are only interested in far field solutions, we can safely ignore this term. By use of the following relations we then find where we use to denote that, due to the fact that we have ignored the boundary terms in the second line of Eq. (53), this expression is only strictly valid inside the future light cone. The term involving an integration in Eq. (50) is given by where we have performed a change of variables fromx + to z + = x + −x + +z + to isolate the integration overz + in the last line. 1 The integral over the Bessel function is given by 1 It is instructive to express the terms in the second line of Eq. (57) with two placeholder functions f (x + −z + ) and g(z + ) and re-formulate the integral bounds in terms of Heaviside functions as +∞ 0 By performing a change of variables, expressing z + = x + −x + +z + such that f (x + −z + ) = f (x + − z + ), one then finds that the Heaviside functions constrain the integration domain to 0 <z + < z + .
which leads to Inserting the results in Eqs. (56) and (59) into the expression Eq. (50), we finally obtain By repeating the same analogous steps for the calculation of a − , we obtain which together with the transverse components of the gauge fields in Eqs. (48) and (49) provide our final expressions for the Glasma fields in the future light cone.

III. NUCLEAR MODEL AND TRANSVERSE PRESSURE
Based on the previous analytical calculation, the longitudinal structure of the Glasma at late times can be obtained by considering a specific model for the color charge distribution inside a nucleus. Within this study, we consider a simple McLerran-Venugopalan (MV) type model [60,61] of a transversally homogeneous nucleus, where fluctuations of the color charge density are given by The constant g 2 µ A/B denotes the color density per unit transverse area which is related to saturation momentum Q s , while the function G characterizes the transverse correlation of color charges inside the nucleus. Similarly, the functions T R and U ξ describe the longitudinal profile and correlations of color charges, and are taken as normalized Gaussians with widths R and ξ identified as the Lorentz contracted size of the nucleus and longitudinal correlation length respectively. In order to enforce color neutrality on average, the one-point function is assumed to be zero.
Using this model, we can use our previous results to investigate a wide range of observables. In particular, we are interested in the various components of the energy-momentum tensor given by where we have made the split into background, mixed and perturbative terms explicit. For the purposes of this paper, we only consider the perturbative part of the energy-momentum tensor. We have derived the perturbative field a µ up to quadratic order, O(ρ A ρ B ), which yields the perturbative energy-momentum tensor up to quartic order, O(ρ 2 A ρ 2 B ). In principle, the mixed terms could also contain quartic contributions, however, the background field strengths F µν are only non-zero along the light cone. Consequently, the mixed terms vanish inside the future light cone to all orders.
Since we are only interested in the Glasma, we can safely ignore the mixed terms. In the following we focus on the transverse pressure which is solely generated during the collision and hence has no contribution from the background and the mixed part outside the space-time region where the colliding nuclei overlap. The transverse pressure is given by where ε E,L and ε B,L are the contributions from the longitudinal electric and longitudinal magnetic field given by

A. Longitudinal magnetic field
To get the longitudinal magnetic field, we first calculate the corresponding field strength f ij with Eq. (48) The square of the above expression contains integrals over four color potentials, arising from the auxiliary fields a d . Since it is quite convenient to solve such integrals in Fourier space, we write the correlation function in Eq. (63) as Exploiting the fact that the nuclear model is diagonal in momentum space, we have where we have evaluated the color factors as f abc f abc = N c (N 2 c − 1). To obtain the last equality, we have used Eqs. (5) and (8), and replaced the gauge field correlators with our nuclear model such that the overall transverse dependence is characterized by which we take as for both nuclei A and B, where, adopting the same conventions as in [55], m and Λ regulate the infrared and ultraviolet modes respectively. Inspecting the coordinate dependence of Eq. (70), it is convenient to perform a change of variables from z ± ,z ± to mean and relative coordinates Z ± , δz ± Since T and U are both Gaussian functions, we can change the limits of integration to The resulting expression for the longitudinal magnetic field is then given by with τ z = 2(Z + + δz + /2)(Z − + δz − /2) and τz = 2(Z + − δz + /2)(Z − − δz − /2).

B. Longitudinal electric field
Similarly, in order to calculate the longitudinal electric field, we start again with the associated field strength f +− by using Eqs. (50) and (51) Since we have already found the derivative of the auxiliary field a d in Eq. (56), we differentiate it again with respect to x + to get the first term of f +− . We find where the in the first line denotes the omission of boundary terms that are not relevant within the future light cone, and we used the Bessel identity to obtain the final equality. Using the Eq. (78) together with Eq. (49), we obtain With this, the expression for the longitudinal electric field takes the following form By combining the results in Eqs. (76) and (81), the resulting expression for the transverse pressure is given by This is the main result of this section which shows the dependence of the transverse pressure on the longitudinal structure of the colliding nuclei. We also note that by regularizing the color potential in the auxiliary source terms (see Eq. (43)) as reduces to the result for 2+1D [23,24,56] IV. NUMERICAL RESULTS AND COMPARISONS TO 3+1D SIMULATIONS Basic features of the reaction dynamics for 3+1D collisions have already been examined in detail using real time lattice simulations [53,55]. In this section we determine the effectiveness of our analytical calculation based on the weak-field approximation by comparing them with full 3+1D simulations. In the 3+1D simulations performed in Ref. [55], we defined the color charge density in Minkowski space as ρ a (t, z, x ⊥ ) = ρ a (2D) (x ⊥ )T Rγ (t + z) (for a left moving nucleus), which, due to the factorized x ⊥ and t + z dependence, leads to a less general model for the color charge densities Since the color charges are assumed to be x − independent, we can write the above two-point function in light cone coordinates as where Now by comparing Eqs. (63) and (85), one finds that the correlators can be matched by equating the factorized longitudinal dependence as By multiplying the two Gaussians on the left, we find that for ξ = 2R, the cross terms cancel and then the resulting relations are given as We can perform an analogous matching for the nuclear model used in Refs. [52,53]. Table I summarizes the parameters for the analytical results obtained from our perturbative expansion (Dilute) and two different 3+1D simulation schemes (3+1D CYM [55] and 3+1D CPIC [52,53]) with which we determine the extent to which the results of our weak-field approximations agree with the fully non-perturbative real time lattice simulations. Note that ξ = 2R, which we call the Dilute 3+1D CYM [55] 3+1D CPIC [52,53]  coherent limit, is the upper physical limit for the correlation length of color structures inside the nucleus.
We note that the non-linearity, which measures the strength of diluteness of a model, can be controlled by the dimensionless ratio of the color charge density g 2 µ and the infrared regulator m.
We vary this dimensionless parameter and compare the transverse pressure as obtained from the analytical result and the result from 3+1D simulations, for different longitudinal extents of the colliding nuclei in Fig. 1. We choose the same transverse lattice discretization for both schemes: where we use g 2 µτ = 2. 3 Before discussing the results of our weak field approximation, we emphasize that the results of the two different 3+1D classical Yang-Mills implementations (3+1D CYM [55] and 3+1D CPIC [52,53]) are in excellent agreement with each other. We find that, as per our expectation, the analytical calculation works remarkably well in the dilute limit g 2 µ/m = 0.5 as seen from the left panel of Fig. 1. By increasing the non-linearity of the model g 2 µ/m ≥ 1, we find that the analytical results in the dilute limit overestimate the transverse pressure; nevertheless the rapidity profiles are still reproduced rather well and the flattening of the rapidity profiles with increasing g 2 µ/m is correctly predicted by the (semi-)analytic calculation. 4 It is further interesting to note that the ratio of the analytical to that of the simulation results are roughly the same for 3 Note that for a proper comparison, we interpolate data on the (t, z) grids of the 3+1D simulations before switching to (τ, η) coordinates. To reduce statistical fluctuations, we use rapidity bins of width ∆η = 0.2. 4 It should be noted that we use the dimensionless length scale g 2 µR for the sake of comparing the weak field approximation to our non-perturbative simulations. However in the dilute limit, g 2 µR is in fact not a particularly useful scale, because g 2 µ only enters as an overall normalization factor of T µν . The more appropriate length scale is given by mR, i.e. when mR and m/Λ are fixed, the shape of the rapidity profile does not change with the non-linearity parameter g 2 µ/m. The widening of the dilute rapidity profiles in Fig. 1 for fixed g 2 µR and increasing g 2 µ/m should be interpreted as widening due to varying mR. different thicknesses of the colliding nuclei, which suggests that the non-linearity could effectively be introduced by re-scaling the pressure profile. Besides the significantly smaller computational cost, another enormous benefit of the semi-analytic calculation is that it is not bounded by lattice size, and therefore, we are able to perform computations at larger rapidities, as is clearly visible from Fig. 1.
Having established that our analytical expressions reproduce the full 3+1D numerical simulations in the dilute limit, we consider the various limits of our nuclear model. We start by looking at the coherence length ξ/R, which accounts for the randomness of color charges across a fixed longitudinal extent of the nucleus. Naturally, the longitudinal extent of the nucleus (R) is greater than the size of a correlated region within the nucleus (ξ) and hence for a physical limit ξ/R 1. In Fig. 2, we plot the transverse pressure for different thicknesses of colliding nuclei while considering two different values of ξ/R = 0.1 (left panels), which roughly corresponds to the ratio A −1/3 of the size of a nucleon and a nucleus in large nuclei, and ξ/R → 0 (right panels), which corresponds to the McLerran-Venugopalan (MV) model (A → ∞). For these plots we use m = g 2 µ and Λ/m = 5.
The top plots show the profiles for different g 2 µR at fixed time g 2 µτ = 5, whereas the bottom plots are evaluated for fixed g 2 µR = 1/8 varying g 2 µτ . We further include a comparison of the results of the 3+1D dilute calculation to the corresponding result in the 2+1D boost-invariant limit, which is obtained by integrating Eq. (83). We observe that the profiles approach the same boost-invariant plateau around mid-rapidity, whereas the flanks at larger rapidities are different and depend on the correlation length ξ/R. By decreasing the thickness of the colliding nuclei g 2 µR → 0, one approaches the boost-invariant limit, where the central plateau extends across larger and larger rapidity intervals. With regards to the proper time dependence, we find that in the limit ξ/R → 0 in Fig. 2 the profiles exhibit a significant time dependence up to very late times, where the flanks continue to move towards larger rapidities, while the central plateau remains time-independent. In contrast, we observe that the time dependence is much milder for ξ/R = 0.1 (Fig. 2 on the bottom left), where a stable profile is reached quickly and the double-peak structure vanishes entirely.
Finally, we investigate the dependence of the rapidity profile on the UV and IR regulators.
In Fig. 3 we plot the transverse pressure normalized to its value at mid-rapidity for different longitudinal extents g 2 µR and ξ/R = 0.1. In the left panel we fix the UV regulator to Λ = 5g 2 µ and vary the infrared regulator m to three different values. Similarly, for the right panel the IR regulator is set to constant m = g 2 µ, and Λ takes three different values. We observe that for a constant proper time g 2 µτ = 5, the profiles are largely insensitive to the variation apart from small deviations in the shoulders and flanks. Similar to Fig. 2, a boost-invariant plateau around mid-rapidity emerges upon decreasing the thickness of the colliding nuclei and the width of the plateau appears to be insensitive to the UV and IR regulators.

V. CONCLUSIONS AND OUTLOOK
We performed the first analytic calculation of the longitudinal profiles of the energy deposition in heavy-ion collisions within the dilute limit of the Color Glass Condensate effective field theory of high-energy QCD. We obtained general analytic expressions for the color fields of the Glasma produced in the future light cone (cf. Eqs. (60) -(62)), and employed them to study the rapidity dependence of the transverse pressure for a simplified nuclear model including non-trivial longitudinal color correlations.
By comparing the (semi-)analytic results in the dilute approximation to non-perturbative 3+1D classical Yang-Mills simulations, we confirm excellent agreement in the dilute regime. Even beyond the dilute limit, our approximation appears to capture the rapidity profiles rather well, while the overall magnitude of energy deposition is overestimated, once non-linear effects become important.
Since our analytic expressions allow for an efficient numerical determination of the energy momentum tensor T µν , the results presented in this paper provide new opportunities to further explore the longitudinal structure of matter produced in high-energy heavy-ion collisions, to study e.g. the interplay of longitudinal and transverse fluctuations and develop new Monte Carlo event generators for the initial state of heavy-ion collisions.