Effect of Longitudinal Fluctuations of 3D Weizsäcker–Williams Field on Pressure Isotropization of Glasma

We investigate the effects of boost invariance breaking on the isotropization of pressure in the glasma, using a 3+1D glasma simulation. The breaking is attributed to spatial fluctuations in the classical color charge density along the collision axis. We present numerical results for pressure and energy density at mid-rapidity and across a wider rapidity region. It is found that, despite varying longitudinal correlation lengths, the behaviors of the pressure isotropizations are qualitatively similar. The numerical results suggest that, in the initial stage, longitudinal color electromagnetic fields develop, similar to those in the boost invariant glasma. Subsequently, these fields evolve into a dilute glasma, expanding longitudinally in a manner akin to a dilute gas. We also show that the energy density at mid-rapidity exhibits a 1/τ decay in the dilute glasma stage.


Introduction
Relativistic heavy-ion collision experiments in nuclear physics aim to study various states governed by quantum chromodynamics (QCD), which describes the strong interaction within the standard model of particle physics.The collision deposits matter at extremely high energy and/or density in the collision region.Subsequently, this matter is believed to evolve into quark-gluon plasma (QGP), where quarks and gluons, normally confined within protons and neutrons in vacuum, are free to move independently.Recently, research has shifted from merely confirming the creation of QGP to exploring its properties through experiments.This includes investigating the effects of magnetic fields [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18] and rotation in the QGP medium [19][20][21][22], as well as the observability of novel phenomena related to quantum anomalies, such as the chiral magnetic effect [6,7] and chiral vortical effect [23][24][25][26].We would like to emphasize that heavy-ion collision experiments provide not only the QGP but also other rich physics landscapes of QCD, including the non-equilibrium state before the formation of QGP and the phase transition from the QGP to the hadron gas.
The pre-equilibrium process between the non-equilibrium state initially produced by the collision and the QGP following relativistic hydrodynamics remains unclear.Since it is not trivial whether the initial matter with high energy density reaches the QGP before undergoing the QCD phase transition, due to the decrease in energy density caused by the rapid expansion of the created matter, a deeper understanding of the pre-equilibrium process reinforces our confidence that QGP is produced in experiments.Not only that, it is also vital for explaining the experimental results.For instance, improved modeling of the pre-equilibrium process can provide more accurate initial conditions for the hydrodynamic simulations used in the analyses of the experimental data.Additionally, photons emitted from the pre-equilibrium matter may carry information about the non-equilibrium state, as they are less affected by the medium.Of course, the non-equilibrium state of QCD matter itself is also an interesting subject of study.
In the study of non-equilibrium processes in heavy-ion collisions, a key issue is understanding how the created matter eventually behaves as a fluid governed by hydrodynamics.If the system were confined within a box, it might be expected to eventually reach equilibrium over time due to the interactions among the quarks and gluons.However, in reality the system expands at nearly the speed of light.Therefore, it is not an obvious question whether or not the created matter can even approach a hydrodynamic state in the first place.
Tracking the real-time evolution of the non-equilibrium state of QCD matter based on first-principle calculation is not practical with our existing computational resources.Therefore, several effective descriptions have been employed to reveal its qualitative characteristics.For example, for the early stage of the collision, the classical Yang-Mills (CYM) field is an excellent description of the system, as explained in the next paragraph.
According to perturbative QCD, in a high-energy collision a nucleus can be considered as a dense state consisting of soft gluons.Here, 'soft' implies that the momentum carried is only a small fraction of the hadron's total momentum.The color glass condensate (CGC) effective theory is a suitable description of such a high-energy nucleus [27].In this theory, a high-energy nucleus is divided into soft partons (gluons) and hard partons (gluons and quarks with higher momentum than the soft gluons) [28][29][30].Soft gluons are described as the CYM fields emitted from hard partons expressed as classical color charge density.The name 'color glass condensate' has its origins as follows: 'Color' in QCD is a property analogous to electric charge.'Glass' is a disordered system and appears static in natural time.Similarly, hard partons are 'frozen' in time due to the time dilation effects of special relativity, compared to the lifetime of soft gluons, yet can be rearranged by interactions.'Condensate' refers to the dense and coherent soft gluons, which are treated as the classical field, akin to a Bose-Einstein condensate.The theoretical description of the CGC is analogous to the Weizsäcker-Williams approximation in electromagnetic dynamics.As a result of the collision between two CGCs, the initial state of the matter is expected to be a dense matter of soft gluons, a state referred to as 'glasma'.The glasma is well described by the CYM field, which is determined by solving the CYM equations of motion with two classical color sources, each representing the hard partons within different nuclei [31].
Here, a natural question arises regarding the extent to which the system approaches the hydrodynamic regime during the glasma stage.Many studies on this issue focus on the isotropization of its pressure because sufficient isotropy in the pressure is crucial for the system to obey hydrodynamics.As mentioned above, the pressure isotropy is driven by the growth of fluctuations that break boost invariance, due to instabilities of the CYM theory.However, previous studies have not demonstrated the pressure isotropy consistent with hydrodynamic fluid.For instance, some studies [53][54][55][56][57][58][59][60]66,67] have simulated the semiclassical time evolution of the glasma.These include fluctuations that break boost invariance due to minimal quantum uncertainty.They use classical statistical approximation (CSA) in a regime with much smaller coupling constants than realistic coupling constants.However, these simulations have not achieved sufficient pressure isotropy.Another study [68], using the JIMWLK equation-which is the evolution equation for the separation scale between soft and hard partons in the CGC effective theory-incorporates the initial condition of the glasma that breaks boost invariance.This initial condition is derived from the momentum distribution of soft gluons in a nucleus with finite momentum.These studies also do not show sufficient pressure isotropy.However, before concluding that the glasma does not exhibit the pressure isotropization under the classical approximation, it remains valuable to consider the effects originating from other physical mechanisms that break boost invariance.
In this study, we investigate the impact of boost invariance breaking, caused by spatial fluctuations of the classical color charge density along the collision axis, on the isotropization of pressure.This is done using the 3+1 dimensional glasma simulation without the boost invariance assumption, which has been studied recently [69][70][71][72][73][74].Boost invariance typically stems from the eikonal approximation for the classical color charge density, which is valid in the infinite momentum limit.The eikonal approximation includes three treatments of the classical color charge density or current: firstly, considering only the dominant component of the classical color current in the infinite momentum limit; secondly, treating the classical color charge density as infinitely thin along the collision axis due to Lorentz contraction; and thirdly, assuming the classical color charge density is not changed by the collision.However, in the 3+1D glasma simulation, the last two treatments are relaxed.The color charge density possesses a finite thickness and evolves dynamically according to the continuity equation.In this model, which accounts for color charge densities with finite thickness, there is scope to include longitudinal spatial fluctuations within the color charge density.The longitudinal spatial fluctuations are taken into account through the longitudinal correlation length in the classical color charge density, a parameter in the model.Although the physical origin of this correlation is still non-trivial at this point, it is interesting to conduct simulations with various correlation lengths and numerically investigate their impact on the evolution of pressure.
The research is significant as it delves into two critical areas: firstly, it is vital for comprehending how closely the glasma reaches a hydrodynamic fluid within the classical approximation, especially in relation to its impact on the evolution of the degree of the pressure isotropy.While pressure isotropy is a reliable indicator for tracking the thermalization of the system, it is crucial to acknowledge that, even when full pressure isotropy is achieved, the classical field is still far from quantum equilibrium.Classical field theory does not reach quantum thermal distribution but exhibits the Rayleigh-Jeans divergence in its equilibrium state.Therefore, it is not a suitable approximation for a quantum system nearing a thermal state, such as a hydrodynamic fluid.Secondly, it provides essential insights for developing more realistic models of the initial stages of the collision.
This article is organized as follows: In Section 2, we introduce the 3+1D glasma simulation method that we have recently developed [74].In Section 3, we present the results of both transverse and longitudinal pressure during the evolution of the 3D glasma, comparing pressure isotropization across various longitudinal correlation lengths inside the colliding nuclei.Section 4 is devoted to providing a summary and perspectives on our findings.

