# Evidence for universal flow and characteristics of early time thermalization
in a scalar field model for heavy ion collisions

Margaret E. Carrington  Department of Physics, Brandon University, Brandon,
Manitoba R7A 6A9, Canada  Winnipeg Institute for Theoretical Physics,
Winnipeg, Manitoba, Canada  Wade N. Cowie  Department of Physics, University
of Manitoba, Winnipeg, Manitoba, R3T 2N2, Canada  Winnipeg Institute for
Theoretical Physics, Winnipeg, Manitoba, Canada  Gabor Kunstatter  Department
of Physics, University of Winnipeg, Winnipeg, Manitoba, R3M 2E9 Canada
Department of Physics, Simon Fraser University, Burnaby, British Columbia, V5A
1S6 Canada  Winnipeg Institute for Theoretical Physics, Winnipeg, Manitoba,
Canada  Christopher D. Phillips  Department of Physics, Brandon University,
Brandon, Manitoba R7A 6A9, Canada

###### Abstract

We study numerically the evolution of an expanding strongly self-coupled real
scalar field. We use a conformally invariant action that gives a traceless
energy-momentum tensor and is better suited to model the early time behaviour
of a system such as QCD, whose action is also conformally invariant. We
consider asymmetric initial conditions and observe that when the system is
initialized with non-zero spatial eccentricity, the eccentricity decreases and
the elliptic flow coefficient increases. We look at a measure of transverse
pressure asymmetry that has been shown to behave similarly to the elliptic
flow coefficient in hydrodynamic systems and show that in our system their
behaviour is strikingly similar. We show that the derivative of the transverse
velocity is proportional to the gradient of the energy in Milne coordinates
and argue that this result means that transverse velocity initially develops
in the same way that it does in hydrodynamic systems.

We conclude that some aspects of the early onset of hydrodynamic behaviour
that has been observed in quark-gluon plasmas are seen in our numerical
simulation of strongly coupled scalar fields.

##  I Introduction

The quark gluon plasma (QGP) has been observed experimentally to be well
described by hydrodynamics at very early times ($1$ fm/c $\sim 3\times
10^{-24}$ seconds) after formation. This implies that some degree of
thermalization must happen on even shorter time scales. Such rapid
thermalization does not seem possible within a standard kinematical
description of the QGP. A field theoretic approach therefore appears to be
required, something that is extremely difficult in full Quantum Chromodynamics
(QCD).

Progress has recently been made in Dusling:2010rm ; Epelbaum:2011pc ;
Dusling:2012ig  (DEGV) in the context of a simple model that has some
properties in common with QCD but is nonetheless considerably more tractable.
DEGV considered a real, self-interacting scalar field in four spacetime
dimensions. Like QCD this theory is scale invariant and exhibits secular
divergences in its perturbative solutions. These secular divergences are due
to instabilities in the classical solutions. DEGV used what is called the
classical statistical approximation to resum the fastest growing divergences
in this theory at each order in perturbation theory and showed that the scalar
field system isotropizes when this resummation is performed Dusling:2012ig .
The only observables considered were the energy and longitudinal/transverse
pressures of the system, and it was found that the difference between the
resummed longitudinal and transverse pressures tends towards zero. In a
previous paper Carrington:2022ccm  we used the the classical statistical
approximation to study angular momentum in a system of expanding scalar
fields. Some issues with the classical statistical approximation are discussed
in Berges:2013lsa ; Epelbaum:2014yja . A completely different technique to
resum the dominant contributions based on a 2-particle irreducible effective
action approach can be found in Gelis:2024iar .

In this paper we study the extent to which an expanding system of scalar
fields behaves hydrodynamically, and how hydrodynamic behaviour develops with
time. We start with a somewhat different theory than DEGV, one that we claim
more closely models gluon behaviour without losing the inherent simplicity of
a real scalar field. Specifically, we consider a real scalar that is non-
minimally coupled to the background geometry so as to produce a conformally
invariant action, i.e. one that is invariant up to a total derivative under
local scale transformations as well as global ones. This action produces an
classical energy momentum tensor (EMT) that is manifestly symmetric, conserved
and traceless on-shell without altering the equations of motion on a Minkowski
background. The theory therefore shares with QCD its conformal symmetry and
the resulting tracelessness of the EMT. In the rest of this paper we refer to
the EMT obtained from the non-minimally coupled conformally invariant action
as the conformal EMT. The conformal and canonical EMTs are almost
indistinguishable at late times111Time is measured in units of the transverse
coordinate lattice spacing. but at early times some components of the
canonical EMT are characterized by rapid fluctuations that are absent in the
conformal EMT.

We use the conformal EMT to study some observables that are related to the
transverse asymmetry of the system of fields. Transverse asymmetry has been
studied for many years in the context of quark-gluon plasmas. In a
relativistic non-central collision there is an initial spatial asymmetry that
rapidly decreases with time, which means anisotropic transverse momentum can
develop only at very early times. Since anisotropic flow is sensitive to the
system’s properties very early in its evolution, it can provide direct
information about the early stages of the system. We find that if our system
of scalar fields is initialized with spatial asymmetry, the asymmetry is
transmitted into the momentum field. This behaviour mimics what is seen in
heavy ion collisions and is commonly considered an indication of the onset of
some kind of hydrodynamic behaviour. We also look at a measure of transverse
pressure asymmetry that has been shown to behave similarly to the elliptic
flow coefficient in quark gluon plasmas. We show that in our system the
behaviour of these two quantities is qualitatively similar. Finally we
investigate a proposal Vredevoogd:2008id  that under very general conditions
the derivative of the transverse velocity will be proportional to the gradient
of the energy, and discuss the significance of the result.

This paper is organized as follows. In section II we discuss the classical
theory and the definition of the energy-momentum tensor (EMT). In section III
we describe the resummation method and give some details of the numerical
procedure. We explain how the discretization of the equations is done, and our
choice of boundary conditions and initial conditions. In section IV we discuss
how the energy density and the transverse and longitudinal pressures are
obtained from the EMT. We derive some equations that characterize the
hydrodynamic behaviour of the system, and we define functions that measure the
eccentricity, elliptic flow and transverse pressure asymmetry of the system.
In section V we present our results and section VI concludes with a summary
and some observations. Some useful equations are collected in Appendix A.

Throughout this paper the spacetime is taken to be Riemannian with signature
$(+,-,-,-)$. In addition to Minkowski coordinates $(t,z,x,y)$ we will also use
Milne coordinates $(\tau,\eta,x,y)$, where $\tau$ is the Bjorken time and
$\eta$ is spacetime rapidity. We choose units such that $c=k_{B}=\hbar=1$,
where $c$ is the speed of light in vacuum, $k_{B}$ is the Boltzmann constant,
and $\hbar$ is the Planck constant divided by $2\pi$.

##  II The classical theory

###  II.1 Canonical EMT

The Lagrangian density used by DEGV is that of a massless real scalar
field222We write everything covariantly since this will be useful when we
switch to Milne coordinates, which are related to Minkowski coordinates by a
non-linear coordinate transformation:

$\displaystyle{\cal
L}=\frac{1}{2}\nabla^{\mu}\varphi\nabla_{\mu}\varphi-V\text{~{}~{}with~{}~{}}V=\frac{g^{2}}{4!}\varphi^{4}\,.$
(1)

The equation of motion is

$\displaystyle\nabla_{\mu}\nabla^{\mu}\varphi+\frac{g^{2}}{3!}\varphi^{3}=0~{}~{}\to~{}~{}\varphi\Box\varphi+4V=0\,.$
(2)

The canonical EMT is

