Physical characteristics of glasma from the earliest stage of relativistic heavy ion collisions

We present analytic results that describe the gluon field, or glasma, at very early times after a collision of relativistic heavy ions at proper time $\tau=0$. We use a Colour Glass Condensate approach, and perform an expansion in $\tau$. The full details of our method are described in our previous paper [1]. In this paper we present an analysis of various physical quantities that can be obtained from the energy-momentum tensor. We show that the expansion to order $\tau^6$ can be trusted to about $\tau=0.05$ fm. For times small enough that the expansion converges, the transverse and longitudinal pressures move towards their equilibrium values of one third of the energy density. The Fourier coefficients of the azimuthal flow are larger than expected, which contradicts the usual assumption that anisotropy is mostly generated during the hydrodynamic evolution of the plasma. We find a significant correlation between the elliptic flow coefficient and the eccentricity, which indicates that the spatial asymmetry introduced by the initial geometry is effectively transmitted to the azimuthal distribution of the gluon momentum field, even at very early times. This result is interesting because correlations of this kind are characteristic of the onset of hydrodynamic behaviour. We show that the angular momentum of the glasma is orders of magnitude smaller than the angular momentum of the initial system of ions colliding with non-zero impact parameter. This indicates that most of the angular momentum carried by the valence quarks is not transmitted to the glasma, and contradicts the picture of a rapidly rotating initial glasma state that has been proposed by several authors, but agrees with the current lack of experimental evidence for a significant polarization effect of the hyperons and vector mesons produced in heavy ion collisions at the highest accessible energies.

assumption that anisotropy is mostly generated during the hydrodynamic evolution of the plasma.
We find a significant correlation between the elliptic flow coefficient and the eccentricity, which indicates that the spatial asymmetry introduced by the initial geometry is effectively transmitted to the azimuthal distribution of the gluon momentum field, even at very early times. This result is interesting because correlations of this kind are characteristic of the onset of hydrodynamic behaviour. Finally, we have calculated the angular momentum of the glasma and obtained results that are many orders of magnitude smaller than the angular momentum of the initial system of colliding ions in a configuration with non-zero impact parameter. This indicates that most of the angular momentum carried by the valence quarks is not transmitted to the glasma. The result is significant because it contradicts the picture of a rapidly rotating initial glasma state that has been proposed by several authors, but agrees with the current lack of experimental evidence for a significant polarization effect of the hyperons and vector mesons produced in heavy ion collisions at the highest accessible energies.