Method
In what follows, we shall explain the description of the 3D glasma that we have recently developed [74].We first provide a brief review of the CGC effective theory for a single relativistic nucleus.Subsequently, we discuss the understanding of the boost-invariant glasma created in the collision of relativistic nuclei using the eikonal approximation.Then, we introduce the 3+1D glasma simulation method, in which several aspects of the eikonal approximation are relaxed.

Cgc Effective Thoery for Single Relativistic Nucleus
In the CGC effective theory, effective for large nuclei at relativistic speeds, we treat 'hard' partons, gluons and quarks with higher momentum than soft gluons, and soft gluons separately.Hard partons are modeled as classical color charge density ρ, while soft gluons are treated as the CYM field A µ radiated from the classical color current J µ , determined by solving the CYM equations of motion with J, [D µ , F µν ] = J ν .A covariant derivative is defined as D µ ≡ ∂ µ − igA µ , and the field strength is defined by the commutator of D, F µν ≡ i[D µ , D ν ].In this context, only the dominant component of the classical color current in the infinite momentum limit is considered, J µ = δ µ+ ρ/g, where light-cone coordinates are introduced as x ∓ = (t ∓ z)/ √ 2. The dominance of the + component in the classical current is evident because the boost transformation in the positive z direction enhance this component, J = (J + , J − , J ⊥ ) → (e ϕ J + , e −ϕ J − , J ⊥ ), where ϕ is a rapidity.Now the continuity equation for the classical color current can be simplified, Under the gauge condition of A + = A − = 0, the continuity equation becomes ∂ + ρ = 0, and thus ρ should be static in the light-cone time The soft gauge field, a solution of the equation of motion, [D µ , F µν ] = J ν , is a functional of ρ and is known as where the index i runs over the transverse components, i = 1, 2, and V is the Wilson line, Here, ρ cov denotes the classical color charge density in the covariant gauge condition and is related to the color charge density in the A − = 0 gauge condition through the gauge transformation, The form of the classical color charge density in the CGC is not uniquely determined; instead, the weight function of the classical color charge density, W Y [ρ], is considered.The expectation value of a given observable that the soft gauge field possesses can be examined as its average under the weight function, which depends on the separation scale Y between hard partons and soft gluons.When the collision energy is not too high, the classical color charge density consists only of valence quarks.In the ideal situation, where a nucleus is extremely dense and extended infinitely in the transverse direction perpendicular to the collision axis, the weight function for such a color charge density consisting of valence quarks is provided as a Gaussian function [28][29][30]75], known as the McLerran-Venugopalan (MV) model, where (g 2 µ) 2 is the averaged square of the color charge density, which is a parameter controlling the magnitude of ρ in the MV model, and its integration, dx − (g 2 µ) 2 , should be taken in the same order of the squared saturation scale, Q 2 s , characterizing the CGC.It should be noted that the weight function in the MV model does not depend on the gauge choice, which is ensured by the trace in Equation ( 6).The expectation value of the single color charge density in the MV model vanishes, where the index a represents the internal degrees of freedom of SU(N), and the expectation value of two color charge densities in the MV model satisfies a relation, suggesting ρ at different spatial points are independent of each other, As the collision energy increases, ρ should incorporate contributions from an increasing number of particles with large momentum, which can be captured by solving the JIMWLK equation [76][77][78][79][80][81][82][83], the evolution equation for the separation scale, Y, with an initial condition based on the MV model.In the actual application for studying observables in the CGC, we randomly generate an event-by-event color charge density according to the weight function, W Y , and estimate the expectation values of observables as event average values.
Later, we consider the color charge density based on the MV model, which is relevant for the energy of Au-Au collisions at RHIC, with a center-of-mass energy per nucleon pair of 200 GeV.Effects beyond the MV model, which can be taken into account by solving the JIMWLK equation, become more important when analyzing larger energy regions at LHC.