$\displaystyle T_{\rm can}^{\mu\nu}=\frac{\delta{\cal
L}}{\delta(\nabla_{\mu}\varphi)}\nabla^{\nu}\varphi-g^{\mu\nu}{\cal
L}=\nabla^{\mu}\phi\nabla^{\nu}\phi-g^{\mu\nu}\left[\frac{1}{2}\nabla^{\alpha}\phi\nabla_{\alpha}\phi-\frac{g^{2}}{4!}\phi^{4}\right]\,.$
(3)

Its divergence is

$\displaystyle\nabla_{\mu}T_{\rm can}^{\mu\nu}$ $\displaystyle=$
$\displaystyle(\Box\varphi)(\nabla^{\nu}\varphi)+(\nabla_{\mu})(\nabla^{\mu}\nabla^{\nu}\varphi)-g^{\mu\nu}(\nabla_{\alpha}\varphi)(\nabla_{\mu}\nabla^{\alpha}\varphi)+g^{\mu\nu}\frac{g^{2}}{3!}\varphi^{3}\nabla_{\mu}\varphi$
(4) $\displaystyle=$
$\displaystyle(\Box\varphi)(\nabla^{\nu}\varphi)+\frac{g^{2}}{3!}\varphi^{3}\nabla^{\nu}\varphi=0$

where we used the equation of motion in the last line. The conservation law
(4) implies the conservation on shell of the momentum four-vector defined by

$\displaystyle{{\mathcal{P}}}^{\mu}=\int d^{3}x\sqrt{-g}T^{\mu 0}$ (5)

provided the metric does not explicitly depend on the spacetime
coordinates333This condition is not satisfied in Milne coordinates. and the
fields vanish on the boundaries.

The trace of the canonical EMT is

$\displaystyle T^{\mu}_{{\rm can}~{}\mu}$ $\displaystyle=$
$\displaystyle-\nabla_{\mu}\varphi\nabla^{\mu}\varphi+4V$ (6) $\displaystyle=$
$\displaystyle-\nabla_{\mu}\varphi\nabla^{\mu}\varphi-\varphi\Box\varphi=-\nabla_{\mu}(\varphi\nabla^{\mu}\varphi)$

using the equation of motion to go from the second line to the third. The
canonical EMT is therefore not traceless. The QCD canonical EMT is also not
traceless, and additionally it is not symmetric or gauge invariant. These
problems can be remedied using several different methods (see Blaschke:2016ohs
for a detailed discussion). In this work we define the EMT from a conformally
invariant scalar field action. It is manifestly traceless due to the conformal
invariance of the action. It is also conserved and symmetric.

###  II.2 Covariant EMT

In this section we derive the conformal EMT for the real massless scalar field
that, like the EMT used for QCD, is traceless. In section IV.1 we explain
further why the definition given in this section is a better choice for our
purposes.

The real scalar field Lagrangian (1) is invariant under constant scale
transformations in which the spacetime metric and scalar field transform as:

$\displaystyle g_{\mu\nu}\to\Omega
g_{\mu\nu}\text{~{}~{}and~{}~{}}\varphi\to\Omega^{-\tfrac{1}{2}}\varphi$ (7)

with constant $\Omega$. It is however not invariant under local scale
transformations (i.e. conformal transformations) for which $\Omega=\Omega(x)$
is an arbitrary function of the spacetime coordinates that goes to unity on
the boundaries. The Lagrangian density of QCD is conformally invariant and to
incorporate as many of the fundamental properties of QCD into our theory as
possible, we start with a conformally invariant theory. The Lagrangian density
is PT :

$\displaystyle{\cal
L}_{\xi}=\frac{1}{2}\sqrt{-g}(g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}\varphi-\xi
R\varphi^{2}-V)$ (8)

where $R$ is the scalar curvature of the spacetime. The choice $\xi=0$ gives
the original Lagrangian (1) and is called the minimally coupled theory. In
four dimensions with $\xi=1/6$ the Lagrange density changes by a total
derivative under the transformation (7) with non-constant $\Omega(x)$.

The equation of motion for the scalar field is

$\displaystyle\left(\Box+\frac{R}{6}+\frac{g^{2}}{3!}\varphi^{2}\right)\varphi=0\,.$
(9)

The metric is treated as a fixed background field and is therefore not
dynamical. When evaluated on a Minkowski background (zero curvature), eq. (9)
reduces to the equation of motion (1) of the original theory.

The requirement that the action be invariant under a general coordinate
transformation $x^{\mu}\to x^{\prime\mu}-\epsilon^{\mu}(x)$ leads to the
definition of a conserved, symmetric EMTPT

$\displaystyle T^{\mu\nu}=-\frac{2}{\sqrt{-g}}\frac{\delta{\cal L}}{\delta
g_{\mu\nu}}\,.$ (10)

One advantage of the definition in (10) is that the tensor $T^{\mu\nu}$ is
always symmetric and preserves the symmetries of the Lagrangian density. This
is not necessarily true for the definition in equation (3). The resulting EMT
for the Lagrange density in (8) is

$\displaystyle T^{\mu\nu}$ $\displaystyle=$
$\displaystyle\nabla^{\mu}\phi\nabla^{\nu}\phi-g^{\mu\nu}\left[\frac{1}{2}\nabla^{\alpha}\phi\nabla_{\alpha}\phi-\frac{g^{2}}{4!}\phi^{4}\right]+\frac{1}{6}\left[g^{\mu\nu}\Box-\nabla^{\mu}\nabla^{\nu}\right]\varphi^{2}$
(11) $\displaystyle-$
$\displaystyle\frac{1}{6}(R^{\mu\nu}-\frac{1}{2}g^{\mu\nu}R)\varphi^{2}\,.$

It is straightforward to verify that $T^{\mu\nu}$ is traceless and
divergenceless in four dimensions. In flat spacetime where $R=0$ we can write
this result in the form

$\displaystyle T^{\mu\nu}=T_{\rm can}^{\mu\nu}+T_{\rm ex}^{\mu\nu}$
$\displaystyle T_{\rm
ex}^{\mu\nu}=\frac{1}{6}\left[g^{\mu\nu}\Box-\nabla^{\mu}\nabla^{\nu}\right]\varphi^{2}$
(12)

where $T_{\rm can}^{\mu\nu}$ is the canonical EMT (3).

###  II.3 Milne coordinates

Milne coordinates are defined as

$\displaystyle\tau=\sqrt{t^{2}-z^{2}}\text{~{}~{}~{}and~{}~{}~{}}\eta=\frac{1}{2}\ln\left(\frac{t+z}{t-z}\right)$
(13)

which gives $t=\tau\cosh(\eta)$ and $z=\tau\sinh(\eta)$. The metric is

$\displaystyle g=(1,-\tau^{2},-1,-1)_{\rm diag}\,.$ (14)

The coordinate $\eta$ is called the rapidity and depends only on the slope
$v_{z}=\tfrac{z}{t}$. A change of $\eta$ corresponds to a boost by velocity
$v_{z}$. It is therefore also referred to as a boost coordinate.

Milne coordinates are particularly well suited to study a high energy
collision of nuclei moving along the $z$-axis. Surfaces of constant $\eta$
represent timelike surfaces moving away from the origin at constant velocity
and a box with boundaries at fixed $\eta$ expands in the $z$-direction with
time. In this sense Milne coordinates mimic the kinematics of a high energy
collision. In our notation a “dot” indicates a derivative with respect to
$\tau$, and $\vec{\nabla}_{\perp}$ is the transverse gradient operator. The
Riemann tensor (and the Ricci scalar) are both zero and therefore the equation
of motion (9) is the same as the one obtained from the minimally coupled
theory (2). In Milne coordinates it takes the form