I. INTRODUCTION
We use a formulation of the Colour Glass Condensate (CGC) effective theory to describe the dynamics of a heavy-ion collision in the early stages after the collision (τ 1 fm).
Many reviews of the CGC effective theory have been published, see for example [2,3].
The details of the evolution of a quark-gluon plasma (QGP) during this early stage are not well understood, but they are important because they provide the initial conditions for subsequent hydrodynamic evolution. The CGC approach is based on a separation of scales between valence partons with large nucleon momentum fraction, and gluon fields with small nucleon momentum fraction. When the separation scale is fixed, the dynamics of the gluon fields can be determined from the classical Yang-Mills (YM) equation with the source provided by the valence partons, by averaging over an ensemble of valence parton colour charge distributions.
This paper is a companion paper to a previous work [1] in which we have explained the strategy of our approach and some details of the calculational method. We use an expansion in the proper time, also called a 'near field' expansion [4][5][6][7][8]. We work with infinitely Lorentz contracted sources, and the description of the system is classical. An advantage of the method however is that our results are analytic, and thus provide a potentially valuable alternative approach to the various numerical methods that are in use. In Ref. [1] we focused on the technical details of the calculation and showed only a few results, which were obtained for the simple case of nuclei that are infinite in the transverse plane and uniform.
In this paper we consider more physically realistic collisions where the nuclear area density function is not assumed constant. We present results from several different calculations and discuss their connection to experimental observables.
In section II we define some notation and give a summary of the results of our previous paper. In section III we describe the structure of the energy-momentum tensor in Milne coordinates and, by exploiting the symmetries of the tensor, we give a fairly compact analytic expression for our result to order τ 2 . The order τ 4 results are given in Appendix B. We also formulate a calculation of the angular momentum of the glasma per unit rapidity, on a hypersurface of constant proper time. In section IV we present some numerical results and discuss their significance in the context of heavy-ion collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC). In section IV A we define our notation for the geometry of the collision and the units we will use. In section IV B we discuss our method to implement non-constant nuclear area density functions, which is based on a gradient expansion of a two dimensional projection of a Woods-Saxon distribution. We follow the method of [7] and discuss carefully the limitations of the expansion and how its convergence depends critically on what quantity is calculated. In sections IV C and IV D we study the energy density and pressure of the glasma. For a system in equilibrium the transverse and longitudinal pressures are equal to each other, and individually equal to a third of the energy density. A calculation of the pressure therefore gives information about how far the glasma is from equilibrium. We show that, for times small enough that the expansion converges, the transverse and longitudinal pressures move towards their equilibrium values. We calculate a quantity that describes the asymmetry between the transverse and longitudinal pressures, and one that characterizes the difference between the two components of the transverse pressure. We show that the first is almost completely insensitive to the gradient expansion, but the second shows strong dependence. In section IV E we study the flow of the energy of the gluon field by calculating the radial projection of the Poynting vector. Our results in this section, and the previous two sections, indicate that the expansion to order τ 6 can be trusted to about τ = 0.05 fm. In our previous paper [1] we have given a simple argument that this time is much greater than the lower bound at which we no longer trust the classical description that is inherent in our approach. In section IV F we study the momentum anisotropy of the glasma by calculating the Fourier coefficients of the azimuthal distribution of the flow. Anisotropy is an important characteristic because it is primarily sensitive to the properties of the system very early in its evolution. It is expected that the spatial anisotropy of the system that is produced at very early times (within the first few fm) will be encoded in the observed anisotropies of the final particle momentum distributions. Our results for the Fourier coefficients are larger than expected, which does not agree with the common assumption that momentum anisotropy is mostly generated during the hydrodynamic evolution of the plasma. In section IV G we consider the spatial azimuthal asymmetry of the glasma by calculating the eccentricity. We look for correlations between spatial asymmetry and momentum anisotropy, which is characterized by the Fourier coeffients calculated in the previous section. Since momentum anisotropy originates in the initial spatial asymmetries in the geometry of the system, such correlations provide information about the effect of the interactions, and the expansion of the system, during glasma formation. We find significant correlation between the eccentricity and the elliptic flow coefficient calculated in the previous section, which indicates that the spatial asymmetry of the initial energy density is converted into the anisotropy of the azimuthal distribution of the gluon momentum field, and this correlation is characteristic of hydrolike behaviour. In section IV H we look at the angular momentum of the glasma, which develops in collisions with non-zero impact parameter, in the direction perpendicular to the reaction plane. It has been proposed that because of spin-orbit coupling effects, the angular momentum of the glasma could lead to the polarization of produced quarks and anti-quarks, which might be detected by measuring the polarization of the Λ hyperon, or various vector mesons [9]. Measurements of this kind could be used to compare different hydrodynamic models and hadronization scenarios [10]. Our results are many orders of magnitude smaller than the initial angular momentum of two ions colliding with non-zero impact parameter.
Our findings therefore disagree with the proposal that the glasma acquires a large fraction of the angular momentum of the participating valence quarks [10], but is not contradicted by any experimental evidence, since attempts to measure the polarization of produced hadrons have found only very small effects [11,12]. In section V we conclude with some discussion and comments on possible future directions of this research.