3D Glasma in Collision of Relativistic Nuclei
We can describe the collision of two relativistic nuclei based on the CGC picture.In the case of two nuclei, the classical color current is described as the sum of the two classical color charge densities, each representing different nuclei, The color charge density ρ (1) represents a nucleus moving in the positive z direction, while the color charge density ρ (2) represents a nucleus moving in the negative z direction.
Let us show the solution of the CYM equation of motion and the continuity equation.Before the two nuclei make contact, due to causality, a solution for the gauge field is simply given as the sum of the solutions for each single nucleus, where However, the problem becomes much more complicated from the beginning time of the collision.After the collision happens, the CYM field includes the glasma emerging as a result of the collision.Amazingly, the solution representing the glasma is known under the eikonal approximation.This approximation includes three treatments of the classical color charge density or current: firstly, we keep the treatment considering only the dominant components of the classical color current, as in the single nucleus case, which is already assumed in Equation ( 9); secondly, we treat the classical color charge density as infinitely thin along the collision axis due to Lorentz contraction, ρ (1/2) ∝ δ(x ∓ ); and thirdly, we assume the classical color charge density is unchanged by the collisions, which means two classical color charge densities are static over the entire time, ∂ ± ρ (1/2) = 0.Under the eikonal approximation, a regular solution of the CYM equation of motion, [D µ , F µν ] = J µ , at the infinitely small proper time, τ = 0+, is explicitly known as [31] A i = lim i , A i ] , where the Fock-Schwinger gauge condition A τ = 0 is imposed, and the electric fields in this gauge are defined as E i ≡ τ∂ τ A i and E η ≡ ∂ τ A η /τ.Here, Milne coordinates are introduced as (τ, η) = x − , which are natural coordinates for describing the matter created in relativistic heavy-ion collisions.The relativistic nucleus exists around τ = 0, and the evolution of the created matter expanding along the collision axis at close to the speed of light can be naturally followed in terms of proper time, τ, instead of ordinary time, t.
Some important points about the solution for τ = 0+ under the eikonal approximation, shown in Equations ( 14) and (15), should be mentioned.First, this solution represents only the glasma, since the gauge fields representing the soft gluon, the Weizsäcker-Williams field, are located solely at τ = 0 in the limit of the infinite thin nucleus.The subsequent evolution of the glasma at a larger τ can be tracked by solving the CYM equation of motion, using the τ = 0+ solution as the initial condition.Second, this solution exhibits exact boost invariance, independent of rapidity, ∂ η A = 0 and ∂ η E = 0. Third, the color electromagnetic fields in this boost-invariant glasma orient along the collision axis, E η , B η ̸ = 0 and E i , B i = 0.These anisotropic color electromagnetic fields lead to the longitudinal pressure and the transverse pressure, defined as P ⊥ ≡ τ 2 T ηη and P ⊥ ≡ (T 11 + T 22 )/2, having the same absolute value but opposite signs, P ⊥ = −P L .
Recently, much attention has been paid to the 3+1D glasma simulation, which goes beyond the eikonal approximation.In the simulation methods proposed in previous papers [69][70][71][72][73][74], the eikonal approximation is relaxed, allowing the classical color charge density to possess a finite thickness and evolve dynamically according to the continuity equation.Among these methods, we employ the 3+1D glasma simulation method that we have recently developed [74].The distinctive feature of our 3+1D glasma simulation method is its use of Milne coordinates to track the evolution of the 3D glasma, whereas other simulation methods employ Minkowski coordinates.In what follows, we review our 3+1D glasma simulation method in Milne coordinates.
In the 3+1D glasma model, the color charge density is no longer assumed to be infinitely thin but has a finite thickness in the direction of the collision axis.Additionally, the color charge density dynamically evolves according to the continuity equation for each nucleus.Since it is not practical to solve a set of CYM equations of motion and continuity equations analytically without any assumptions, we numerically solve these equations in Milne coordinates with the initial collision conditions set for the proper time before the collision occurs.It is important to note that, in this method, we set the center positions of each nucleus in light-cone coordinates not around x ∓ = 0 but around positive finite values, c can be any value as long as they are sufficiently large.Consequently, the two nuclei do not yet make contact at the proper time of 0+.To distinguish these Milne coordinates from the usual ones, which are defined so that the centers of mass for each nucleus coincide at τ = 0, we introduce new notation, ( τ, η) = Accordingly, the definition of the usual Milne .
The initial condition of the modified proper time used in the 3+1D glasma simulation, τini , is set to be significantly earlier than the moment when the centers of mass for each nucleus coincide.Causality implies that the initial condition of the classical color charge density and CYM field at τ = τini is given by the sum of the solutions of the evolution equations for each individual nucleus.As shown in Equation ( 1), the classical color charge density before the collision is static in the light-cone time.Therefore, the initial classical color charge densities are the same as those in the infinite past, lim The way of giving the classical color charge density in the infinite past is explained in the next section.The initial CYM field is given by the sum of the solutions of the CYM equation of motion for each individual nucleus, as shown in Equation (2), where (A τ , A η ) are related to (A − , A + ) through the general coordinate transformation.This solution satisfies the Fock-Schwinger gauge condition for the modified proper time, A τ = 0. Accordingly, the electric fields in this modified Fock-Schwinger gauge condition are given by i i , However, in the actual setup for calculations performed later, we assume that the color charge densities before the collision have a Gaussian shape.Since the Gaussian function extends to infinity, two nuclei are correlated with each other, even in the infinite past.Therefore, the simple sum of the solutions for each single nucleus, as described above, is no longer a solution of the evolution equations and violates Gauss's law, In fact, we can make the electric fields given above comply with Gauss's law by adding the longitudinal electric field, the expression of which is the same as that for the boost-invariant glasma, The additional longitudinal electric field is negligibly small, E η ≪ E i , because the overlap of the tails of the Gaussian-shaped color charge densities is very minor at τ = τini .Finally, the initial condition of the classical color charge density and CYM field we employ is given by the sum of the solutions of the evolution equations for each individual nucleus with the small modification for the longitudinall electric field, i , A i ] .
In actual calculations, we place the CYM field and classical color current on a spatial lattice and track their evolution using their respective evolution equations.The detailed formulation of these is provided in our previous paper [74].