$\displaystyle\ddot{\phi}(\tau,\eta,\vec{x}_{\perp})+\frac{1}{\tau}\dot{\phi}-\frac{1}{\tau^{2}}\partial^{2}_{\eta}\phi-\nabla^{2}_{\perp}\phi+\frac{g^{2}}{6}\phi^{3}=0\,.$
(15)

In Appendix A we write all components of the two EMTs (3) and (11) in Milne
coordinates.

##  III Resummation of quantum fluctuations

Observables in the scalar theory exhibit secular divergences at next-to-
leading order if they are calculated in a loop expansion. These divergences
originate from instabilities of the classical solutions. The problem can be
cured using the classical statistical approximation, which amounts to
averaging over a Gaussian ensemble of initial conditions Dusling:2010rm ;
Dusling:2012ig . This procedure implements a resummation scheme that collects
the leading secular terms at each order of an expansion in the coupling
constant444 The resulting EMT has an ultraviolet divergence corresponding to a
vacuum contribution but this can be removed by repeating the calculation with
the background field set to zero and subtracting the results. This vacuum
subtraction has been done for all the calculations presented in this paper..
We briefly describe the structure of the calculation. Further details can be
found in Dusling:2012ig ; Carrington:2022ccm .

###  III.1 Fluctuations

We write the initial field at some small but non-zero555The initial time
$\tau_{0}$ must be small in order to describe a system expanding from
essentially zero volume ($z=0$ at $t=0$) but it cannot be exactly zero due to
the coordinate singularity at $\tau=0$. One can check that the value chosen
for this small initial time does not change the results at finite times.
$\tau=\tau_{0}$ as the sum of a background field contribution,
$\varphi(\tau_{0},\vec{x}_{\perp})$, and a fluctuation,
$\alpha^{(\gamma)}(\tau_{0},\eta,\vec{x}_{\perp})$

$\displaystyle\phi^{(\gamma)}(\tau_{0},\eta,\vec{x}_{\perp})=\varphi(\tau_{0},\vec{x}_{\perp})+\alpha^{(\gamma)}(\tau_{0},\eta,\vec{x}_{\perp})\,.$
(16)

The field $\phi^{(\gamma)}(\tau,\eta,\vec{x}_{\perp})$ at finite proper time
is obtained by solving equation (15) with the initial condition (16). We then
obtain resummed values of observable quantities by averaging over the
ensemble:

$\displaystyle\langle
O(\tau,\eta,x_{\perp})\rangle=\frac{1}{N_{\gamma}}\sum^{N_{\gamma}}_{\gamma=1}O[\phi^{(\gamma)}(\tau,\eta,x_{\perp})]\,.$

We explain below how the two terms in (16) are calculated and the meaning of
the index $\gamma$.

The field $\varphi(\tau_{0},\vec{x}_{\perp})$ is constructed from solutions of
the classical equation of motion (15) at $\tau_{0}$. It is assumed independent
of the spatial rapidity $\eta$. The motivation behind this assumption is the
physical picture of a heavy ion collision in which the nuclei pass through
each other without significant slowing. The resulting velocity distribution
has a property called Bjorken boost-invariance - which is that the
longitudinal velocity $v_{z}$ of frames locally comoving with the fluid is
related to their spacetime position by $v_{z}=z/t$. A fluid with this velocity
distribution will look the same in all longitudinally comoving fluid elements.
We also assume that at very early times the dynamics of the system is
dominated by expansion and we therefore drop the interaction term in the
equation of motion for the purposes of determing
$\varphi(\tau_{0},\vec{x}_{\perp})$. The resulting equation is linear in
$\varphi$ and separating variables we have that at fixed $\tau_{0}$ the
background field can be written as a sum over plane wave mode functions. The
specific combinations of plane waves that we use is discussed in section
III.3.

The fluctuation $\alpha$ that is added to the background field in equation
(16) carries an index $\gamma$ that indicates a Gaussian ensemble of
$N_{\gamma}$ initial fluctuations defined as

$\displaystyle\alpha^{(\gamma)}(\tau_{0},\eta,\vec{x}_{\perp})=\int
dK\big{[}c^{(\gamma)}_{K}a_{K}+c^{(\gamma)*}_{K}a_{K}^{*}\big{]}\,$ (17)

where the index $K$ labels the momentum variables $(\nu,\vec{k}_{\perp})$ that
are conjugate to the coordinate-space variables $(\eta,\vec{x}_{\perp})$,
respectively. The notation $c^{(\gamma)}_{K}$ denotes an element in a Gaussian
distributed ensemble of $N_{\gamma}$ random numbers, with variance

$\displaystyle\langle c^{*}_{K}c_{L}\rangle$ $\displaystyle:=$
$\displaystyle\sum_{\gamma=1}^{N_{\Gamma}}c^{(\gamma)*}_{K}c^{(\gamma)}_{L}=\frac{1}{2}\delta_{KL}\,.$
(18)

The momentum space integration measure is $dK=d\nu
d\vec{k}_{\perp}/(2\pi)^{3}$ and the delta function in equation (18) is
defined so that $\int dK\delta_{KL}=1$. The initial mode functions
$a_{K}\equiv a_{\nu\vec{k}_{\perp}}(\tau_{0},\eta,\vec{x}_{\perp})$ are
obtained from the linearized equations of motion

$\displaystyle\ddot{a}_{K}+\frac{1}{\tau}\dot{a}_{K}-\frac{1}{\tau^{2}}\partial^{2}_{\eta}a_{K}-\Delta_{\perp}a_{K}+\frac{g^{2}}{2}\varphi^{2}(\tau_{0},\vec{x}_{\perp})a_{K}=0\,$
(19)

and normalized so that $\int dK(a_{K},a_{L})=1$ with

$\displaystyle(a_{K},a_{L})=i\tau\int d\eta\int
d^{2}\vec{x}_{\perp}\,\big{(}a_{K}^{*}\partial_{\tau}a_{L}-(\partial_{\tau}a^{*}_{K})a_{L}\big{)}\,.$
(20)

Separating variables and performing the normalization one finds 666The time
dependent part of the equation is second order and has two independent
solutions. We use only the one that has positive frequency behaviour at large
times which is called the second Hankel function.

$\displaystyle a_{K}$ $\displaystyle\equiv$ $\displaystyle
a_{\nu\vec{k}_{\perp}}(\tau_{0},\eta,\vec{x}_{\perp})=\frac{1}{2}\sqrt{\pi}e^{\pi\nu/2}\,e^{i\nu\eta}\chi_{\vec{k}_{\perp}}(\vec{x}_{\perp})H_{i\nu}(\lambda_{\vec{k}_{\perp}}\tau_{0})\,$
(21)

where the functions $\chi_{\vec{k}_{\perp}}(\vec{x}_{\perp})$ are solutions of
the eigenvalue equation

$\displaystyle\left[-\Delta_{\perp}+\frac{g^{2}}{2}\varphi^{2}(\tau_{0},\vec{x}_{\perp})\right]\chi_{\vec{k}_{\perp}}(\vec{x}_{\perp})=\lambda^{2}_{\vec{k}_{\perp}}\chi_{\vec{k}_{\perp}}(\vec{x}_{\perp})\,.$
(22)

When $\tau\to 0$ the Hankel function in (21) oscillates like $e^{\pm
i\tau\nu}$ and the derivative diverges. As explained at the beginning of this
section, we start the evolution at a small positive time $\tau_{0}=10^{-2}$.
To calculate the Hankel function at this initial time we use an asymptotic
series to find it at the smallest time for which the series converges to the
desired accuracy (numerically we truncate the series when successive terms are
less than $10^{-9}$), and then use adaptive fifth order Runge-Kutta to evolve
it to the initial time $\tau_{0}$.

###  III.2 Discretization