II. SUMMARY OF PREVIOUS WORK
A. Energy-momentum tensor In a previous paper [1] we calculated the energy-momentum tensor using the CGC effective theory, to sixth order in an expansion in the proper time. In this section, we review the main elements of that calculation.
We consider a collision of two ions moving towards each other along the z-axis and colliding at t = z = 0. In the pre-collision region of space-time, the system is most naturally described using light-cone coordinates, but in the post-collision region Milne coordinates are more convenient. At the end of the calculation, when we look at physical quantities, we usually want to use Minkowski coordinates. The calculation is done most efficiently by using the three different coordinate systems at different stages, and transforming between them as needed. We will use greek letters for 4-indices, and latin letters for 2-indices denoting coordinates in the transverse plane. For example, the components of a position 4-vector in Minkowski coordinates (t, z, x, y) are written x µ with µ ∈ (0, 1, 2, 3), and the transverse components are denoted x i with i ∈ (2, 3). In light-cone and Milne coordinates, the time and longitudinal coordinates are defined as functions of the Minkowski variables (t, z), but the transverse coordinates are the same in all three bases. Our notation for light-cone and Milne variables is standard, and is reviewed in Appendix A.
The vector potential of the gluon field is described with the ansatz [13] A The functions β i 1 (x − , x ⊥ ) and β i 2 (x + , x ⊥ ) represent the pre-collision potentials, and the functions α(τ, x ⊥ ) and α i ⊥ (τ, x ⊥ ) give the post-collision potentials. In the forward light-cone the vector potential satisfies the sourceless YM equation. We find solutions valid for early postcollision times by expanding in the proper time τ [4][5][6][7][8]. Using these solutions we can write the post-collision field-strength tensor, and energy-momentum tensor, in terms of the initial potentials α(0, x ⊥ ) and α ⊥ (0, x ⊥ ). The initial potentials are related to the physical properties of the ions and the geometry of the collision using boundary conditions that connect the pre-collision and post-collision potentials. These conditions are found by integrating the YM equation across the light-cone [1,13] where the notation lim w→0 indicates that the width of the sources across the light-cone is taken to zero (the pre-collision potentials depend only on transverse co-ordinates in this limit). Using equations (1, 2) the energy-momentum tensor can be written in terms of the pre-collision potentials β 1 (x − , x ⊥ ) and β 2 (x + , x ⊥ ) and their derivatives.
The next step is to use the YM equation to write the pre-collision potentials in terms of the colour charge distributions of the incoming ions. One then averages over a Gaussian distribution of colour charges within each nucleus. The average of a product of colour charges can be written as a sum of terms that combine the averages of all possible pairs, which is called Wick's theorem. The average of a product of pre-collision potentials, which depend on the colour charges in a nontrivial way, is much more difficult to calculate, and the calculation becomes more and more complicated as the number of potentials increases [14][15][16][17][18]. We use the Glasma Graph approximation [18], which was also used in other near field expanded calculations, and is equivalent to the application of Wick's theorem to lightcone potentials directly. In our previous work [1] we showed that for the simple case of homogeneous ions that are infinite in the transverse plane, there is some evidence that the effect of this approximation is small. We stress however that the range of validity of the Glasma Graph approximation has not been carefully studied, and this is an open and important issue.
In the context of our calculation, the use of the Glasma Graph approximation means that the energy-momentum tensor can be written in terms of the 2-point correlators of the precollision potentials. This correlator was originally calculated in Ref. [19] and generalized to include some effects of nuclear structure in Ref. [7]. We give below the result in our notation, see [1] for further details.
The correlator of two potentials from different ions is assumed to be zero. The 2-point correlator for two potentials from the same ion is written in terms of the colour charge surface density for that ion, which we denote µ 1 ( x ⊥ ) and µ 2 ( x ⊥ ). These functions are not determined by the CGC model, rather one assumes some form, constrained by experimental knowledge. We use a 2-dimensional projection of a Woods-Saxon distribution which is characterized by three parameters that correspond to the surface thickness, radius, and displacement of the centre of the ion relative to the beam axis (for details see Section II B).
The 2-point correlators for the two ions have the form and using the index n ∈ {1, 2} to represent the two ions we have where The function K 0 is a modified Bessel function of the second kind, and m is an infrared regulator. In the limit m → 0 the Green's function G( x ⊥ ) goes like ∼ ln(m| x ⊥ |). Since valence parton sources come from individual nucleons, the Green's function should go to zero when | x ⊥ | approaches the confinement scale, and therefore we choose m ∼ Λ QCD .

B. Parton surface density
To make further progress we must specify the form of the colour charge surface density of the nuclei. We will use a 2-dimensional projection of a Woods-Saxon potential of the form The parameters R A and a give the radius and skin thickness of a nucleus of mass number A, and their numerical values are discussed in section IV B. The integral in (8) is normalized so that for a lead nucleus µ( 0) =μ, which is sometimes called the McLerran-Venugopalan (MV) scale. This parameter is related to the saturation scale Q s , but its exact value cannot be determined within the CGC approach (for a discussion see [20]). We useμ = Q 2 s /g 4 , and some motivation for this choice can be found in [1]. Due to the ambiguity associated with the value of the MV scale, our numerical results for quantities like the energy density and pressure should be regarded as order of magnitude estimates. Quantities that depend on ratios of different elements of the energy momentum tensor, like Fourier coefficients of the azimuthal flow, will have much weaker dependence on the MV scale.
We obtain an analytic result for the energy-momentum tensor by substituting equation (8) into equation (6) and performing a gradient expansion, using the method developed in Ref. [7]. The coordinates x ⊥ and y ⊥ are rewritten in terms of relative and average coordinates. To consider collisions with non-zero impact parameter we expand the distribution µ 1 ( z ⊥ ) around We will keep terms up to second order in gradients of the distribution. The parameter that we assume to be small is where the gradient operator indicates differentiation with respect to the argument of the function. The region of validity of this expansion is discussed in section IV B.
We remind the reader that for a realistic nucleus, which is made up of individual partons, the transverse charge distribution is not a very smooth function. It is possible that the transverse charge distribution of a real nucleus could be sufficiently irregular that a Woods-Saxon distribution is not a good representation.