3D Color Charge Density
We assume the classical color charge density in the infinite past as the initial condition of the color charge density, ρ ini , under the premise that each nucleus is still largely unaffected by the other at the initial simulation time, τini .To simulate the effect of non-trivial longitudinal correlation in the classical color charge density, the initial condition, ρ ini , is modeled using the modified MV model.In this model, ρ ini is given to satisfy the following relation, the generalization of Equation ( 8) by replacing the delta function for x ∓ with a Gaussian function, where N(x ∓ , l ) denotes the normal distribution with a variance of l , which represents the longitudinal shape of nucleus in the x ∓ coordinate, and g 2 µ (1/2) becomes a function of the center position of the two color charge densities.The variances of the Gaussian functions, l (1/2) L , act as longitudinal correlation lengths.For simplicity, we further assume that the averaged square of the color charge density also has a Gaussian shape, where r (1/2) L are parameters that act as longitudinal extensions of a nucleus in the lightcone coordinates.In the actual simulations, we consider the collision of the same kinds of a nucleus, L = r L , and take x  28), the event-by-event color charge density satisfying the relation shown in Equation ( 27) as an event average can be obtained by using the event-by-event random number Γ, The random number Γ (1/2) satisfies the following event average, with To make σ L positive, l L should be smaller than 2r L .This upper bound of the correlation length seems to be consistent with the spatial expansion of the color charge density.The detailed derivation of the relation between Γ and ρ ini is given in our previous paper [74].
To avoid the infrared singularity of the Wilson line stemming from ∂ −2 ⊥ ρ, we introduce an infrared cutoff, Λ IR , by adding the following factor in front of the Fourier transform of the classical color charge density, where k ⊥ is the transverse momentum.The ultraviolet cutoff, Λ UV , is also introduced.The UV cutoff suppresses the high transverse momentum modes of the color charge density that exceed the limit of the CGC description.The color charge density represents the aggregate of color charges from hard partons within a transverse resolution scale, which is approximately the inverse of the transverse momentum.Therefore, UV modes originating from color charges in a small transverse spatial region are not so dense that the classical treatment is justifiable.Furthermore, the UV cutoff plays a role in facilitating the taking of the limit in continuous space.It is shown by an analytic calculation using the MV model that, without UV cutoff, a local operators of two gauge fields diverge in continuous space [84].The ultraviolet cutoff should be chosen as a larger value than the saturation scale, Q s , which is the scale focused on in the CGC.A detailed analysis of the UV cutoff dependence is beyond the scope of this paper.