We discretize in both directions in the transverse plane with $L$ grid points
and lattice spacing set to one, which means we define all dimensionful
quantities in terms of the transverse lattice grid spacing. The rapidity
variable $\eta$ is discretized with $N$ grid points and lattice spacing $h$.
We consider a unit slice of rapidity with $\eta\in(-1/2,1/2)$ and $h=1/N$. The
discretized version of equation (22) is

$\displaystyle D_{ij;kl}\;\chi_{kl}=\lambda^{2}\chi_{ij}$ (23)

with

$\displaystyle
D_{ij;kl}=(4+V^{\prime\prime}_{ij})\delta_{ik}\delta_{jl}-(\delta_{i+1~{}k}+\delta_{i-1~{}k})\delta_{jl}-\delta_{ik}(\delta_{j+1~{}l}+\delta_{j-1~{}l})\,.$
(24)

Since $D$ is a rank 4 tensor with $L^{4}$ components, we obtain $L$
eigenvalues $\lambda^{e}$ with $e\in(1,L^{2})$ and $L^{2}$ eigenfunctions
$\chi_{ij}^{e}$ which are normalized
$\sum_{ij}\chi^{*e}_{ij}\chi^{\bar{e}}_{ij}=L^{2}\,\delta^{e\bar{e}}$. To
discretize the longitudinal variables we note that the constraint
$\partial_{\eta}^{2}e^{i\nu\eta}=-\nu^{2}e^{i\nu\eta}$ gives

$\displaystyle\varepsilon_{v}:=\nu=\frac{2}{h}\sin\left(\frac{\pi
v}{N}\right)\,$ (25)

and we replace $\nu\to\varepsilon_{v}$ in every factor $e^{\pi\nu/2}$. For the
complex exponential we use $e^{i\nu\eta}\to e^{\frac{2\pi ivn}{N}}$ and the
integral over $\nu$ becomes a sum over $v$ using $\int d\nu/(2\pi)\to
1/(Nh)\sum_{v}^{N}$.

Combining these results we find the discretized versions of equations (16, 18,
21):

$\displaystyle\alpha_{nij}(\tau)=\frac{1}{NL^{2}h}\sum_{v=1}^{N}\sum_{p=1}^{L^{2}}\left[c_{vp}a^{vp}_{nij}(\tau)\,+\text{c.c.}\right]\,$
$\displaystyle a^{vp}_{nij}(\tau)=\frac{1}{2}\sqrt{\pi}e^{\frac{2\pi
ivn}{N}}\,\chi_{ij}^{p}\,e^{\pi\nu/2}H^{(2)}_{i\nu}(\lambda_{\vec{k}_{\perp}}\tau)\,$
$\displaystyle\langle
c_{ve}c^{*}_{u\tilde{e}}\rangle=\frac{1}{2}NL^{2}h\delta_{vu}\delta_{e\tilde{e}}\,.$
(26)

To verify that the discretization is done correctly we have checked the
discretized version of the normalization condition (20).

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

$\displaystyle\nabla_{F}\phi(x)\equiv\phi(i+1)-\phi(i)$
$\displaystyle\nabla_{B}\phi(x)\equiv\phi(i)-\phi(i-1)$
$\displaystyle\sum_{i}f(i)\big{(}\nabla_{F}g(i)\big{)}=-\sum_{i}\big{(}\nabla_{B}f(i)\big{)}g(i)\,.$

We use forward derivatives, defined as $\partial_{x}f(x)\to f(i+1)-f(i)$, so
that the integral of a total derivative term like the last term in equation
(6) will be zero.

###  III.3 Initialization

As explained in section III.1, the initial background field
$\varphi(\tau_{0},\vec{x}_{\perp})$ can be written as a sum of plane wave mode
functions of the form $\cos(\vec{k}_{\perp}\cdot\vec{x}_{\perp})$. In this
section we discuss how to choose these functions and construct the sum.
Ideally we would like to consider initializations that correspond to
collisions of sources with specified radius colliding with a specified impact
parameter, and we could do this by constructing a wave-packet of transverse
plane waves that corresponds to a collision with some specific geometry
777Since the sources that describe the colliding projectiles have support only
on the light cone (15), they do not directly drive the field evolution, but
instead provide information on how to construct this wave-packet. In a Colour
Glass Condensate (CGC) description of a heavy ion collision, the background
field would contain a distribution of momentum modes up to the saturation
momentum.. This calculation is considerably more difficult than what has been
done so far with the method of DEGV. For the purposes of studying the
isotropization of the longitudinal and transverse pressures, it is sufficient
to consider a classical background field that consists of only one mode, which
corresponds to sources with constant surface charge densities with infinite
spatial extent888 In the CGC approach this is called the McLerran-Venugopalan
(MV) model.. In this paper we restrict ourselves to the consideration of
several different initializations involving asymmetric combinations of two
transverse plane waves.

The initial condition that we use for the background field is

$\displaystyle\varphi(\tau_{0},\vec{x}_{\perp})=\varphi_{0}\,\big{(}\cos(k_{x}x+k_{y}y)+A\cos(k_{y}y)\big{)}\,$
(27)

and the time derivative of the background field at $\tau=\tau_{0}$ is set to
zero (numerically $10^{-3}$). The momenta $k_{x}$ and $k_{y}$ are chosen to
take a range of values below the largest lattice eigenvalue, which is
$\sqrt{8}$. The reader will note that our initial classical field is not
always periodic and therefore does not respect our boundary conditions. The
reason is that we want 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
initial boundary conditions is not important.

##  IV Physical observables

###  IV.1 Canonical and conformal EMTs in Minkowski and Milne coordinates

The EMT is conserved

$\displaystyle\nabla_{\mu}T^{\mu\nu}$ $\displaystyle=$
$\displaystyle\partial_{\mu}T^{\mu\nu}+\Gamma^{\mu}_{\mu\sigma}T^{\sigma\nu}+\Gamma^{\nu}_{\mu\sigma}T^{\mu\sigma}=0\,$
(28)

which means that there is a momentum four vector

$\displaystyle{\mathcal{P}}^{\nu}=\int d\sigma_{\mu}T^{\mu\nu},$ (29)

where $d\sigma_{\mu}$ is an element of a spacelike surface $\sigma$, whose
components are conserved by the time evolution providing that the metric
contains no explicit dependence on the coordinates $x^{\mu}$ and the fields
vanish on the boundaries of $\sigma$. This can be seen by noting that in
Minkowski coordinates the integral of the $\nu=0$ component of the
conservation law (28) takes the form:

$\displaystyle\int d^{3}x\frac{\partial}{\partial t}T^{00}=-\int
d^{3}x\left(\partial_{x}T^{x0}+\partial_{y}T^{y0}+\partial_{z}T^{z0}\right)\,.$
(30)

The right side is the integral of a total derivative so assuming the fields
vanish at spatial infinity one finds that the energy defined as $E=\int
d^{3}xT^{00}$ is conserved ($dE/dt=0$).

In Milne coordinates the situation is more complicated due to the $\tau$
dependence of the metric so that $\int d^{3}x\sqrt{-g}T^{00}$ integrated along
a surface of constant $\tau$ is not in general conserved, even for an isolated
system. In Milne coordinates the only non-zero components of the connection
are $\Gamma^{0}_{\eta\eta}=\tau$ and
$\Gamma^{\eta}_{0\eta}=\Gamma^{\eta}_{\eta 0}=1/\tau$. The $\nu=0$ component
of the divergence equation (28) gives

$\displaystyle\frac{\partial T^{00}}{\partial\tau}=-\partial_{\eta}T^{\eta
0}-\frac{1}{\tau}T^{00}-\tau
T^{\eta\eta}-\partial_{x}T^{x0}-\partial_{y}T^{y0}\,.$ (31)

