Isotropization of a Rotating and Longitudinally Expanding ϕ4 Scalar System

We study numerically the evolution of an expanding system of scalar fields. The initial configuration is non-isotropic and rotating. We calculate the energy–momentum tensor and angular momentum vector of the system. We compare the time scales associated with the isotropization of the transverse and longitudinal pressures, and the decay of the initial angular momentum. We show that even a fairly large initial angular momentum decays significantly faster than the pressure anisotropy.


Introduction
In this paper, we study the time evolution of an expanding system of rotating massless real scalar fields with quartic coupling. Our calculation is based on the method developed in [1,2]. Observables calculated in a loop expansion exhibit divergences at next-to-leading order, which originate from instabilities in the classical solutions. The effect is seen in a calculation of the energy-momentum tensor at the next-to-leading order, where the energy density and pressures of the system diverge rapidly with increasing time. Gelis and his collaborators have shown that this problem can be cured using a resummation scheme that collects the leading secular terms at each order of an expansion in the coupling constant. This resummation can be performed by allowing the initial condition for the classical field to fluctuate, and averaging over these fluctuations. They have shown that a system of scalar fields isotropizes when this resummation is performed [2].
The motivation behind the development of this approach is to study the thermalization of the glasma phase of the matter created in a relativistic heavy ion collision. It is known that a hydrodynamic description, which is valid when the system is fairly close to thermal equilibrium, works well at very early times (∼1 fm/c). Approaches that are based on kinetic theory descriptions of the scattering of quasi-particles cannot explain this rapid thermalization. Another possibility that has been studied extensively is that the system is strongly coupled, even at very high energies. The proposal of Gelis et al. is that rapid thermalization could be achieved by a resummation of quantum fluctuations. The color glass condensate (CGC) effective theory provides a natural framework for this formulation [3][4][5]. At very early times, the system is best described as a system of strong classical fields that can be obtained from solutions of the Yang-Mills equation using a CGC approach. The spectrum of quantum fluctuations was derived in [6]. The success of the resummation method was demonstrated in [7], where the authors showed that pressure isotropiztion occurs in an SU(2) analog of QCD.
Our ultimate goal is to use the Gelis et al. approach to study the creation and evolution of angular momentum in a glasma. This is interesting in the context of recent proposals that the glasma is produced in a rapidly rotating state, which could be detected by looking for the polarization of produced hyperons. There have been calculations that predict very large values for the initial angular momentum of the system [8][9][10], but significant final state polarization effects have not been observed [11,12]. In this paper, we develop a formulation to calculate the angular momentum of a system of real scalar fields. We present preliminary results that indicate the angular momentum relaxes to a small value on a time scale significantly smaller than the time scale for pressure isotropization. If a similar result is obtained in a QCD glasma, it would be consistent with the observations in [11,12]. We also comment that a calculation of angular momentum in glasma was conducted in [13], using a CGC approach with a proper time expansion, and it was found also that large amounts of angular momentum were not produced.
Since computations in a gauge theory are considerably more complicated, we will work with a scalar theory. While it is true that QCD and scalar φ 4 theory are different in many ways, they have important similarities in the context of this calculation because they both have unstable modes and are scale invariant at the classical level. In addition, we will mimic the kinematics of a relativistic nuclear collision by working in Milne coordinates with a rapidity independent background field. Milne coordinates are suitable because in a nuclear collision, there is a preferred spatial direction provided by the collision axis, and in the high energy limit, one expects invariance under Lorentz boosts in the z-direction.
This paper is organized as follows. In Section 2, we describe the method, and in Section 3, we formulate the calculation of the energy-momentum tensor and angular momentum. Some details of our numerical procedure are discussed in Section 4. In Section 5, we present our results, and in Section 6 we make some concluding remarks.
Throughout this paper, the spacetime is always taken to be Minkowski, with the signature (+, −, −, −). In addition to standard inertial coordinates (t, x, y, z), we also use Milne coordinates (τ, x, y, η), where τ is the proper time and η is the spacetime rapidity. Finally, we choose units such that c = k B =h = 1, where c is the speed of light in a vacuum, k B is the Boltzmann constant, andh is the Planck constant divided by 2π.

Preliminaries
We consider a massless self-interacting real scalar field φ with quartic coupling. The Lagrangian density is given by where g is the coupling constant. To mimic the kinematics of a high energy nuclear collision, we work in Milne coordinates (τ, η, x ⊥ ) with Under a Lorentz boost in the z-direction, the proper time is unchanged and η is shifted by a constant. The metric in Milne coordinates is (2) Figure 1 shows the curves of constant τ and η.