Numerical Results
We show the numerical results for the 3D glasma, computed using the modified Milne coordinates.Our simulations are conducted using a lattice framework with a link variable, which is a gauge covariant and ensures that Gauss's law holds under their evolution equations [74].We utilize the leap-frog method to solve those equations set on the lattice.As for the boundary conditions, we employ periodic boundaries in the transverse directions, while ensuring that both the CYM field and classical color current are set to zero at the longitudinal direction's boundaries.The number of colors is 2, not 3.While calculations based on SU(2) cannot be used to make quantitative predictions for relativistic heavy-ion experiments, they are expected to lead to results qualitatively similar to SU (3).
In what follows, we first present the set of parameters involved in the model.Then, we provide the numerical results for the pressure and energy density at mid-rapidity, as well as the pressure results in the broader rapidity region.The energy density, transverse pressure, and longitudinal pressure are defined by the energy-momentum tensor in the 'usual' Milne coordinates, Here, they are integrated values over the transverse plane since we focus on their dependence on rapidity.Reflecting the fact that the energy-momentum tensor of the CYM field is traceless due to its conformal symmetry, the relationship between energy and pressure, ε = 2P ⊥ + P L , is held.Since the actual simulation is performed in modified Milne coordinates, the energy-momentum tensor to be calculated directly is in modified Milne coordinates T τ τ , T τ η and T η η .Therefore, the energy-momentum tensor in usual Milne coordinates T ττ and T ηη can be obtained by a general coordinate transformation, where xc = √ 2x c .To estimate the energy-momentum tensor from lattice simulations, we use the definition of the energy-momentum tensor on a lattice given in our previous paper [74].