C. 2-particle correlators
In this subsection we drop the subscript that indicates which ion is being considered, and we set b = 0. Performing the derivative expansion and keeping terms up to second order in gradients of µ, equation (6) becomes We can rewrite equation (10) in the form where we have made the replacementr irj → δ ij /2 because in the limit r → 0 we knowγ must be independent of the direction of the vectorr. We note we are able to make this replacement before performing any derivatives with respect to x ⊥ and y ⊥ , since lim r→0 ∂ i x . . . ∂ j y . . .r krl = 0, where the dots indicate any number of derivatives.
The correlator B ij ( x ⊥ , y ⊥ ) and its derivatives have ultra-violet divergences that must be regulated. We use a modified version of the method proposed in Ref. [6]. To illustrate how we make use of equation (11) we consider, for example, the calculation of which appears in the expression for B ij in equation (4). The integration over angular variables gives k i k j → δ ij k 2 /2. The second term in (12) is finite, but the first term is logarithmically divergent and we regulate it using an ultraviolet momentum cutoff Λ. This cutoff will be set to the saturation scale, see our previous paper [1] for some discussion of this point. In section IV D 1 we test the dependence of some of our results on the saturation scale, and show that they are not very sensitive to its numerical value. Now we consider the contribution from the factor in round brackets in equation (4). Ex-panding this factor we have When we calculate derivatives of the correlator B ij ( x ⊥ , y ⊥ ), the derivatives operate on all terms in the expansion in equation (13). At sixth order in the τ expansion the energymomentum tensor includes terms with six derivatives acting on the correlator in equation (4). Naively it would seem that we need to expand the exponent in equation (13) to seventh order, since each of the six derivative operators will have a piece proportional to ∂/∂r i which could act separately on each of the six factors in the term (∆γ( x ⊥ , y ⊥ )) 6 . For example, if we differentiate six times with respect to r 1 we obtain an expression of the form where the dots represent additional terms that give zero when r is taken to zero. However, it is easy to see from equation (11) that if we differentiateγ an odd number of times with respect to r 1 or r 2 , the integration over momentum variables gives zero. This means terms with more than three factors of ∆γ, which are operated on with a maximum of six derivatives with respect to components of r, can be set to zero. Equivalently, we have to expand the exponential only to fourth order.
All correlators and their derivatives can be obtained using the method described above.
We give one example

A. Structure of the energy-momentum tensor
In this section we present our analytic result for the energy-momentum tensor, to order τ 6 . For simplicity of notation we give our results for the tensor in Milne coordinates, where there is no dependence on rapidity. All elements in the energy-momentum tensor have either even or odd powers of τ . We summarize this information in the symbolic equation The top left element of the tensor shows that the element T 00 has contributions of order τ 0 , τ 2 , τ 4 and τ 6 . To give one other example, the entry in the top right corner shows that the element T 03 has contributions of order τ , τ 3 and τ 5 .
We do not need to give all components of the energy-momentum tensor because of its symmetry properties. We will give our results for the six elements in equation (17) that are written in boldface. All of the other elements in the tensor can be generated from these using symmetries, as explained below. We emphasize that we have calculated all elements of the energy-momentum tensor, and the symmetries summarized in equation (17) have been verified, and not assumed.
All elements below the diagonal can be related to elements above the diagonal using the fact that the tensor is symmetric. In equation (17) we write, for example, T 10 = T 01 . There are also pairs of elements where one can be obtained from the other using the transformation ∇ x ↔ ∇ y where we have defined ∇ x ≡ ∂/∂R x and ∇ y ≡ ∂/∂R y . We will also use ∇ 2 ≡ ∇ 2 x + ∇ 2 y . For example, we write T 13 = FT 12 to indicate that the element T 13 can be obtained from T 12 by interchanging the derivative operators ∇ x and ∇ y . Finally, since the energy-momentum tensor is traceless it satisfies g µν T µν = 0, which means that we do not have to give all elements on the diagonal. We replace the element T 11 with the symbol Tr to indicate that this matrix element can be constructed from the other diagonal elements of the tensor. Combining this information we write In summary, equation (17) tells us that we need to give only the six elements of the energymomentum tensor that are written in boldface.