The Resummation Procedure
As explained in [1,2], observables calculated in a loop expansion exhibit secular divergences at next-to-leading order that originate from instabilities of the classical solutions. Gelis et al. proposed to cure this problem using a resummation scheme that collects the leading secular terms at each order of an expansion in the coupling constant, by averaging over an ensemble of initial conditions. The energy-momentum tensor is ultraviolet divergent, but the divergence corresponds to a vacuum contribution and can be removed by repeating the calculation with the background field set to zero, and subtracting the results. This vacuum subtraction was performed for all the calculations presented in this paper.
The equation of motion for the scalar field obtained from the Lagrangian (1) is where the "dot" indicates a derivative with respect to τ, and ∆ ⊥ is the transverse Laplacian operator. The initial field is written as the sum of a background field contribution, ϕ, which is assumed to be boost invariant and therefore independent of η, and an η-dependent fluctuation, which we call α The initial time τ 0 is chosen to be small but nonzero (see Section 4.3 for further discussion). The initial background field ϕ(τ 0 , x ⊥ ) is discussed in Section 4.5. The index χ in Equation (4) indicates that we have a Gaussian ensemble of initial conditions defined as The index K labels the momentum variables (ν, k ⊥ ) that are conjugate to the coordinatespace variables (η, x ⊥ ), respectively. The notation c χ K indicates an element in a Gaussiandistributed ensemble of N χ random numbers, with variance We use the momentum space integration measure and the delta function in Equation (6) is defined so that dKδ KL = 1. The mode functions a K ≡ a ν k ⊥ (τ 0 , η, x ⊥ ) are obtained from the linearized equations of motion and normalized so that dK(a K , a L ) = 1 with Separating variables and performing the normalization, one finds where the χ k ⊥ is the solution of the eigenvalue equation The field φ χ (τ, η, x ⊥ ) at finite proper time is obtained by solving Equation (3) with the initial condition φ χ (τ 0 , η, x ⊥ ) obtained from Equations (4), (5)-(7), (10) and (11). From this point on, we drop the subscript χ.

Energy Momentum Tensor
The energy-momentum tensor of theory (1) is The invariance of the Lagrangian under the conformal transformation implies that T µν is traceless on shell.
The expressions for the energy and pressure are where V(φ) = g 2 φ 4 /4!. In terms of the energy and pressure, the trace condition is

Angular Momentum
We use the standard Pauli-Lubanski formalism [14,15] to obtain an expression for the angular momentum in terms of the energy-momentum tensor. We define the tensor field where R µ is the coordinate vector. Using Stokes' theorem, one obtains a set of six conserved quantities where n µ is a unit vector perpendicular to the hypersurface Σ, γ ij 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. Equations (16)-(18) will give where the energy-momentum tensor is given in Equation (12).
To find the angular momentum on a surface of constant τ, we define In Minkowski coordinates, this gives n µ = (cosh(η), 0, 0, −sinh(η)), and it is easy to verify that n Milne µ = (1, 0, 0, 0), as expected. The fluid velocity is the local rest frame in comoving coordinates, which is written u ρ Milne = (1, 0, 0, 0). In Minkowski coordinates, this becomes u ρ = (cosh(η), 0, 0, sinh(η)). We could calculate the angular momentum directly in Minkowski coordinates, or alternatively, we could perform the calculation in Milne coordinates and perform a coordinate transformation to obtain the Minkowski space result. We checked our computations by verifying that both calculations give the same result. The components of the angular momenta about each of the Minkowski coordinate axes are We note that all components of the angular momentum are dimensionless (in natural units, withh = 1).

Lattice Discretization
We discretize in both directions in the transverse plane with L grid points and lattice spacing set to 1, which effectively means that we define all dimensionful quantities in terms of the transverse lattice grid spacing. The rapidity variable η is discretized with N grid points and lattice spacing h. We consider a unit slice of rapidity, and therefore take h = 1/N.
The discretization of the transverse variables is straightforward. The discretized version of Equation (11) is Since D is a rank 4 tensor with L 4 components, we obtain L 2 eigenfunctions χ e ij , and L eigenvalues (λ 2 ) e , with e ∈ (1, L 2 ). The normalized eigenfunctions are and the momentum integration is discretized as Since the spatial lattice spacing is set to 1, an integral over transverse coordinates is discretized as The discretization of the longitudinal variables is a little more subtle. The constraint gives and we replace ν → ε v in every factor e πν/2 and in the Hankel functions. For the complex exponential, we use e iνη → e 2πivn N . The integral over ν becomes a sum over v using Combining these expressions, we find the discretized versions of Equations (4), (6) and (10): To verify that discretization is performed correctly, we checked the discretized version of the normalization condition (9).

Boundary Conditions
We use periodic boundary conditions, which means that the indices (i, j) that correspond to the transverse spatial coordinates are defined as modulo L, and the index n for the rapidity is modulo N. The boundary conditions satisfy the self-adjointness condition