Model Parameters
The set of model parameters for the initial 3D color charge density includes the following: the longitudinal extension of the classical color charge density, r L , the longitudinal correlation length in the classical color charge density, l L , the parameter controlling the averaged square of the color charge density, g 2 μ, and the infrared and ultraviolet cutoffs, Λ IR and Λ UV .To simulate the Au-Au collision with a center-of-mass energy per nucleon pair of 200 GeV at RHIC, we adjust √ 2r L to roughly match the radius of a gold nucleus divided by the gamma factor γ = 108 and the factor √ 2 stemming from the definition of the light-cone variables, r L = 0.1/Q s , where the saturation scale is set at Q s = 1 GeV.The parameter g 2 μ is set to ensure that the integration of the averaged square of the color charge density over the light-cone variable coincides with the squared saturation scale, z g 2 µ(x, x ⊥ ) 2 = Q 2 s .This setup for g 2 μ indicates that the event average of a given observable has a translational invariance in the transverse plane.Such a translational invariance is approximately established in the central region of large nuclei in the transverse plane, where the thickness of the nucleus does not vary largely.Therefore, our setup simulates the glasma in the central region, as depicted in Figure 1.The infrared cutoff is set at Λ IR = 0.2Q s , and the ultraviolet cutoff is taken to be double the value of Q s , Λ UV = 2Q s .It is important to note that in the classical approximation the choice of the coupling constant only provides the trivial change in the field variables, as A, E, J ∝ 1/g.This is because the coupling constant can be excluded in both the classical equations of motion and the continuity equation by redefining the scaled fields, A ′ = gA, E ′ = gE, and J ′ = gJ.The numerical results shown in the later sections are normalized by g to make them independent of g.We explore the effect of longitudinal correlation by comparing the 3+1D glasma under various correlation lengths, l L , ranging from 0.01r L to 2r L .The smallest value, l L = 0.01r L , which is significantly less than the nucleus radius, is thought to arise from the intricate structure within a nucleon.The largest value, l L = 2r L , which is twice the nucleus radius, is the maximum limit discussed in Section 2.3.It is crucial to emphasize that the physical origin of the longitudinal correlation is beyond the scope of this paper.
In addition, the following parameters should be inputted for the 3D glasma simulation.These parameters are chosen to be sufficiently safe values to minimize their impact on the physical results.The transverse lattice grid number and lattice spacing are set to N ⊥ = 56 and a ⊥ = 0.4/Q s , respectively.For the lattice in the modified rapidity direction η, the grid number and spacing are L η = 896, 7168 and a η = 2.6/L η , respectively.It is noteworthy that simulations with smaller l L require smaller longitudinal lattice spacings to accurately capture the structure of the color charge density in the longitudinal direction.The beginning modified proper time in the simulation and the center positions of the nuclei are taken as τini = 32r L and x c = 20 √ 2r L , respectively.All results in this section are estimated using 48 events.The error in the results is estimated by dividing the unbiased variance by the square root of the number of events.According to the central limit theorem, the unbiased variance aligns with the variance of the mean value in the infinitely large number of events.