B. Coefficients of the energy-momentum tensor
We give generic equations for these elements organized by powers of τ and numbers of derivatives with respect to each of the transverse coordinates. Our result for the energymomentum tensor has the form For each of the greek letter variables in equation (18) the subscript gives power of τ that multiplies the variable when it is not acted on by any external derivatives, and the superscripts give the number of internal derivatives with respect to R x and R y . For example, the second term in T 12 contains the variable β 10 2 which has coefficient τ 2 and is defined (see equation (20)) as a sum of terms each of which has one derivative with respect to R x .
We use X to denote any of the greek letters {E, β, γ, δ, ξ}. We give below our results for X lm n with 0 ≤ n ≤ 2. In Appendix B we give our results for X lm n with n = 4. Each symbol is either even or odd under the transformation µ 1 ↔ µ 2 , and therefore we only need to give half of the terms, and the sign of the symmetry. We set the ultraviolet cutoff Λ equal to the saturation scale Q s and introduce the notation L ≡ ln(Q s /m). To save space we give only the contributions that are leading-order in an expansion in the infrared mass m, and this power counting is done with the assumption (∇ x ) n i (∇ y ) n j µ 1 ( R)/µ 1 ( R) ∼ m n i +n j (and similarly for µ 2 ( R)). We will also factor the constantμ ≡ Q 2 s /g 4 out of the source density functions and defineμ 1 ( R) ≡ µ 1 ( R)/μ, and similarly for ion 2. To save space we will use the shorthand notationμ 10 1 ≡ ∇ xμ1 ( R),μ 11 1 ≡ ∇ x ∇ yμ1 ( R), etc. We emphasize that in all of our numerical calculations there is no expansion in m and we include all contributions from the gradient expansion up to second order.
Using the notation defined above, the coefficients of the energy-momentum tensor at order τ 0 have the simple form At order τ 2 we have

C. Angular momentum
In this section we derive an expression for the angular momentum of the glasma per unit rapidity. Our method is similar to that of Ref. [8]. We define the tensor where R µ denotes a component of the position vector. The energy momentum tensor is divergenceless, and therefore the tensor in equation (21) satisfies the tensor equation ∇ µ M µνλ = 0.
Using Stokes' theorem one obtains a set of six conserved quantities where n µ is a unit vector perpendicular to the hypersurface Σ, γ is the induced metric on this hypersurface, and d 3 y is the corresponding volume element. The angular momentum is obtained from the Pauli-Lubanski vector where u γ is the vector that denotes the rest frame of the system. One can easily verify that equations (21,22,23) reduce to the usual definition of angular momentum in Minkowski space. We denote indices for spatial variables in Minkowski space with underscored latin letters, for example i ∈ (1, 2, 3) and x i is a component of the vector (x, y, z). We use n µ = (1, 0, 0, 0) so that Σ is a hypersurface of constant t, and with u γ = (1, 0, 0, 0) equation (23) becomes where d 3 x represents the spatial volume element in Minkowski space, and we have written the Poynting vector P i ≡ T 0i .
We work in Milne coordinates and use n µ = (1, 0, 0, 0) so that is defined on a hypersurface of constant τ . Using u γ = (1, 0, 0, 0) gives In our calculation angular momentum is not conserved, because the currents on the lightcone act as sources, and therefore we consider instead the angular momentum per unit rapidity. From equation (26) we obtain 1 We note that although the right side of equation (27) is independent of rapidity, our calculation is only meaningful close to mid-rapidity where boost invariance is a good approximation.
The integral over the transverse plane in equation (27) can be simplified using symmetry 1 We comment that our result is different from that of Ref. [8] where a slightly inconsistent procedure is used.
In that paper the authors define angular momentum in Minkowski space, on a surface of constant time, and then enforce separately that the integral should be calculated with τ held fixed.
considerations. The source charge distributions that we use (see Section II B) are even under the transformation R y → −R y . We will consider symmetric displacements of the ions relative to the collision axis ( b 1 = − b 2 ) so that the transformation R x → −R x interchanges the distributions for the first and second ions. Using these symmetries one can show that each component of the energy momentum tensor in Milne coordinates is either even or odd under the R x and R y parity transformations. We summarize these symmetries in Table I.
T 00 even even T 01 odd even T 02 odd even T 03 even odd T 11 even even T 12 even even T 13 odd odd T 22 even even T 23 odd odd T 33 even even Using the symmetry relation in Table I it is easy to see that the only non-zero component of the angular momentum per unit rapidity is

A. Notation and units
We remind the reader of the geometry of the collision we are considering. The two ions approach each other along the z-axis and collide at the origin, at time t = 0. Post collision, the first ion moves outward along the positive z-axis, and the second ion moves along the negative z-axis. We will consider collisions with non-zero impact parameter, which we denote GeV, Q s = 2 GeV and g = 1, unless stated otherwise. We consider lead-lead collisions, which corresponds to mass numbers A 1 = A 2 = 207, except for a few situations where we will explicitly specify different mass numbers. We will show below that our results to order τ 6 are valid to τ ≈ 0.05 fm, which corresponds toτ ≈ 0.51.