Hankel Functions
The differential equation for the mode function was solved by separating variables, which gives the solution in (10). The time-dependent part of the equation is of the second order, and has two independent solutions, which are the Hankel functions H (1) iν (λτ) and H (2) iν (λτ). We use only the second because it has positive frequency behavior at large times From now on, we suppress the superscript (2) on the Hankel function. When τ → 0, the Hankel function oscillates like e ±iτν and the derivative diverges. Numerically, we must start the evolution at a small positive time, which we choose as τ 0 = 10 −2 . One can check that the value chosen for this small initial time does not change the results at finite times.
We describe below our method to calculate the Hankel functions. First, we define the scaled function h iν (λτ) = e πν/2 H iν (λτ) (32) which is easier to calculate numerically. At large times, one can obtain the scaled Hankel function for given values of ν and λ from the asymptotic series This expression must be used carefully because the series does not converge for arbitrarily large values n. We proceed as follows. For a given value of ν and λ, choose some value of τ and look for a value of k max so that t k max +1 < 10 −9 and Max(t k≤k max ) < 10 6 . If this k max can be found, use Equation (33) with n = k max . If k max does not exist, then increase the chosen value of τ and try again. Using this procedure, we can find h iν (λτ) and its first derivative for each value of ν and λ, for some (possibly very large) time. We then use adaptive fifth-order Runge-Kutta to find each Hankel function at the initial time τ 0 .

Discretized Derivatives
The conservation equation is an exact equation that should be satisfied whether or not the system is in equilibrium. Additionally, we should have that the trace of the energy-momentum tensor is zero, so that Equation (15) is satisfied. It is easy to show analytically that these conditions are satisfied for background fields if we use forward derivatives: The point is that while centered derivatives are not wrong, much larger lattices must be used to achieve the same numerical accuracy. For angular momentum, the situation is different. All contributions to the angular momentum have an integral of the form dxφ ∂ x φ. If the initial value ofφ is constant, the integrand is a total derivative and therefore the integral will give zero. However, this is not well satisfied numerically with forward derivatives. In the calculation of angular momentum, it is therefore better to use centered derivatives:

Initial Conditions
The initial conditions that we use for the background field and its derivative are The argument of the sine function is ∓π/2 at i = 1 and i = L, and zero at i = (L + 1)/2, so the field has negativeφ 0 on the left side of the lattice and positivė ϕ 0 on the right side.
The astute reader will note that our initial classical field is not periodic, and therefore does not respect our boundary conditions. The reason is that we wish to avoid problems that may arise when resonant modes are considered, which in the present model would correspond to the normal modes of the finite spatial lattice. For a large enough lattice, all modes are effectively periodic, and it is therefore expected that the precise form of the initialization is not important.

Results and Discussion
All of our results are obtained with L = 41 spatial grid points, N = 120 points for the rapidity coordinate, and N χ = 256 configurations. The initial conditions for the background field are obtained from (35) with ϕ 0 = 15, k x = k y = 1/ √ 2 andφ 0 = 10. To investigate if the system obeys Equation (15), we compare the energy density and the sum of the pressures. This is shown in Figure 2. One sees that after some initial oscillations have damped out, the condition = 2p T + p L is well satisfied. To see if the system approaches an isotropic state, and if it obeys an equation of state, we look at the transverse and longitudinal pressures. The left panel of Figure 3 shows that, after some initial oscillations have disappeared, the transverse and longitudinal pressures approach each other up to a time of about τ ≈ 160. The right panel shows the two pressures normalized by the energy density, both approaching 1/3, again up to τ ≈ 160. For large times, the simulation breaks down, which is not unexpected when one studies the dynamics of an expanding system inside a box of finite size. In Figure 4, we show the three components of the angular momentum in Equation (21). The z component, which depends weakly on the rapidity, is averaged over the unit slice of rapidity that we consider. In comparison with the energy and pressure, the oscillatory behavior is more severe and does not completely disappear. To obtain a better idea of the overall behavior, we also plot the accumulated average for each component, which is shown in Figure 4 with the thick lines. In each case, the darker color corresponds to the average of the component with the same but lighter color. The figure shows that even a fairly large initial angular momentum decays very quickly. We want to compare the time scales for the isotropization of the pressures, and the decay of the initial angular momentum. In Figure 5, we show in blue the curve in the left panel of Figure 3 over the range of τ for which the decay is strongest. To produce the light green points, we took the data for | L| versus τ with τ > 12.0, where the large initial fluctuations are mostly gone, and shifted the first point (which was (12.0, 9.10)) so that it sits on top of the first point of the data that made the blue curve. The dark green line is a fit obtained for these data using the function A + B/τ + Ce −Dτ . The plot shows clearly that the initial angular momentum decays much more quickly than the pressure anisotropy, although the dispersion in the data is large. To quantify this dispersion, we calculated where the sum is over the points shown in green in Figure 5 and L fit (τ) is the fitted function that gives the green curve. The result is σ = 0.22.

Conclusions
In this paper, we presented some preliminary results from our study of the angular momentum in an expanding system of rotating massless scalar fields. Our results indicate that even when a large amount of angular momentum is put into the system, it decays very rapidly. Future work will include an investigation of how much these results depend on the exact form of the initialization and the boundary conditions, and possibly the extension of the calculation to physical theories, such as QCD.