We define the energy density and longitudinal and transverse pressures in
terms of the EMT in Milne coordinates as:

$\displaystyle{\cal E}$ $\displaystyle=$ $\displaystyle\langle T^{00}\rangle$
$\displaystyle p_{L}$ $\displaystyle=$ $\displaystyle\tau^{2}\langle
T^{\eta\eta}\rangle$ $\displaystyle p_{T}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\langle T^{xx}+T^{yy}\rangle\,$ (32)

where the angle brackets indicate an average over rapidity and the transverse
coordinates. Since the average of a total derivative vanishes eq. (31) gives

$\displaystyle\frac{\partial{\cal E}}{\partial\tau}=-\frac{1}{\tau}({\cal
E}+p_{L})\,.$ (33)

We note that when the EMT is written in Milne coordinates $T^{00}$ does not
correspond to the physical energy density as measured in the lab frame except
at zero rapidity. In our calculation we include only a narrow slice of
rapidity centered around mid-rapidity ($\eta=0$). The motivation is that we
are trying to model the gluon fields produced in a relativistic collision of
heavy ions for which the system of fields is largely rapidity independent at
mid-rapidity. The definitions in (32) also involve an average over this thin
slice of rapidity. We therefore take the energy density and pressures as
defined in (32) as the physically relevant ones.

Now we address the fact that the canonical and conformal EMTs are not the same
(even when they are calculated in the same coordinate system). Equation (60)
shows that the difference between the two EMTs, for all components, is a sum
of terms that are either a total derivative in transverse coordinates, a total
derivative in rapidity, or contain the product of the field and its time
derivative (below we call the product $\phi\dot{\phi}$ a decoherence factor).
As explained previously, the total derivative terms do not affect the bulk
properties of the medium because they give zero when they are averaged. The
important point to understand when comparing the two EMTs is that the
decoherence factor goes rapidly to zero at very early times. At
$\tau=\tau_{0}$ we start with a fairly narrow Gaussian distribution of initial
conditions centered on the chosen initial values
$\varphi(\tau_{0},\vec{x}_{\perp})$ and
$\dot{\varphi}(\tau_{0},\vec{x}_{\perp})$. Each initial condition evolves
independently using the equation of motion so that after some characteristic
decoherence time the ensemble average of the product $\phi\dot{\phi}$ goes to
zero. This means that the integrated EMT in (3) is the same as the conformal
form (10) except at very early times. If we only wanted to study late time
properties of the system, like transverse-longitudinal pressure
isotropization, both EMTs would be equally good. Since we are particularly
interested in early time behaviour however, it is important to know which
formalism is better for our purposes.

To understand this we consider the trace of the two EMTs. The trace of the
conformal EMT is zero by construction. The trace of the canonical EMT is

$\displaystyle T^{\mu}_{~{}\mu}$ $\displaystyle=$
$\displaystyle-\nabla^{\mu}\left(\phi\partial_{\mu}\phi\right)=-\frac{d}{d\tau}(\varphi\dot{\varphi})-\frac{\varphi\dot{\varphi}}{\tau}+\vec{\nabla}_{\perp}\cdot(\varphi\vec{\nabla}_{\perp}\varphi)\,.$
(34)

This expression has terms that are proportional to the decoherence factor and
its derivative, in addition to total derivatives, which means that the trace
of the averaged canonical EMT is not zero at early times. We are interested in
assessing the extent to which our system of fields behaves hydrodynamically,
and the traditional formulation of hydrodynamics based on the assumption of
local equilibrium assumes a traceless EMT (see section IV.3). This suggests we
want to work with the conformal definition of the EMT. Furthermore, in section
V we present some results that show that at early times the time dependence of
the canonical EMT is characterized by rapid oscillations created by
decoherence terms. These oscillations disappear at larger times but make it
impossible to study early time dynamics. They are not present when the
conformal EMT is used.

Finally we make the formal argument that we should work with a conformally
invariant theory since we hope to use the scalar system to understand some
properties of glasma.

###  IV.2 Spatial and momentum asymmetries

A correlation between initial spatial transverse asymmetry and azimuthal
asymmetry in the momentum field is generally considered characteristic of the
onset of hydrodynamic behaviour. Physically the idea is that if pressure
gradients are converted into fluid velocities, then momentum anisotropy will
grow as the spatial anisotropy decreases. Spatial deviations from azimuthal
symmetry can be characterized with the quantity Kolb:1999it ; Blaizot:2014nia

$\displaystyle\varepsilon=-\frac{\int
d^{2}x_{\perp}\left(\frac{x^{2}-y^{2}}{\sqrt{x^{2}+y^{2}}}\right)T^{00}}{\int
d^{2}x_{\perp}{\sqrt{x^{2}+y^{2}}}T^{00}}$ (35)

where we have used $\vec{x}_{\perp}=(x,y)$. In a system of fields, momentum
anisotropy can be described by an elliptic flow coefficient Carrington:2021qvi

$\displaystyle v_{2}=\frac{\int
d^{2}x_{\perp}\left(\frac{(p^{x})^{2}-(p^{y})^{2}}{\sqrt{(p^{x})^{2}+(p^{y})^{2}}}\right)}{\int
d^{2}x_{\perp}{\sqrt{(p^{x})^{2}+(p^{y})^{2}}}}$ (36)

where $p^{x}=T^{0x}$ and $p^{y}=T^{0y}$. We comment that this definition is
not directly related to what is normally measured in heavy-ion collisions,
where one looks at the momentum distribution of produced particles. A
different quantity that is commonly calculated in heavy-ion physics is the
transverse pressure anisotropy defined

$\displaystyle A_{xy}=\frac{\int d^{2}\vec{x}_{\perp}\,(T^{xx}-T^{yy})}{\int
d^{2}\vec{x}_{\perp}\,(T^{yy}+T^{xx})}\,.$ (37)

Hydrodynamic calculations show that $A_{xy}$ is closely related to the
elliptic flow coefficient of produced particles Kolb:1999it  and the
transverse pressure anisotropy has been used extensively in the literature to
characterize momentum anisotropy Luzum:2008cw ; Luzum:2009sb .

###  IV.3 Universal flow

The authors of ref. Vredevoogd:2008id  derived an equation that describes the
way that transverse flow develops for any system for which the longitudinal
flow is boost invariant, the EMT is traceless, and the anisotropy of the
spatial components of the EMT depends more strongly on $\tau$ than on the
transverse coordinates at early times. In our calculation the first condition
is satisfied by construction, the second is always true if we use the
conformal EMT, and the third is well satisfied in our calculation.

We define the $x$-component of a velocity and an effective acceleration

$\displaystyle V_{x}=\frac{T^{0x}}{T^{00}}$
$\displaystyle\alpha=\frac{\partial
V_{x}}{\partial\tau}=\frac{\partial}{\partial\tau}\left(\frac{T^{0x}}{T^{00}}\right)=\frac{\partial_{\tau}T^{0x}}{T^{00}}-\frac{T^{0x}}{T^{00}}\frac{\partial_{\tau}T^{00}}{T^{00}}\,.$
(38)

We can rewrite this in a more convenient form as follows. The zeroth component
of the equation $\nabla_{\mu}T^{\mu\nu}$ in Milne coordinates is given in
equation (31). The $x$-component is

$\displaystyle\frac{\partial
T^{0x}}{\partial\tau}=-\frac{1}{\tau}T^{0x}-\partial_{x}T^{xx}-\partial_{y}T^{yx}\,$
(39)

where we have dropped the derivative with respect to $\eta$ since we are
interested in a system that it at least approximately boost invariant. We will
assume that $T^{yx}$ is small and set it to zero. We rewrite $\alpha$ using
equations (31, 39) to obtain