B. Physical observables and the gradient expansion
The form of the charge density function µ( x ⊥ ) that we use is discussed in section II B.
We use r 0 = 1.25 fm and a = 0.5 fm so that the radius of a nucleus with A = 207 is As explained in section II B, we allow for non-homogeneous nuclear density functions by performing a gradient expansion around the coordinate that gives the position of the center of each nucleus in the transverse plane. In figure 1 we show the density function with x ⊥ = (R x , 0), its first and second derivatives with respect to R x , and the quantity δ in equation (9) which must be small for the gradient expansion to converge.
The condition δ < 0.75 is satisfied in the region to the left of the vertical line in the figure.  We show that the former is almost completely insensitive to the gradient expansion, whereas the latter depends strongly on the behaviour of the nuclear density function µ( x ⊥ ) close to the nuclear radii. In section IV E we look at various projections of the Poynting vector which describes the flow of the energy of the gluon field. The leading-order contribution to the Poynting vector comes solely from the first order term in the gradient expansion, which means that we can restrict to the region where the expansion parameter δ is small, and still see clearly the contribution of the gradient terms. In section IV F we study the momentum anisotropy of the glasma by calculating the Fourier coefficients of the azimuthal distribution of the flow, and in section IV G we look at the spatial azimuthal asymmetry of the glasma by calculating the eccentricity. These calculations involve integration over the transverse plane, and are therefore potentially sensitive to the gradient expansion. One must show that results are largely insensitive to the choice of the integration limits, and we will find that this condition restricts us to the consideration of small impact parameters. This is because when the centers of the two ions are separated, the inner edge of the first/second ion, where the density changes rapidly, will be closer to the center of the second/first ion, where integrand can be large. In section IV H we look at the angular momentum of the glasma. In this case the gradient expansion severely limits the accuracy of the calculation, but we are able to see that the angular momentum carried by the glasma is many orders of magnitude smaller than the total angular momentum of the participant nucleons of the colliding nuclei [21,22].

C. Energy density
We look at the initial energy density E = T 00 mink at mid-spatial-rapidity (η = 0) for four different configurations of the colliding ions which are defined in Table II. The last row of the table shows the maximum initial energy density.  In figure 2 we show, for case B, the initial energy density, and the difference between the energy density atτ = 0.42 and the initial energy density. The energy density drops fastest at the centre and more slowly at the edges of the almond shaped interaction region.

Transverse and longitudinal pressures
We define the normalized longitudinal and transverse pressures as For a system in equilibrium p T /E = p L /E = 1/3. The glasma energy-momentum tensor at τ = 0 + has the diagonal form (with both indices raised) We remind the reader that in our notation the components of a position 4-vector in Minkowski coordinates are (t, z, x, y). Equation (30) shows that the initial longitudinal pressure is large and negative. The initial system is therefore not only far from equilibrium, but also far from the regime where a quasi-particle picture would be valid. As τ increases the longitudinal pressure grows and, because the energy-momentum tensor is traceless, the transverse pressure decreases.
In figure 3 we show the vector (p L /E, p T /E) in the transverse plane with b = η = 0. In the left panel we see −p L = p T at τ = 0. In the next two panels we use the biggest value of τ for which we trust the τ expansion at that order (these times are determined from additional results in this and the following sections, see for example figures 4 and 10). In the middle figure we include terms to order τ 4 and set τ = 0.04 fm. The figure shows that the vector has straightened slightly across the transverse plane. In the right figure we include all terms to order τ 6 and use τ = 0.052 fm. We see that the vector has straightened even more, but not uniformly. The authors of Ref. [23] suggested that the anisotropy of the transverse and longitudinal pressures should be characterized using the quantity which takes the value A T L = 6 at τ = 0 (using equation (30)) and would be zero in an  We also study how the behaviour of A T L depends on azimuthal angle (denoted φ), spatial rapidity, and impact parameter. As expected, A T L moves towards the equilibrium value more quickly when the impact parameter is smaller, and the region where the two ions overlap is greater. In figure 5 we show the quantity A T L at sixth order in the proper time expansion as a function ofτ , for different values of η and φ. We consider φ = 0, which corresponds to R in the reaction plane, and φ = π/2, where R is perpendicular to the reaction plane. The graph shows that A T L drops more quickly when either the azimuthal angle or the spatial rapidity increases. In figure 6 we show contour plots of A T L in the transverse plane for η = 0, b = 0 and τ = 0.045 fm at order τ 4 and τ 6 . One sees that when the order τ 6 terms are included, the region of the transverse plane where A T L is small is significantly broader. In figure 7 we show contour plots of A T L in the transverse plane for η = 0, b = 0 at order τ 6 for three different times. One sees that the value of A T L shrinks across the transverse plane as τ increases. We can also use the quantity A T L to demonstrate that our results are not strongly dependent on the UV and IR scales that enter the calculation (Q s and m in our notation). This is important because the exact values of these scales are not known, and also because the way that these scales enter the calculation depends on the method chosen to perform the

