1. 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 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 , we also use Milne coordinates , where is the proper time and is the spacetime rapidity. Finally, we choose units such that , where c is the speed of light in a vacuum, is the Boltzmann constant, and ℏ is the Planck constant divided by .
4. Numerical Implementation
4.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 .
The discretization of the transverse variables is straightforward. The discretized version of Equation (
11) is
with
Since
D is a rank 4 tensor with
components, we obtain
eigenfunctions
, and
L eigenvalues
, with
. 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
in every factor
and in the Hankel functions. For the complex exponential, we use
. 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).
4.2. Boundary Conditions
We use periodic boundary conditions, which means that the indices
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
4.3. 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
and
. We use only the second because it has positive frequency behavior at large times
From now on, we suppress the superscript on the Hankel function. When , the Hankel function oscillates like and the derivative diverges. Numerically, we must start the evolution at a small positive time, which we choose as . 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
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
so that
and Max
. If this
can be found, use Equation (
33) with
. If
does not exist, then increase the chosen value of
and try again. Using this procedure, we can find
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
.
4.4. 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 . 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: .
4.5. Initial Conditions
The initial conditions that we use for the background field and its derivative are
The argument of the sine function is at and , and zero at , so the field has negative on the left side of the lattice and positive 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.
5. Results and Discussion
All of our results are obtained with
spatial grid points,
points for the rapidity coordinate, and
configurations. The initial conditions for the background field are obtained from (
35) with
,
and
.
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
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
. The right panel shows the two pressures normalized by the energy density, both approaching 1/3, again up to
. 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
versus
with
, 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
. 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
is the fitted function that gives the green curve. The result is
.