$\displaystyle T^{00}\alpha$ $\displaystyle=$
$\displaystyle\partial_{\tau}T^{0x}-\frac{T^{0x}}{T^{00}}\partial_{\tau}T^{00}$
(40) $\displaystyle=$
$\displaystyle-\left(\frac{1}{\tau}T^{0x}+\partial_{x}T^{xx}\right)+\frac{T^{0x}}{T^{00}}\left(\frac{1}{\tau}T^{00}+\tau
T^{\eta\eta}+\partial_{x}T^{x0}+\partial_{y}T^{y0}\right)\,.$

We assume that at early times the field $\phi$ has a power series expansion in
$\tau$ and define

$\displaystyle\phi=\sum_{i=0}\phi_{2i}\tau^{2i}$ (41)

where only even powers of $\tau$ are included because the equation of motion
does not allow solutions with odd powers. From the explicit expressions for
the components of the EMT (59, 60) we see that the diagonal components are
series with only even powers of $\tau$ and both $T^{x0}$ and $T^{y0}$ have
only odd powers. We define $T^{00}=T^{00}_{0}+T^{00}_{2}\tau^{2}+\dots$ where
the subscript indicates the power of $\tau$ associated with each coefficient,
and similarly for each component of the EMT.

We write $\alpha=\alpha_{0}+\dots$ where $\alpha_{0}=T^{0x}_{1}/T^{00}_{0}$
and keep terms in (40) to order $\tau^{0}$ which gives

$\displaystyle
T^{0x}_{1}=-\partial_{x}T^{xx}_{0}-T_{1}^{0x}+\frac{T_{1}^{0x}}{T_{0}^{00}}(T^{00}_{0}+T^{\eta\eta}_{-2})\,.$
(42)

Rearranging we have

$\displaystyle
T_{1}^{0x}(T^{00}_{0}-T^{\eta\eta}_{-2})=-T^{00}_{0}\partial_{x}T^{xx}_{0}\,.$
(43)

The condition that the trace is zero, to zeroth order in $\tau$, is
$T^{00}_{0}-T^{\eta\eta}_{-2}-T^{xx}-T^{yy}=0$. Equation (43) therefore takes
the form

$\displaystyle
T_{1}^{0x}=-\frac{T^{00}_{0}}{T^{xx}_{0}+T^{yy}_{0}}\partial_{x}T^{xx}_{0}\,.$
(44)

We assume that at early times $T^{xx}\approx T^{yy}$ and that $T^{00}/T^{xx}$
is approximately independent of $\vec{x}_{\perp}$. Using these assumptions and
multiplying by $\tau$ equation (44) gives

$\displaystyle T^{0x}=-\frac{\tau}{2}\partial_{x}T^{00}+{\cal O}(\tau^{2})\,$
(45)

or equivalently

$\displaystyle\frac{\partial
v_{x}}{\partial\tau}=-\frac{1}{2}\frac{\partial_{x}T^{00}}{T^{00}}+{\cal
O}(\tau^{2})\,.$ (46)

This result tells us that at early times the impulse to the transverse
collective flow is connected to the gradient of the energy density. It means
that different systems with the same energy density but very different
pressures, ranging from hydrodynamic $p_{T}=p_{L}={\cal E}/3$ to systems with
highly asymmetric initial pressures like ours, will develop transverse flow in
the same way. In section (V) we show that (45) is well satisfied in our
calculation, at early times.

For comparison we discuss how (45) can be derived for a system governed by
hydrodynamics. The hydrodynamic definition of the EMT in Minkowski coordinates
is

$\displaystyle T^{\mu\nu}_{\rm
hydro}=(\epsilon+p)u^{\mu}u^{\nu}-pg^{\mu\nu}\,$ (47)

where $\epsilon$ is the energy density, $p$ is the thermodynamic pressure, and
$u^{\mu}=\gamma(1,\vec{v})$ is the fluid velocity. The tracelessness of the
EMT gives the equation of state $p=\epsilon/3$. If we assume that each
velocity component can be written as a power series in $t$

$\displaystyle v^{i}=v^{i}_{0}+v^{i}_{1}t+v^{i}_{2}t^{2}+\dots$ (48)

then at leading order we find

$\displaystyle T^{0x}=-t\frac{\partial T^{xx}}{\partial x}+{\cal O}(t^{2})\,$
(49)

and since $T^{xx}=p=\epsilon/3$ we have

$\displaystyle T^{0x}=-\frac{t}{3}\frac{\partial T^{00}}{\partial x}+{\cal
O}(t^{2})\,.$ (50)

The same calculation in two spatial dimensions gives

$\displaystyle T^{x0}=-\frac{t}{2}\frac{\partial T^{00}}{\partial x}+{\cal
O}(t^{2})\,$ (51)

which agrees with (45) at mid-rapidity.

##  V Results

Our goal is to study how the system isotropizes. We are particularly
interested in understanding the extent to which the system behaves
hydrodynamically, and how hydrodynamic behaviour develops with time. Unless
stated otherwise, all calculations are done with
$(N,L,N_{\gamma})=(121,41,256)$. We have checked that results are not
sensitive to these numbers. The initialization is specified by giving the
values of $(\varphi_{0},A,k_{x},k_{y})$ that are used in equation (27). The
coupling is $g=4$ in all calculations.

###  V.1 Trace of the canonical EMT

The trace of the averaged canonical EMT is not zero, as shown in equation (6),
but it goes rapidly to very small values as the system evolves because the
decoherence factor goes to zero. In fig. 1 we show the energy and the sum of
the pressures for the canonical EMT. The blue curve, which shows the sum of
the pressures, oscillates around the red curve, which shows the energy, until
$\tau\approx 20$. At larger times the trace (the difference between the blue
and red lines) is numerically very close to zero. For the conformal EMT the
sum of the pressures and the energy are right on top of the canonical energy
curve. The oscillations in the pressures from the canonical EMT are typical of
the early time behaviour and caused by the decoherence terms.

Figure 1: The energy (red) and the sum $p_{L}+2p_{T}$ (blue) for the canonical
EMT. For the conformal EMT the sum of the pressures and the energy are right
on top of the canonical energy curve. The calculation is done with
$(\varphi_{0},A)=(15,0)$ and $k_{x}=k_{y}=0.41$.

###  V.2 Transverse-Longitudinal Isotropization

A simple way to study how the system comes to equilibrium is to see how
closely it satisfies the isotropic equation of state: $p_{L}=p_{T}={\cal
E}/3$. This was studied by DEGV who showed that $(p_{T}-p_{L})/{\cal E}$
decreases as a function of time999 Note that because we work numerically with
a box of finite size in coordinate space, the calculation will break down at
large times because there will be highly occupied momentum modes that are not
supported.. The result is important because it proves that the resummation
method developed by DEGV is successful in capturing the dominant physics of
the expanding plasma.

Figure 2 shows the transverse and longitudinal pressures normalized by the
energy from the canonical EMT, for two initializations that have the same
initial field amplitude. The figure shows that the isotropization of the
transverse and longitudinal pressures is not significantly affected when the
second cosine is included in the initialization.

Figure 2: The transverse and longitudinal pressures normalized by the energy
from the canonical EMT for the initializations $(\varphi_{0},A)=(15,0)$ and
$(\varphi_{0},A)=(15/\sqrt{2},1)$. Both calculations are done with
$k_{x}=k_{y}=0.41$. In both cases the upper curves are $p_{T}/{\cal E}$ and
the lower lines are $p_{L}/{\cal E}$.