Midrapidity Region
First, we investigate the behavior of pressure and energy density in the mid-rapidity region, η = 0. Figure 2 illustrates the evolution of glasma pressure at η = 0, with the pressure normalized to the energy density.The time range 0 ≤ Q s τ ≤ 9.7 covers the time before the onset time of the hydrodynamics, τ = 0.6 ∼ 1.0 [85], as predicted by the experimental analysis.The blue and red lines indicate pressures resulting from collisions of nuclei with two different lengths of longitudinal correlation in color charge densities, l L = 0.01r L and 2r L .The solid and dotted lines correspond to transverse and longitudinal pressures over the energy density, P ⊥ /ε and P L /ε, respectively.
Evolution of the transverse and longitudinal pressure normalized by the energy density, P ⊥ /ε and P L /ε, in the mid-rapidity, η = 0.The red and blue lines show the results for l L /r L = 2 and 0.01, respectively.The proper time is normalized by the saturation scale as Q s τ.
Despite exploring two distinct longitudinal correlations in the color charge densities, the observed pressures were qualitatively similar, though they showed quantitative differences, particularly in the early time region.Notably, at Q s τ = 0, transverse pressure is minimal compared to longitudinal pressure, reflecting the overlap of two nuclei centers at (τ, η) = (0, 0) and the dominance of the Weizsäcker-Williams field, a transverse electromagnetic field representing nuclei, which significantly contributes to longitudinal but not transverse pressure.Between Q s τ = 0 and 0.6, both transverse and longitudinal pressures undergo substantial changes, eventually acquiring opposite signs.The realization of the relation P ⊥ > 0 > P L indicates the formation of a strong longitudinal color electromagnetic field between separating nuclei, like the boost invariant glasma.From Q s τ = 0.6 to 2, the transverse pressure trends towards 0.5, while longitudinal pressure approaches zero, with no significant changes observed beyond τ = 2.These findings are in qualitative agreement with the previous results from the glasma calculations, where other effects of boost invariance breaking are taken into account [53][54][55][56][57][58][59][60]68].
Figure 3 illustrates the evolution of the glasma energy density ε at η = 0, normalized by that for Q s τ = 1, ε 0 , where the solid blue and red lines depict the energy density arising from collisions of nuclei with two different lengths of longitudinal correlation in color charge densities, l L = 0.01r L and 2r L .The time range is taken to be the same as in Figure 2. As in the pressure isotropization in Figure 2, despite the choice of l L , the results of the 3D glasma simulation are the same qualitatively at late time.Between Q s τ = 0 and 0.6, the energy density is dominated by the Weizsäcker-Williams field rather than the glasma.After Q s τ = 0.6, as the nuclei move past η = 0 the evolution of the glasma's energy density becomes apparent.A power function fitting with fitting parameters A and B, f fit (τ) = A/τ B , was employed to accurately determine the decay rate of the glasma energy density in the later stages.The obtained fitting parameter B, B = 1.06 ± 0.001 for l L = 0.01r L and B = 1.04 ± 0.0004 for l L = 2r L , indicates that the energy of the glasma at η = 0 decreases at almost the inverse of the proper time, 1/τ.This 1/τ decrease in the energy density at late time has been observed by previous 2D and 3D glasma calculations [53,55,61].The negligible longitudinal pressure, P ⊥ ≫ P L , and the decay of energy density, ϵ ∝ 1/τ, in the late stage, as observed above, are consistent with the system expanding longitudinally like a dilute gas.These observations significantly deviate from the behavior expected in an ideal gas, where pressure is isotropic, P ⊥ = P L = P, and the energy density decay is predicted by Bjorken hydrodynamics using the ideal gas equation of state, ε = 3P, ε ∝ 1/τ 4/3 .

Broader Rapidity Region
To deepen the understanding of the pressure isotropization in the 3D glasma, we examine the pressure results across a broader rapidity region.The Figure 4 provides both transverse and longitudinal pressures normalized by the energy density across the broader rapidity region, −2 ≤ η ≤ 2. The upper left and right panels display the transverse and longitudinal pressures for l L = 2r L , respectively, and the lower panels show results for l = 0.01r L in a similar arrangement.The plots with different colors show the results for the different proper times, Q s τ = 1, 2, 3, 4, and 5.It is found that around η = 0 both transverse and longitudinal pressures over the energy density are relatively flat, undergoing drastic changes in a larger rapidity range.To understand what contributes to the drastic change in behavior, we consider the evolution of the classical color current by examining its amplitude, defined as It can be seen in Figure 5 that ρ amp.vanishes in the rapidity region for flat P ⊥ /ϵ and P L /ϵ, shown in Figure 4. On the other hand, in the larger rapidity region, where P ⊥ /ϵ and P L /ϵ change significantly, ρ amp. is also found to increase significantly.This correspondence indicates that the former and latter regions should be identified as the glasma and the mixture of the glasma and Weizsäcker-Williams field.Therefore, the flat region in Figure 4 indicates that the pressure isotropization of the glasma at non-zero rapidity is almost the same as those at mid-rapidity.