Transverse pressure anisotropy
In this section we look at a quantity that can be used to characterize the asymmetry of the transverse pressure. We define [25] {A xy } ≡ T yy − T xx T xx + T yy (32) where the angular brackets indicate integration over the transverse plane. For comparison we will also calculate The leading-order contribution to {A xy } comes from the first order term in the gradient expansion and therefore this quantity, in contrast to {A T L }, will be sensitive to the region of the transverse plane that is close to the edges of the nuclei. We must verify that the integral is largely independent of the upper limit that is used to perform the two dimensional integration over the transverse plane, which we call R max . In figure 9 we

E. Radial flow
To describe the radial flow of the expanding glasma in the transverse plane we look at the radial projection of the transverse Poynting vector P ≡R i T i0 . In the left panel of figure   10 we show this quantity at η = 0, b = 6 fm, R = 5 fm and φ = π/2 at different orders in the τ expansion. At lowest order P increases linearly with time. Including higher order contributions we see that P increases more slowly with time as the system expands. If we keep only terms at order τ 3 , it appears that P actually starts to decrease whenτ 0.4, however the near field expansion is not valid for these times when only up to cubic terms are included. The result at order τ 5 shows a less pronounced flattening up to aboutτ = 0.5, which again indicates the limit of validity of the near field expansion. In the right panel of figure 10 we show P at η = 0, b = 6 fm and R = 5 fm for φ = 0, which corresponds to R in the reaction plane, and φ = π/2, where R is perpendicular to the reaction plane. One sees that the flattening is more pronounced when the azimuthal angle is smaller.  In figure 12 we show the Poynting vector at R y = 0 in the (η,R x ) plane, in arbitrary units.
the right moving ion that has been displaced in the positive x direction, the sign of the Poynting vector is predominantly negative. This corresponds to the negative value of the directed flow coefficient discussed in section IV F. The second panel shows that the system expands more strongly in the wake of the larger nucleus.  plane. In our calculation the reaction plane is known, and fluctuations are not included. One consequence is that the Fourier coefficients in our calculation exhibit specific symmetries: the coefficients v n with n odd are rapidity odd, and those with n even are rapidity even.
To estimate the error in the Fourier coefficients, we perform the integrals using, for each value of b, a set of 15 evenly spaced values of R max between 4.5 fm and 7.0 fm. We average the results and calculate the error bar from the standard deviation. In figure 14 we show the first three Fourier coefficients at η = 0.5 and τ = 0.04 fm. As in section IV D 2, we find that the calculation is reasonably insensitive to the upper limit of the integration for impact parameters b 2.5 fm.
In figure 15 we look at the first three Fourier coefficients as functions of rapidity and impact parameter, at τ = 0.04 fm. We use R max = 5.9 fm in all calculations. In the left panel of  This flattening of all three curves inside the region where the τ expansion converges provides some evidence that the Fourier coefficients will not change rapidly immediately outside the region of validity of the near field expansion. We also comment that the radial flow shows similar behaviour (see figure 10).  Our results for the second and third Fourier coefficients are of the same order as experimental values [26,27], although our result for |v 1 | is much bigger than expected [28]. However, it is usually assumed that anisotropy develops mostly during the hydrodynamic evolution that follows the glasma phase, and in this context our results are surprisingly large for all three Fourier coefficients. Azimuthal correlations in glasma have been investigated by several other groups, using different approaches and looking at different regimes, which makes comparison difficult. In Ref. [25] the authors use a different method and consider impact parameters that are larger than the maximum value for which the gradient expansion can be trusted, and times that are beyond the validity of the near field expansion, but still obtain smaller values of v 2 . In Ref. [29,30] the origins of azimuthal correlations in the CGC approach are studied, but quantitative results for a glasma system are not obtained.

G. Eccentricity
Spatial deviations from azimuthal symmetry can be characterized with the quantity [31,32] where E( R) denotes the energy density.  In figure 18 we show v 2 , ε 2 and v 2 /ε 2 at τ = 0.04 fm and η = 0 using R max = 5.9 fm.
In figure 19 we show the same curves normalized so that the value at b = 0.5 fm is set to one. The results show a clear correlation between the spatial asymmetry introduced by the initial geometry, and the anisotropy of the azimuthal distribution of the gluon momentum. Correlations of this kind are considered characteristic of the onset of hydrodynamic behaviour.