In fig. 3 we show the transverse and longitudinal pressures normalized by the
energy for the canonical and conformal EMTs. The blue curves show the
oscillations that are characteristic of the canonical EMT, which die off at
larger times when the decoherence terms $\sim\phi\dot{\phi}$ go to zero. The
red curves are the results from the conformal EMT and are much smoother at
early times. At $\tau\gtrsim 40$ the results from the two EMTs are
indistinguishable.

Figure 3: The normalized transverse and longitudinal pressures for the
canonical EMT (blue) and the conformal EMT (red). The initialization is
$(\varphi_{0},A)=(15,0)$ and $k_{x}=k_{y}=0.41$. The upper and lower lines are
respectively $p_{T}/{\cal E}$ and $p_{L}/{\cal E}$.

###  V.3 Azimuthal asymmetry

When we use an initialization with $A\neq 0$ the initial fields are
azimuthally asymmetric. We want to study how the azimuthal asymmetry of the
energy density and momentum distribution develop in time, and how they are
related to each other. The asymmetry in the transverse plane is orders of
magnitude smaller than the transverse-longitudinal asymmetry discussed in the
previous subsection, which makes it hard to study numerically. In fig. 4 we
show the eccentricity (left panel) and elliptic flow coefficient (right panel)
as functions of $\tau$ for initializations $(\varphi_{0},A,k_{x})=(15,1,2.83)$
and $k_{y}\in(1.30,1.47,2.83)$. One sees that the initial spatial eccentricity
decreases with time, and that the size of the initial eccentricity is
correlated with the size of the momentum anisotropy that is produced. The
elliptic flow coefficient initially increases and eventually decays. In fig. 5
we show the similarity in the shape of the two measures $v_{2}$ and $A_{xy}$
defined in equations (36, 37). The left panel shows $A_{xy}$ for the same
three initializations as in fig. 4. The right panel shows the elliptic flow
coefficient, and the transverse pressure asymmetry shifted left by 1.7 time
units and multiplied by 2.4.

Figure 4: The left panel shows the eccentricity as a function of $\tau$ and
the right panel is the elliptic flow coefficient $v_{2}$. The initialization
is $(\varphi_{0},A,k_{x})=(15,1,1.47)$ and from bottom to top
$k_{y}=(1.12\text{~{}(red)},1.30\text{~{}(green)},1.47\text{~{}(blue)})$. The
calculations are done with $(N,L,N_{\gamma})=(81,31,256)$.

Figure 5: The left panel is the transverse pressure asymmetry in equation (37)
for the same initializations as in fig. 4. The right panel is $v_{2}$
superimposed on a plot of $A_{xy}$ that has been shifted left 1.7 time units
and multiplied by a factor of 2.4. The lighter coloured lines show the shifted
results for $A_{xy}$.

###  V.4 Universal flow

We can check how well equation (45) is satisfied for our system of scalar
fields. The calculation is numerically difficult because the functions on both
sides of the equation fluctuate in transverse position at early times, and
decay rapidly as time increases. To compare them we calculate the root of the
average of the squares of both sides of the equation. We define

$\displaystyle h_{1}=\sqrt{\left\langle\left(\frac{\tau}{2}\frac{\partial
T^{00}}{\partial
x}\right)^{2}\right\rangle}\text{~{}~{}and~{}~{}}h_{2}=\sqrt{\langle(T^{0x})^{2}\rangle}$
(52)

where the angle brackets indicate an average over spatial rapidity and the
transverse coordinates101010The factor $\tau$ in $h_{1}$ is multiplied by the
length scale $1/L$ because we have set the lattice spacing in the transverse
direction to one.. Figure 6 shows there is good agreement of the two
quantities in (52) at early times.

Figure 6: The two quantities in equation (52) at early times. The left panel
is a linear plot and the right is a log-log graph that shows the differences
at later times more clearly. The calculation is done with
$(\varphi_{0},A)=(15,0)$ and $k_{x}=k_{y}=0.41$.

##  VI Conclusions

In this paper we have studied the dynamics of a self-interacting real scalar
field. We work with a conformally invariant action and calculate the energy
momentum tensor using the classical statistical approximation. We have studied
the time evolution of the azimuthal asymmetry of the energy density and the
momentum field and shown that the two quantities are correlated. The
eccentricity decreases and the momentum elliptic flow coefficient initially
increases and then later decays. We have studied the asymmetry of the
transverse pressures and shown that its behaviour is strikingly similar to
that of the elliptic flow coefficient. Both of these behaviours are seen in
hydrodynamic calculations that describe the physics of a relativistic heavy-
ion collision at much later times. We have compared the derivative of the
transverse velocity and the gradient of the energy density and shown that the
characteristics of universal flow are present in our system.

Our results provide further support for the use of the classical statistical
approximation to describe a system of strongly coupled fields and suggest that
the method is worthy of further investigation using more realistic initial
configurations, and physical theories that involve additional fields. Such
calculations could lead to a deeper understanding of the early onset of
hydrodynamic behaviour that is seen in heavy ion collisions.

##  Appendix A Useful equations

In Milne coordinates the only non-zero components of the connection are
$\Gamma^{0}_{\eta\eta}=\tau$ and $\Gamma^{\eta}_{0\eta}=\Gamma^{\eta}_{\eta
0}=1/\tau$. The transformation from Minkowski to Milne coordinates is

$\displaystyle A^{\mu}{}_{\nu}\equiv\frac{\partial x^{\mu}}{\partial x^{\nu}}$
$\displaystyle=$ $\displaystyle\left[\begin{tabular}[]{c c c
c}$\cosh{\eta}$&$-\sinh(\eta)$&0&0\\\
$-\frac{1}{\tau}\sinh(\eta)$&$\frac{1}{\tau}\cosh(\eta)$&0&0\\\ 0&0&-1&0\\\
0&0&0&-1\\\ \end{tabular}\right]\,.$ (57)

Some results for the gradient operator are

$\displaystyle\partial_{\mu}=(\partial_{\tau},\partial_{\eta},\vec{\nabla})$
$\displaystyle\partial^{\mu}=(\partial_{\tau},-\frac{1}{\tau^{2}}\partial_{\eta},-\vec{\nabla})$
$\displaystyle(\partial_{\mu}\varphi)(\partial^{\mu}\varphi)=\dot{\varphi}^{2}-\frac{1}{\tau^{2}}(\partial_{\eta}\varphi)^{2}-(\vec{\nabla}\varphi)^{2}$
$\displaystyle\Box=\partial_{\tau}^{2}+\frac{1}{\tau}\partial_{\tau}-\frac{1}{\tau^{2}}\partial_{\eta}^{2}-\vec{\nabla}^{2}\,.$
(58)

In equations (59, 60) we give the components of the canonical EMT (3) and the
additional piece that appears in the conformal EMT (12). For the diagonal
terms we give two ways to write the components of $T^{\mu\nu}_{\rm ex}$. The
label ‘alt’ indicates that the equation of motion has been used. These
expressions are easier to calculate numerically. The components of the
canonical EMT are