Summary
We have investigated the effects of boost invariance breaking on the isotropization of pressure in the glasma.This breaking is caused by spatial fluctuations of the classical color charge density along the collision axis.Our calculation is performed using the 3+1D glasma simulation that we have recently developed.In this method, the longitudinal spatial fluctuations are taken into account through the longitudinal correlation length in the color charge density, a model parameter.While the physical origin of this correlation is still a subject of ongoing research, simulations with various correlation lengths offer useful insights into both the pressure evolution of the glasma and a deeper understanding of the 3D glasma model itself.
We have provided numerical results for pressure and energy density at mid-rapidity and in a broader rapidity region.Our findings show that, despite variations in longitudinal correlation lengths, the evolutions of the pressure isotropy observed are qualitatively similar.At Q s τ = 0, the transverse pressure is minimal compared to longitudinal pressure, indicative of the dominance of the Weizsäcker-Williams field.Then, at the early time, the positive transverse pressure and negative longitudinal pressure appear, which indicates the formation of strong longitudinal color electromagnetic fields between separating nuclei like the boost invariant glasma.Finally, the transverse pressure becomes almost half of the energy density, and longitudinal pressure becomes negligible at the late time.Further, it is found that, at this late time, the glasma's energy density at mid-rapidity decreases as 1/τ.Both the negligible longitudinal pressure and the 1/τ decay of energy density in the late stage are consistent with the system expanding longitudinally like a dilute gas, deviating significantly from the Bjorken hydrodynamics with the ideal gas equation of state.By examining the pressure results across a broader rapidity region, we further enhance the understanding of pressure isotropization.The pressure isotropization of the glasma in the non-zero rapidity region is almost the same as that observed in mid-rapidity.
Within the scope of this study, fluctuations originating from the correlation length of the color charge density do not facilitate the isotropization of pressure in the glasma.This suggests that a framework incorporating effects beyond the classical approximations, such as kinetic theory, is required to explain the hydrodynamization of the matter produced in heavy-ion collisions.However, note that our results do not imply that the setting of the correlation length is insignificant in the 3D glasma model.In fact, the degree of pressure isotropization in the early stage and the magnitude of the normalized energy density are found to be dependent on the correlation length.Therefore, selecting an appropriate correlation length is likely essential for conducting more realistic 3D glasma simulations.
In addition, as a potential area for future research, it would be intriguing to explore how the longitudinal correlation length of the color charge density is translated into rapidity correlations in experimental observables, such as elliptic flow.
We discuss the interpretation of the longitudinal correlation length in our 3D nucleus model.It is roughly estimated as the inverse of the typical longitudinal momentum of the hard partons contributing to the color charge density.This length may be associated with the Bjorken variable x, which in turn implies that the correlation length is also related to the saturation scale Q s (x).Addressing the challenge of advancing beyond a basic Gaussian model to more precisely depict the 3D color charge density for a specific Bjorken x is important.The JIMWLK equation offers a framework for calculating observables from the color charge density over a wide range of Bjorken variable x but assumes an extremely thin spatial distribution of color charge density.Then, to accurately determine the spatial rapidity dependence of an observable, a conversion from momentum rapidity to spatial rapidity is essential.Therefore, the feasibility of constructing a spatial 3D color charge density based on the JIMWLK equation presents a complex challenge.
Finally, we point out the gap between the current 3D glasma model and phenomenology.To enhance the realism of the initial conditions, it would be advantageous to utilize the methodology from the 2D IP-glasma model [86,87], a phenomenological model of the glasma in the mid-rapidity.This approach allow us to account for the event-by-event fluctuation of nucleon positions, the saturation scale based on the IP-sat model [88,89], and the expansion of the system in the transverse plane, all of which are not taken into account in the current model.For the transition of glasma into initial hydrodynamic conditions, separating the glasma from the Weizsäcker-Williams(WW) field is crucial.One proposed method is to employ new "usual Milne coordinates" for observation, setting the proper time to 0 at the nucleus's rear surface, unlike the current usual Milne coordinates, which link proper time of 0 with the nucleus's central longitudinal position.This alteration would place the majority of the WW field outside the forward light cone, facilitating the separation from the glasma.Considering the energy density and pressure in the local rest frame, defined as eigenvalues of the energy-momentum tensor, could potentially exclude the contribution of the WW field on the EM tensor [90].Since the EM tensor's determinant for the WW field is zero, energy density and pressure in the local rest frame should effectively disappear in regions heavily influenced by the WW field.To further align the 3D glasma model with real-world phenomena, transitioning from SU(2) to SU(3) is necessary.

c
= x c for simplicity.Under the setup shown in Equation (

Figure 1 .
Figure 1.Our setup simulates the system in the central region of nuclei in the transverse plane.

Figure 3 .
Figure 3. Evolution of the energy density ε in the mid-rapidity, η = 0, normalized by that for Q s τ = 1, ε 0 .The red and blue lines show the results for l L /r L = 2 and 0.01, respectively.The proper time is normalized by the saturation scale as Q s τ.

Figure 4 .Figure 5 .
Figure 4.The transverse and longitudinal pressures in the rapidity region, −2 ≤ η ≤ 2, for Q s = 1, 2, 3, 4, and 5.The upper left and right panels show the transverse and longitudinal pressures for l L /r L = 2, respectively.The lower left and right panels show the transverse and longitudinal pressures for l L /r L = 0.01, respectively.The plots with different colors show the results for the different proper times, Q s τ = 1, 2, 3, 4, and 5.