H. Angular Momentum
The angular momentum of the glasma can be calculated from equation (28). In the left panel of figure 20 we show our results at three different times, using R max = 5.9 fm. The error bars in the right panel are obtained by calculating the standard deviation of the results using a set of evenly spaced values for R max between 4.5 fm and 7.0 fm. The error bars that are produced by this procedure are large, even for small impact parameters (for comparison see figures 14 and 17). This happens because the problem discussed in section IV B is much more serious in the calculation of angular momentum than it was in sections IV F and IV G.
The dominant contribution to the angular momentum comes from the parts of the nuclei that are farthest from the collision centre, with respect to which angular momentum is calculated, but these are the regions where the gradient expansion is least to be trusted. The angular momentum is negative, which is expected when the ions moving in the positive/negative z directions are displaced in the positive/negative x directions. We note that in spite of the size of the error bars in figure 20, the general shape of the curves is consistently reproduced when the value of R max is changed. Furthermore, this shape matches the basic form of the results in Refs. [21,22]. It is especially interesting that in all three calculations, the peak occurs at ≈ 2.0 − 2.5 fm, which is within the range of impact parameters where the gradient expansion that we use can be trusted. However, it is important to note that the values for angular momentum obtained in [21,22] are ≈ 10 5 at RHIC energies, and even larger at LHC energies, and are thus five to six orders of magnitude larger than our results.
In figure 21 we show the angular momentum as a function of proper time, at fixed impact parameter b = 2 fm at different orders in the τ expansion. We have used the value R max = 5.9 fm in all cases. The figure shows that the τ expansion appears to converge for much larger times than the other quantities we have calculated. The authors of Ref. [8] observed that the near field expansion method can produce very large values of angular momentum, if one considers large enough times. We have verified that our calculation reproduces this behaviour, but the sign of the angular momentum changes, and the numerical value of the result depends strongly on the order of the τ expansion. Both of these properties indicate that there is no reason to believe the result is physical. Our results indicate that the glasma carries only a very small imprint of the primordial angular momentum, which means that the majority of the angular momentum is carried by valence quarks, for times at which the calculation is valid. This result casts doubt on the idea of a rapidly rotating initial glasma state that could be observed via polarization of final state hadrons.

V. CONCLUSIONS
We have used a CGC approach and an expansion in proper time to derive an analytic result for the glasma energy-momentum tensor to sixth order in the proper time. We have taken into account some aspects of nuclear structure using spatially dependent nuclear density functions. In our previous paper [1] we gave a detailed description of the steps involved in our calculation. In this paper, we have concentrated on the physical results that can be obtained from our final expression for the energy-momentum tensor.
For most of the quantities that we calculated, the proper time expansion can be trusted to about τ = 0.05 fm, and we have shown that the glasma moves towards equilibrium until the time at which the near field expansion breaks down. Using simple arguments based on the uncertainty principle [1], one can argue that this upper bound of the region where the near field expansion converges reaches far beyond the lower bound at which we no longer trust the classical description that is inherent in our method. Our calculation also requires an expansion in gradients of the nuclear density function. Some of the quantities we have calculated are almost entirely insensitive to this expansion, and in other cases the gradient expansion restricts us to the consideration of small impact parameters.
We have calculated the first three Fourier coefficients of the azimuthal distribution of the momentum field of the glasma. Our results are larger than is generally expected, which is interesting because it contradicts the usual assumption that azimuthal anisotropy is mostly generated during the hydrodynamic evolution of the plasma. We have also calculated the eccentricity of the glasma, which describes the spatial azimuthal asymmetry of the system.
We have found a sizable correlation between the elliptic flow coefficient and the eccentricity, which indicates that the spatial asymmetry introduced by the initial geometry is effectively transmitted to the azimuthal distribution of the gluon momentum field. This result is interesting because a correlation of this kind is an indication of the onset of hydrodynamics. We have formulated a calculation of the angular momentum per unit rapidity, on a hypersurface of constant proper time. Our results are much smaller than the total angular momentum of the participating nucleons [21,22,28], which shows that most of the angular momentum carried by valence quarks is not transmitted to the glasma. These results are significant because they contradict the picture of a rapidly rotating initial glasma state that has been proposed by several authors [9,10] and led to experimental searches for a polarization effect in the hyperons and vector mesons produced in heavy ion collisions, which have been, so far, largely unsuccessful [11,12].
Finally, we comment that all of the calculation done in this paper make use of analytic solutions of the near field expanded YM equation which we have obtained to order τ 6 , and these results are also useful in other contexts. We are currently using these solutions to perform a calculation of the transport properties of heavy quarks in glasma. conventional definitions We define the relative and average transverse coordinates We will write unit vectors asr = r/| r| = r/r andR = R/| R| = R/R and use standard notation for derivatives, for example and similarly for y ⊥ , r and R. In this section we define the azimuthal distribution of the flow vector T i0 ( x ⊥ ). The azimuthal angle is measured with respect to the x-axis and is written We define the distribution where we have introduced the weighting function and the normalization factor The distribution P (φ) can be decomposed into Fourier harmonics as where coefficients v n are given by the relation v n = 2π 0 dφ cos(nφ) P (φ) . (C6)