$\displaystyle T_{\rm
can}^{00}=\frac{1}{2}\left((\dot{\phi})^{2}+\frac{(\partial_{\eta}\phi)^{2}}{\tau^{2}}+(\partial_{x}\phi)^{2}+(\partial_{y}\phi)^{2}\right)+\frac{g^{2}}{4!}\phi^{4}$
$\displaystyle T_{\rm
can}^{\eta\eta}=\frac{1}{2}\tau^{-4}(\partial_{\eta}\phi)^{2}+\frac{1}{\tau^{2}}\left[\frac{1}{2}\left((\dot{\phi})^{2}-(\partial_{x}\phi)^{2}-(\partial_{y}\phi)^{2}\right)-\frac{g^{2}}{4!}\phi^{4}\right]$
$\displaystyle T_{\rm
can}^{xx}=\left[\frac{1}{2}\left((\dot{\phi})^{2}-\frac{(\partial_{\eta}\phi)^{2}}{\tau^{2}}+(\partial_{x}\phi)^{2}-(\partial_{y}\phi)^{2}\right)-\frac{g^{2}}{4!}\phi^{4}\right]$
$\displaystyle T_{\rm
can}^{yy}=\left[\frac{1}{2}\left((\dot{\phi})^{2}-\frac{(\partial_{\eta}\phi)^{2}}{\tau^{2}}-(\partial_{x}\phi)^{2}+(\partial_{y}\phi)^{2}\right)-\frac{\lambda}{4!}\phi^{4}\right]$
$\displaystyle T_{\rm can}^{xy}=(\partial_{x}\phi)(\partial_{y}\phi)$
$\displaystyle T_{\rm
can}^{0\eta}=-\frac{1}{\tau^{2}}\dot{\phi}(\partial_{\eta}\phi)$
$\displaystyle T_{\rm can}^{0y}=-\dot{\phi}(\partial_{y}\phi)$ $\displaystyle
T_{\rm can}^{0x}=-\dot{\phi}(\partial_{x}\phi)\,.$ (59)

The components of $T^{\mu\nu}_{\rm ex}$ are

$\displaystyle T_{\rm
ex}^{00}=\frac{1}{3}\frac{1}{\tau}(\varphi\dot{\varphi})-\frac{1}{3}\frac{1}{\tau^{2}}\big{(}(\partial_{\eta}\varphi)^{2}+\varphi\partial_{\eta}^{2}\varphi\big{)}-\frac{1}{3}\vec{\nabla}_{\perp}\cdot(\varphi\vec{\nabla}_{\perp}\varphi)$
$\displaystyle T_{\rm ex-
alt}^{00}=-\frac{1}{3}\left[\phi\ddot{\phi}+\frac{g^{2}}{6}\phi^{4}+\frac{(\partial_{\eta}\phi)^{2}}{\tau^{2}}+|\nabla_{\perp}\phi|^{2}\right]$
$\displaystyle T^{\eta\eta}_{\rm
ex}=-\frac{1}{3\tau^{2}}\frac{\partial}{\partial\tau}(\varphi\dot{\varphi})+\frac{1}{3\tau^{2}}\vec{\nabla}_{\perp}\cdot(\varphi\vec{\nabla}_{\perp}\varphi)$
$\displaystyle T_{\rm ex-
alt}^{\eta\eta}=-\frac{1}{3\tau^{4}}\left[\phi\partial_{\eta}^{2}\phi-\tau\phi\dot{\phi}-\frac{g^{2}}{6}\phi^{4}\tau^{2}+\tau^{2}\left(\dot{\phi}^{2}-|\nabla_{\perp}\phi|^{2}\right)\right]$
$\displaystyle T^{xx}_{\rm
ex}=-\frac{1}{3}\left(\frac{\partial}{\partial\tau}(\varphi\dot{\varphi})+\frac{1}{\tau}(\varphi\dot{\varphi})-\frac{1}{\tau^{2}}\partial_{\eta}(\varphi\partial_{\eta}\varphi)\right)+\frac{1}{3}\frac{\partial}{\partial
y}\left(\phi\frac{\partial\phi}{\partial y}\right)\,$ $\displaystyle T_{\rm
ex-
alt}^{xx}=-\frac{1}{3}\left[\phi\partial_{x}^{2}\phi-\frac{g^{2}}{6}\phi^{4}+(\dot{\phi})^{2}-\frac{(\partial_{\eta}\phi)^{2}}{\tau^{2}}-(\partial_{y}\phi)^{2}\right]$
$\displaystyle T^{yy}_{\rm
ex}=-\frac{1}{3}\left(\frac{\partial}{\partial\tau}(\varphi\dot{\varphi})+\frac{1}{\tau}(\varphi\dot{\varphi})-\frac{1}{\tau^{2}}\partial_{\eta}(\varphi\partial_{\eta}\varphi)\right)+\frac{1}{3}\frac{\partial}{\partial
x}\left(\phi\frac{\partial\phi}{\partial x}\right)\,$ $\displaystyle T_{\rm
ex-
alt}^{yy}=-\frac{1}{3}\left[\phi\partial_{y}^{2}\phi-\frac{g^{2}}{6}\phi^{4}+(\dot{\phi})^{2}-\frac{(\partial_{\eta}\phi)^{2}}{\tau^{2}}-(\partial_{x}\phi)^{2}\right]$
$\displaystyle T_{\rm
ex}^{0\eta}=\frac{1}{3\tau^{2}}\left(\phi\partial_{\eta}\dot{\phi}+\dot{\phi}\partial_{\eta}\phi-\frac{1}{\tau}\phi\partial_{\eta}\phi\right)$
$\displaystyle T_{\rm
ex}^{0x}=\frac{1}{3}(\phi\partial_{x}\dot{\phi}+\dot{\phi}\partial_{x}\phi)$
$\displaystyle T_{\rm
ex}^{0y}=\frac{1}{3}(\phi\partial_{y}\dot{\phi}+\dot{\phi}\partial_{y}\phi)$
$\displaystyle T_{\rm
ex}^{xy}=-\frac{1}{3}\left[\phi(\partial_{x}\partial_{y}\phi)+(\partial_{x}\phi)(\partial_{y}\phi)\right]\,.$
(60)

###### Acknowledgements.

Margaret Carrington gratefully acknowledges helpful discussions with François
Gelis. This work has been supported by the Natural Sciences and Engineering
Research Council of Canada Discovery Grant program from grants 2023-00023 and
2018-04090. This research was enabled in part by support provided by WestGrid
(www.westgrid.ca) and the Digital Research Alliance of Canada
(alliancecan.ca).

## References

  * (1) K. Dusling, T. Epelbaum, F. Gelis and R. Venugopalan, Nucl. Phys. A 850, 69, (2011).
  * (2) T. Epelbaum and F. Gelis, Nucl. Phys. A 872, 210 (2011).
  * (3) K. Dusling, T. Epelbaum, F. Gelis and R. Venugopalan, Phys. Rev. D 86, 085040, (2012).
  * (4) M. E. Carrington, G. Kunstatter, C. D. Phillips and M. E. Rubio, Entropy 24, 1612 (2022).
  * (5) J. Berges, K. Boguslavski, S. Schlichting and R. Venugopalan, JHEP 05, 054, (2014).
  * (6) T. Epelbaum, F. Gelis and B. Wu, Phys. Rev. D 90, 065029, (2014).
  * (7) F. Gelis and S. Hauksson, [arXiv:2403.11908 [hep-ph]].
  * (8) J. Vredevoogd and S. Pratt, Phys. Rev. C 79, 044915, (2009).
  * (9) D. N. Blaschke, F. Gieres, M. Reboud and M. Schweda, Nucl. Phys. B 912, 192, (2016).
  * (10) Leonard Parker and David Toms, Quantum Field Theory in Curved Spacetime, Cambridge University Press, 2009.
  * (11) P. F. Kolb, J. Sollfrank and U. W. Heinz, Phys. Lett. B 459, 667 (1999).
  * (12) J. P. Blaizot, W. Broniowski and J. Y. Ollitrault, Phys. Lett. B 738, 166 (2014).
  * (13) M. E. Carrington, A. Czajka and S. Mrówczyński, Phys. Rev. C 106, 034904, (2022).
  * (14) J.Y. Ollitrault, Phys. Rev. D 46, 229, (1992).
  * (15) M. Luzum and P. Romatschke, Phys. Rev. C 78, 034915, (2008).
  * (16) M. Luzum and P. Romatschke, Phys. Rev. Lett. 103, 262302, (2009).