From the colour glass condensate to filamentation: systematics of classical Yang–Mills theory

The non-equilibrium early time evolution of an ultra-relativistic heavy ion collision is often described by classical lattice Yang–Mills theory, starting from the colour glass condensate (CGC) effective theory with an anisotropic energy momentum tensor as initial condition. In this work we investigate the systematics associated with such studies and their dependence on various model parameters (IR, UV cutoffs and the amplitude of quantum fluctuations) which are not yet fixed by experiment. We perform calculations for SU(2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2$$\end{document}) and SU(3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3$$\end{document}), both in a static box and in an expanding geometry. Generally, the dependence on model parameters is found to be much larger than that on technical parameters like the number of colours, boundary conditions or the lattice spacing. In a static box, all setups lead to isotropisation through chromo-Weibel instabilities, which is illustrated by the accompanying filamentation of the energy density. However, the associated time scale depends strongly on the model parameters and in all cases is longer than the phenomenologically expected one. In the expanding system, no isotropisation is observed for any parameter choice. We show how investigations at fixed initial energy density can be used to better constrain some of the model parameters.


Introduction
The medium created by ultra-relativistic heavy-ion collisions is characterised by strong collective behaviour. It is generally accepted that a quark-gluon plasma (QGP) is formed and the effective theory describing the multiparticle correlations of this nearly-perfect fluid is relativistic viscous hydrodynamics. For a long time it was believed that the application of a e-mail: philipsen@th.physik.uni-frankfurt.de b e-mail: wagenbach@th.physik.uni-frankfurt.de c e-mail: zafeiropoulos@thphys.uni-heidelberg.de hydrodynamic models requires the thermalisation time from the initial non-equilibrium stage of the collision to the QGP to be very short [1,2] compared to the lifetime of the QGP. More recently it was argued that hydrodynamics also applies to a not-yet equilibrated system [3].
From a theoretical point of view, a heavy-ion collision has different stages. As an initial condition, one assumes the colour glass condensate (CGC), i.e. an effective field theory description of boosted, saturated gluons [4]. The resulting strong gauge field dynamics constitutes the first stage of the evolution. A following later stage is then governed by hydrodynamic equations until the medium becomes too dilute for this long wavelength description. The precise duration of the early stage is not yet known for realistic values of the coupling. Models of the hydrodynamical stage constrain it to be around or less than 1 fm [5].
In this work, we focus on the early time dynamics of the gauge fields out of equilibrium, where we pursue a purely classical treatment of Yang-Mills theory. This approach is justified for the infrared modes of gauge fields with a high occupation number, see, e.g., [29][30][31][32].
Our goal is to initiate a systematic study of the dependence on a variety of parameters entering through the CGC initial condition as well as the systematics of the classical evolution itself. In particular, we compare a treatment of the realistic SU(3) gauge group with the more economical SU (2), monitor a gauge-invariant definition of the occupation number of field modes to address the validity of the classical approximation, and compare the evolution in a static box with the one in an expanding medium. We also attempt to quantify the dependence of our results on various model parameters introduced in the literature, like the amplitude of initial boost non-invariant fluctuations, an IR cutoff to emulate colour neutrality on the scale of nucleons as well as a UV cutoff on the initial momentum distribution. Many of these issues have already been addressed one by one when they were introduced, as indicated in the following sections, but not in their interplay, as we attempt to do here.
In the next section we summarise the theoretical framework of our approach and give the CGC initial conditions this work is based on. In Sect. 3, we present the numerical results of our simulations, where we extensively elaborate on the underlying parameter space of the CGC. We will see that the system is highly sensitive to the model parameters and suggest a method to reduce the number of free parameters by keeping the system's physical energy density fixed. We also present depictions of the filamentation of the energy density in position space, which results from initial quantum fluctuations and indicates the occurrence of chromo-Weibel instabilities. Section 4 contains our conclusions and an outlook. Some very early stages of this work appeared as a conference proceeding [33].
The anisotropy parameter ξ = a σ /a t is the ratio of spatial and temporal lattice spacings which does not renormalise in the classical limit, and β = 2N c /g 2 is the lattice gauge coupling (we choose N c = 2 and N c = 3 colours).
In the expanding geometry, where we use comoving coordinates τ = √ t 2 − z 2 and η = atanh(z/t), the lattice action reads We introduced the transverse lattice spacing a ⊥ and the dimensionless rapidity discretisation a η . Inserting the link variables , and expanding around small values of the lattice spacing one recovers the classical Yang-Mills action in the continuum limit, a μ → 0. In order to choose canonical field variables and construct a Hamiltonian, we set i.e., we are using temporal gauge. The field variables are then the spatial (and rapidity) links and the rescaled dimensionless chromo-electric fields, For the situation in a static box this results in the standard Hamiltonian with corresponding classical field equations and Gauss constraint For the expanding case we have, in comoving coordinates, Re Tr 2τ 1 − U xy with field equations and Gauss constraint We then consider the time evolution of the classical statistical system whose equilibrium states are determined by the classical partition function For simulations in equilibrium, initial configurations are generated with a thermal distribution governed by this partition function, and then evolved in t by solving (10) or (13), respectively. For a system out of equilibrium, by definition there is no partition function. Rather, specific field configurations satisfying the Gauss constraint have to be given by some initial conditions, and are then evolved using the field equations.

Non-equilibrium initial conditions (CGC)
Heavy-ion collisions at high energy density can be described in terms of deep inelastic scattering of partons. The corresponding parton distribution functions are dominated by gluonic contributions, which motivates the description in terms of a colour glass effective theory [4,34]. The gluonic contribution to the parton distribution is limited by a saturation momentum Q s , which is proportional to the collision energy. When the saturation scale Q s becomes large there is a time frame where soft and hard modes get separated [35]. The colliding nuclei constitute hard colour sources, which can be seen as static. Due to time dilatation, they are described as thin sheets of colour charge.
Choosing z as the direction of the collision, this is usually described in light cone coordinates, The colour charges are distributed randomly from collision to collision. In the McLerran-Venugopalan (MV) model [36] the distribution is taken to be Gaussian, with charge densities Here μ 2 ∼ A 1/3 fm −2 is the colour charge squared per unit area in one colliding nucleus with atomic number A. It is non-trivially related to the saturation scale [37], with Q s ≈ Q := g 2 μ. For Pb-Pb or Au-Au collisions, this is larger than the fundamental QCD scale Λ QCD . We choose a value in the range of expectations for ultra-relativistic heavy-ion collision at the Large Hadron Collider (Q s ≈ 2 GeV [38]) and fix Q = 2 GeV for our simulations throughout this paper. Originally the MV model was formulated for a fixed time slice. Later it was realised that, in order to maintain gaugecovariance in the longitudinal direction, this initial time slice has to be viewed as a short-time limit of a construction using N l time slices, containing Wilson lines in the longitudinal direction [37,39]. In the literature the designation "N y " is also frequently used for the number of longitudinal sheets, but in order to distinguish it from the lattice extent in ydirection we use N l instead.
The colour charge densities produce the non-Abelian current and the corresponding classical gluon fields are then obtained by solving the Yang-Mills equations in the presence of those sources, For the lattice implementation of this initial condition, we follow [37] and solve with the lattice Laplacian in the transverse plane, The two nuclei are labelled by k = 1, 2, the index v = 1, . . . , N l indicates the transverse slice under consideration and m is an IR regulator. For m = 0, a finite lattice volume acts as an effective IR cutoff. However, a finite m ∼ Λ QCD is expected to exist, since correlators of colour sources are screened over distances of Λ −1 QCD , as was initially proposed in [37]. Of course, a determination of this screening length requires the full quantum theory and thus is beyond a classical treatment. We shall investigate the dependence of our results by varying m between zero and some value of the expected order of magnitude. Physically, the parameter m indicates the inverse length scale over which objects are colour neutral in our description, and hence m = 0.1 Q ≈ 200 MeV ≈ 1 fm −1 ≈ 1 R p , with R p being the proton radius, is a sensible choice.
Although we already have a UV cutoff ∼ 1/a ⊥ from the lattice discretisation, often an additional UV momentum cutoff Λ is used in the literature [19,29,40,41]. It is implemented by neglecting all modes larger than Λ while solving Poisson's equation (20a) in momentum space. There are two ways to interpret this additional UV cutoff. It is sometimes used as a technical trick to maintain an initial spectrum in the IR while allowing to make a ⊥ smaller, in order to reduce discretisation effects. As we shall see, this is only consistent in the expanding scenario. Alternatively, it can be interpreted simply as an additional model parameter of the CGC, which restricts the colour sources in Fourier space to modes in the IR. Again, we shall investigate how results depend on the presence and size of this parameter.
To get the transverse components of the collective initial lattice gauge fields U k = exp(iα a T a ), α a ∈ R, we have to solve N g equations at each point on the transverse plane, Tr T a (U (1) For the case of N c = 3 we do this numerically using multidimensional root finding methods of the GSL library [42]. For the case of N c = 2, one can find a closed-form expression and circumvent this procedure, i.e. (22) reduces to The remaining field components are with the index convention introduced in Sect. 2.1.
To make the initial conditions more realistic, fluctuations can be added on top of this background [15,43], which are supposed to represent quantum corrections to the purely classical fields. They are low momentum modes constructed to satisfy the Gauss constraints (11) and (14), respectively, where χ k (x ⊥ ) are standard Gaussian distributed random variables on the transverse plane. The amplitude of the fluctuations is parametrised by Δ. So far there is no theoretical prediction for its value, which is yet another model parameter we shall vary in order to study its effect on the physical results. Note that, in principle, these modelled fluctuations could be replaced by the spectrum calculated at NLO from the initial conditions, without additional parameters [44]. To our knowledge, this has not been implemented so far, and we first assess the relative importance of fluctuations before attempting such a task.

Setting the lattice scale and size
In a non-equilibrium problem, a scale is introduced by the physical quantity specifying the initial condition. In our case this is the magnitude of the initial colour charge distribution defined in (17) and we follow again [19] in setting the dimensionless combination QL = 120, where L corresponds to the transversal box length in physical units. It is chosen to correspond to the diameter of an Au atom with A = 197, R A = 1.2 A 1/3 fm ≈ 7 fm. In the LHC literature it is conventional to define the transverse section of the box by π R 2 A = L 2 , which then sets the transverse lattice spacing through L = N ⊥ a ⊥ . Together with Q = 2 GeV we thus have As long as we do not add any term describing quantum fluctuations, the system reduces to a 2D problem and thus the results are independent of a ζ . For non-vanishing fluctuations in the static box we work with an isotropic spatial lattice, i.e. a z = a ⊥ , whereas our 3D simulations in comoving coordinates are performed at a η N η = 2.0 as proposed, e.g., in [45].

Observables
Energy density and pressure are convenient observables to investigate the early isotropisation process of the plasma. The system's energy density is the 0th diagonal element of the energy-momentum tensor, T 00 , and can be separated into its evolving chromo-magnetic and chromo-electric components, B and E , respectively, and further into transverse and longitudinal components, On the lattice, the chromo-electric and chromo-magnetic contributions to the Hamiltonian density in Cartesian coordinates, H ≡ T tt , are The contributions to the lattice Hamiltonian density in comoving coordinates, H ≡ τ T τ τ , read Summing the transverse and longitudinal components over the lattice then gives the averaged energy density contributions, with the lattice volume V = N 2 ⊥ N ζ . A suitable measure for isotropisation is given by the ratio of longitudinal and transverse pressure. These are given by the spatial diagonal elements of the energy momentum tensor, Note that at early times the field component of the longitudinal pressure is negative. This is due to the leading order of the CGC initial condition which sets P L to exactly the negative value of P T [46], and reflects the force of the colliding nuclei. In complete equilibrium both pressures are equal.

Validity of the classical approximation
One necessary requirement for a continuum quantum field to behave effectively classically is a high occupation number N of its field modes. In addition, for a classical description to be a good approximation, the IR sector should dominate the total energy of the system, since the classical theory breaks down in the UV. Occupation numbers can be defined unambiguously for free fields, and only then. In the framework of canonical quantisation, and in a fixed gauge, the Fourier modes of the gauge and chromo-electric fields correspond to annihilation and creation operators of field quanta (or gluons) with energy ω(p) = |p|. With proper normalisation, these combine to a number density operator n(p), returning the number of gluons with momenta in (p, p + dp) when acting on an arbitrary Fock state. The (vacuum subtracted) Hamiltonian of the free theory can then be expressed simply by counting the excited field quanta of momentum p, For interacting fields, the interpretation of their Fourier modes is changed and occupation number cannot be defined rigorously. It is thus a valid concept only for sufficiently weak coupling and weak fields. In this case it has been shown that the energy contribution of the gauge-fixed field modes according to the last equation agrees well with the gaugeinvariant energy of the system, see e.g. [49].
In order to study the population of different momentum modes and their contribution to the energy, it is therefore customary to compute the Fourier components of the gauge and/or chromo-electric field, e.g. [47][48][49][50], or, close to equilibrium, those of field correlators [40,51]. However, besides the gauge-dependence, which gets amplified for interacting and strong-fields and causes ambiguities in the interpretation of the momentum distribution, this also introduces a significant computational overhead for the process of gauge fixing, especially for SU (3). For this reason, we turn the procedure around and consider the spectral decomposition of the manifestly gauge-invariant Fourier transform of the total energy density, whose average over equal absolute values of momenta 2 normalised on that momentum, provides a measure for the population of momentum modes. That is, we define an alternative occupation number density with the physical volume V and (p) ≡ H(p) in the static case and (p) ≡ H(p)/τ in the expanding one, respectively. We used the eigenfrequency ω corresponding to the free massless dispersion relation ω( p) ≈ p, as is appropriate for p 0.1 Q [47]. In the non-interacting limit and Coulomb gauge, this definition results in the same energy density as the gauge-fixed ones used in the literature [47][48][49][50]. In the interacting and strong field case, when occupancy becomes ambiguous, our definition removes any gauge dependence while retaining its physical interpretation based on Fourier modes of the energy density.
p ≡ |p| = p 2 x + p 2 y + p 2 z . This is indicated by the notation · p .
Another question is up to which energy level the modes of a classical theory provide a good approximation: because of the Rayleigh-Jeans divergence, the UV sector of the classical theory in equilibrium increasingly deviates from that of the full quantum theory, irrespective of occupation numbers. In thermal equilibrium, a UV cutoff is usually fixed by matching a thermodynamical observable between the full and an effective theory. In a non-equilibrium situation, however, it is difficult to identify a scale up to which the classical theory is valid. A common self-consistent procedure then is to demand that the total energy of the system under study is "dominated by infrared modes".

Ordering of scales and parameters
We wish to study the dependence of the classical Yang-Mills system on the lattice spacing and volume, as well as of the various parameters introduced through the CGC initial conditions. For the classical description of the CGC model to be self-consistent, the parameters representing various scales of the problem have to satisfy The original MV model without additional IR and UV cutoffs corresponds to the special case m = L −1 and Λ = a −1 ⊥ . The dimensionless version of these relations to be satisfied by our lattice simulation is obtained by dividing everything by Q.

Numerical results
Our numerical implementation is based on the well-tested and versatile QDP++ framework [52], which allows for dataparallel programming on high performance clusters. Unless stated differently, we will use QL = 120 throughout this section. Furthermore, as introduced in Sect. 2.2, the initial conditions in the boost invariant scenario, i.e. the one without longitudinal fluctuations, are identical in both frameworks. We will therefore present corresponding results for the energy density solely in the expanding formulation, since the counterparts in the static box can easily be derived therefrom due to energy conservation.

SU(2) vs. SU(3)
Performing the calculations for the realistic SU(3) rather than SU(2) gauge theory introduces roughly an additional factor of 3 in terms of computational time, depending on the studied observables. Comparing physical results between the groups is non-trivial, since the ratio Q s /Q depends on the number of colours, as well as our observables like the energy den- Fig. 1 Total energy density and its chromo-magnetic and chromoelectric components for SU(2) and SU (3) sity. For the saturation scale we have Q s ∼ √ N c Q [37] and for the initial energy density g 2 (t|τ = 0) ∼ N c N g [39]. A physically meaningful, dimensionless combination with the leading N c -behaviour scaled out is thus g 2 /(Q 4 N c N g ) plotted vs.
√ N c Qτ . In Fig. 1, where we applied this rescaling, 3 we clearly see that there is no significant difference in the observables we are studying. In particular, the sub-leading N c -dependence appears to be much weaker than the sensitivity to the parameters of CGC initial conditions, which will be discussed in Sect. 3.4. These results support early findings on the N c -scaling of classical simulations [53]. We checked this observation for several parameter settings with the same outcome and will therefore focus mostly on SU(2) in the following, in order to reduce the numerical cost.

Boundary effects
In the MV model, the nucleus is usually "spread" over the whole lattice. This introduces a systematic error when using periodic boundary conditions. However, for our choice of parameters the total diameter of the plane representing the nucleus is about 12 fm, which should be large enough to suppress boundary effects. In Fig. 2 we show the total energy density (times the proper time τ ) in comoving coordinates for three different scenarios: first, the nucleus is "spread" over the whole 400 2 points on the transverse lattice plane, second, the nucleus is represented by 400 2 lattice points within a 600 2 lattice and third, the same nucleus is embedded in an 800 2 lattice. We observe an effect at the 5%-level. We have explicitly checked that the size of finite volume effects does not change when additional model parameters are introduced, as in the following subsections. 3 In the following, we will keep the scaling factor for the energy density, but we will drop the √ N c normalisation factor in front of Qτ in order to ease the comparison with other works, where this is almost always neglected, too.

Discretisation effects
Ideally, the non-physical scales a σ or a ⊥ entering our calculations because of the lattice discretisation should have no effect on our results. On the other hand, a continuum limit does not exist for a classical theory and one has to investigate which values of the lattice spacing are appropriate and to which extent observables are affected by it.
For our problem at hand, the transverse lattice spacing is set by the number of lattice points spanning the size of the nucleus, cf. (26). On a coarser lattice less momentum modes are available, which translates into lower initial energy density for a fixed colour charge density Q, as shown in Fig. 3  (top). For a non-expanding system the energy density stays constant, thus implying large discretisation effects. In the expanding system, these differences are quickly diminished below percent level, which in the literature is often interpreted as a sign for continuum-like behaviour. However, this behaviour should not be confused with a proper continuum limit, which does not exist for the classical theory. Rather, the expansion adds more and more infrared modes to the system, thus "diluting" the initial UV modes affected by the lattice cutoff and maintaining the apparent classicality. Note also, that the apparent freedom to choose a lattice spacing results from our ignorance of the detailed physics. While yet unknown, there must be a relation (Q) between energy density and colour charge density for given nuclei and collision energy. The lattice spacing would then be fixed by matching the energy density of the classical system to the physical one, similar to the situation in equilibrium.
For our further investigations we will choose a 400 2 lattice, since it is a reasonable compromise between small discretisation effects and computation time. As can be seen in Fig. 3 (top), with this choice the discretisation effects are negligible for Qτ 0.3.
We also have to be sure that there are no discretisation effects coming from the numerical integration over the time variable. To this end we vary the anisotropy parameter ξ , with the results for the transverse and longitudinal energy density shown in Fig. 3 (bottom). We used ξ = 20 ⇔ a t|τ = 0.05 a σ |⊥ for all the results presented in this work, since this choice leads to negligible systematic errors coming from our time discretisation.

Investigation of the parameters of the CGC initial conditions
In the following we elaborate on the different parameters entering the system's description through the CGC initial conditions.

Number of longitudinal sheets N l
As shown in [39], the originally proposed initial conditions of the MV model lack randomness within the longitudinal dimension. Fukushima proposed to use N l sheets of the nucleus rather than only a single one. This is a merely technical parameter coming from the numerical implementation and thus vanishes in continuous time, where N l → ∞. Figure 4 shows that the total energy density depends strongly on N l for small values 30 and then saturates. This effect is amplified by adding an IR cutoff m, leading to a faster saturation for m/Q = 0.1 than for m/Q = 0. This has also been observed in [37] and can be expected: the IR cutoff introduces an additional screening of the colour sources and hence reduces the correlation length also in the rapidity direction. The computation time of the system's initialisa-tion grows linearly with N l and hence a reasonable choice is N l = 30, which we set for most of our simulations.

IR cutoff m
As explained in the last section, the IR parameter m provides a simple way to incorporate the colour neutrality phenomenon studied in [54]. While m = 0.1 Q ≈ 1 R p , with R p being the proton radius, is a physically motivated choice, the precise value of m/Q has a large effect on the initial energy density which can be seen in Fig. 5 (top). With a higher cutoff, less modes are populated to contribute to the energy density. As studied in [37], the parameter m also affects the ratio Q/Q s : at N l = 30 the physical saturation scale Q s is around 0.85 Q for m/Q = 0.1 and around 1.03 Q for m/Q = 0. Since the energy density is normalised by Q 4 , this difference amounts to about a factor of 2 in the dimensionless quantity /Q 4 s . Since the effect of m is in the infrared, it does not get washed out by the expansion of the system, in contrast to the discretisation effects. Hence a careful understanding to fix this parameter is important. For example, one might wonder whether this inverse length scale should not also be anisotropic in the initial geometry. In what follows we will either use m = 0, as in the initial MV model, or the physically motivated choice m/Q = 0.1.

UV cutoff Λ
As discussed in Sect. 2.2, one can apply a UV cutoff Λ while solving Poisson's equation (20a), in addition to the existing lattice UV cutoff. This is an additional model parameter limiting the initial mode population to an infrared sector determined by Λ. Figure 5 (bottom) shows the influence of this parameter on the energy density, which gets reduced because of the missing higher modes in the Poisson equation. This is similar to the observation we made on the IR cutoff m, but with the important difference that the ratio Q/Q s is independent of Λ [55]. We are not aware of a unique argument or procedure to set this parameter, for the sake of comparison with the literature we choose Λ/Q = 1.7 [19] in some of our later investigations. As a welcome side effect, with the emphasis of the infrared modes strengthened, the dependence of the total energy density on the lattice spacing is reduced and the expanding system saturates even faster towards a ⊥independent values, cf. Fig. 6 and the previous Fig. 3 (top).

The energy density mode spectrum
The occupation number of field modes in Fourier space is the most direct and often applied criterion to judge the validity of the classical approximation during the time evolution of the system. It is well-established that, starting from CGC with an additional UV cutoff of Λ = 1.7 Q initial conditions, simulations in a static box quickly populate higher modes, implying a breakdown of the classical description beyond some time. In the expanding system this process is considerably slowed down [19,38,56,57]. We confirm these earlier findings by plotting our generalised occupation number as a function of the momentum modes defined via (35).  full range of the additional UV cutoff, we deliberately chose Λ = Q as its smallest value, cf. (36). One clearly sees that the additional UV cutoff causes a strong suppression of higher modes, thus strengthening the validity of the classical approximation. This is also consistent with the observation from Sect. 3.4.3, that the additional cutoff can be used to weaken discretisation effects. Another observation is that the distribution is rather independent of the IR cutoff value. In Fig. 7 (bottom) we present the evolution of the same initial configuration in the static and expanding framework. While without an additional UV cutoff the distributions nearly reach a plateau in the static box, the occupation of the higher modes in the expanding system stays considerably lower, thus extending the validity of the classical approximation.
One can now try to get a quantitative measure of the supposed dominance of infrared modes. By integrating the Fourier modes of the energy density up to some momentum scale, one can infer the energy fraction of the system contained in the modes below that scale, thus assessing the classicality of the mix (see for example [29]). For example, without applying any cutoffs, integrating modes up to 2Q ≈ 4 GeV contains 65% of the total energy of the system at initial time. At Qt|Qτ = 150, this changes to 60% or 77% in the static and expanding cases, respectively. Hence, the quality of the classical approximation deteriorates only slowly or not at all. Nevertheless, a significant systematic error should be expected when several 10% of the energy is in the UV sector, where a running coupling and other quantum effects should be taken into account. This must certainly be the case when modes 5Q ≈ 10 GeV get significantly populated, as in Fig. 7. At this stage of the evolution a better description might be obtained by an effective kinetic theory [26][27][28], where quantum effects are already included.
Finally, we remark that the Fourier mode distribution of energy density, like occupation number in a free field theory, is also sensitive to the homogeneity of the system in coordinate space: a plane wave with only one momentum mode occupied corresponds to a (finite) delta peak in occupation number, whereas wave packets have broader distributions.

Isotropisation
In this section we add small quantum fluctuations on the initial conditions, as described by Eq. (25). These initial fluctuations lead to an eventual isotropisation of the system, which can be studied by the evolution of the ratio of the pressure components P L /P T . To include their effects, we have to extend our two-dimensional analysis by an additional longitudinal direction N z|η , increasing the computation time linearly with N z|η . Within our computational budget, this forces us to use smaller lattices (200 3 ) for this section, thus inevitably increasing the cutoff and finite volume effects we have discussed so far. However, as we shall see, the effects of the model parameters are by an order of magnitude larger.

Static box
We begin with the static box. The general behaviour of the pressure ratio P L /P T has been known for a while and is shown in Fig. 8. After a peak at around Qt ≈ 0.6 follows an oscillating stage until the system isotropises. The oscillating  Table 1 The initial total energy density and its relative increase due to the fluctuations for different cutoff setups. First row: no additional cutoff, second row: m/Q = 0.1, third row: m/Q = 0.1 and Λ/Q = 1.7. The statistical errors are all below the 1 %-level  [18,19] and precludes a hydrodynamical description. We see a strong finite size effect in N z , Fig. 8 (top), which decreases for larger values and should vanish in the limit N z → ∞. For very small values of N z ≤ 10, the fluctuations cannot evolve and the system behaves as in the unperturbed Δ = 0 case. The dependence on the fluctuation amplitude Δ is studied in Fig. 8 (bottom). In accord with expectation, increasing the fluctuation amplitude Δ reduces the isotropisation time. Note the interesting dynamics associated with this: while for larger initial amplitudes the onset towards isotropisation occurs earlier, the eventual growth of the longitudinal pressure appears to be faster for the smaller amplitudes. The initial fluctuation amplitude Δ also significantly affects the early behaviour of the system, causing a strong change of the pressure ratio and a significant increase of the energy density (∼ Δ 2 ), as shown in Table 1. Also the frequencies of the plasma oscillations are affected. Of course, increasing the quantum fluctuation amplitude weakens the classicality of the initial condition: for Δ ≥ 0.1 the fluctuations already make up ≥ 20% of the initial energy density. On the other hand, for Δ 10 −2 there is no visible effect on the pressure ratio at early times (Qt 20), and also the energy remains the same within numerical fluctuations. The hydrodynamisation time of a heavy ion collision is the time, after which hydrodynamics is applicable to describe the dynamics of the system. This is commonly believed to be the case once the pressure ratio P L /P T ≥ 0.7. For an initial amplitude of Δ = 10 −2 and without further model cutoffs, this happens at t ≈ 770/Q ≈ 76 fm in our simulations. This value is considerably larger than experimentally expected ones, but it is in line with earlier numerical results in a static box, e.g. [19].
The pressure ratio is highly sensitive both to the additional IR and to the UV cutoff introduced in the initial condition, cf. Fig. 9 (top). Especially the UV cutoff changes the qualitative shape of the curve at early times significantly. Furthermore, both cutoffs considerably slow down the process of isotropisation as shown in Table 2. The hydrodynamisa-tion time grows by factors of 2-6 for cutoff values as chosen before. Hence, a better understanding and fixing of those model parameters is mandatory for any quantitative investigation.
Note that our tabulated hydrodynamisation times have been obtained by extrapolation. In principle it would be possible to simulate the late stage of the P L /P T evolution and compare its details to the predicted power law behaviour observed in other studies [18,41]. However, our discussion of the mode distribution in Sect. 3.5 suggests that at such late times a purely classical evolution might no longer be selfconsistent, so we refrained from this computationally very expensive investigation.
In accord with Sect. 3.1, we see no significant change in the isotropisation time when using N c = 3 colours instead of 2, cf. Fig. 9 (bottom). By contrast, the details of the oscillatory behaviour at early times differ. This implies that for the investigation of the properties of collective excitations as in [58], the correct gauge group will eventually be important for quantitative results.

Chromo-Weibel instabilities
It has been suggested that the apparent rapid thermalisation during heavy ion collisions might be caused by chromo-Weibel instabilities [7,8]. Indeed, the final increase of the pressure ratio towards isotropisation, as observed in Fig. 8, may be attributed to such an instability, as we now show. Firstly, our anisotropic initial conditions imply a fluctuating current, which is a necessary ingredient for the occurrence of a Weibel instability. Secondly, an instability causes a rapid population of harder longitudinal modes, which during the evolution in time spreads to others, as suggested by the mode spectrum in Fig. 7. The most striking illustration of the presence of a chromo-Weibel instability is obtained by observing the chromo-electric and chromo-magnetic energy densities in position space, where filaments caused by the instability are clearly visible. Figure   The pattern at Qt = 0.3 (first row of Fig. 10) for Δ = 10 −2 and Δ = 10 −3 represents the initial fluctuations which are independent of the longitudinal direction z. At a later time Qt = 90 the chromo-Weibel instability is visible with filaments that are more pronounced for higher fluctuation seeds. At very late times Qt = 300 the filaments dissolve again. Note how the detailed timing of the growth and decay of the filaments crucially depends on the value of Δ. It is interesting to compare these plots with Fig. 8 (bottom): apparently the dynamical instabilities arise late, after the oscillatory period around the onset to isotropisation.
For consistency, we checked that indeed no filamentation arises in the transverse plane, as expected. This holds for all components of both the chromo-magnetic and for the chromo-electric energy density. Instead, the average values of the energy densities are random with large fluctuations at early stages, which get smoothed during the time evolution.

Expanding system
By contrast, in an expanding system, as realised in heavy ion collisions, the pressure ratio does not appear to isotropise after the oscillatory stage but settles at a small or zero value, The red solid lines indicate the constant energy densities corresponding to Λ ∈ {Q, 2Q, 3Q, 4Q} at QL = 120. The yellow solid line refers to the constant energy density contour obtained for QL = 120 without an additional UV cutoff, which is the choice of the majority of previous studies. The gray horizontal line represents the lattice UV cutoff above which the additional UV cutoff no longer affects the system as shown in Fig. 11. This is in accord with the findings in [14,15,41] and robust under variation of all model parameters. In particular, it also holds for the largest fluctuation seed considered, cf. Fig. 11 (bottom). Correspondingly, in the expanding system no dynamic filamentation takes place either. Only for fluctuation amplitudes 10 −1 filaments are forced right from the beginning, since the initial configuration is equivalent to the one we have shown for the static box scenario. The conclusion is that an expanding gluonic system dominated by classical fields according to the CGC does not appear to isotropise and thermalise. For future work it would now be interesting to check whether adding light quark degrees of freedom helps towards thermalisation, as one might expect.
Note that the non-thermalisation of the expanding classical system is in marked contrast to simulations of an effective kinetic theory, which predict hydrodynamisation times consistent with phenomenological expectations [28,59] (see, however, [23]). Further studies of the systematics of both approaches are necessary to see whether quantum effects and/or the role of the UV sector are the reason for this discrepancy.

Initial condition at fixed energy density
Altogether the numerical results of classical simulations show a large dependence on the various model parameters of the CGC initial condition. This creates a difficult situation, because the initial condition and the early stages of the evolution until freeze-out are so far not accessible experimentally. We now propose a different analysis of the simulation data which should be useful in constraining model parameters such as Λ, m and Δ. In a physical heavy ion collision the initial state is characterised by a colour charge density, an energy density and some effective values of Λ, m and Δ. However, these cannot all be independent, rather we must have = (Q, Λ, . . .), where the detailed relation is fixed by the type of nuclei and their collision energy. We should thus analyse computations with fixed initial energy density L 4 , while varying the model parameters. The outcome of such an investigation for Λ and Q are the contour plots shown in Fig. 12. We consider Qτ = 0.3 as well, since then even without an additional UV cutoff the discretisation effects are negligible for N ⊥ = 400, cf. Fig. 3. In the same figure we also compare the situation with an additional IR cutoff as discussed earlier. Thus, to the extent that the energy density as a function of time can be determined experimentally, it should be possible to establish relations between the parameters Q, Λ and m to further constrain the initial state.
The same consideration can be applied to study the fluctuation amplitude. Figure 13 shows contours of fixed energy density /Q 4 in the (Δ, Λ/Q) plane, where Δ = 0 represents the classical MV initial conditions, i.e., the tree-level CGC description without any quantum fluctuations, and we have chosen QL = 120. Clearly, similar studies can be made for any pairing of the model parameters at any desired time during the evolution and should help in establishing relations between them in order to constrain the initial conditions.

Conclusions
We presented a systematic investigation of the dependence of the energy density and the pressure on the parameters entering the lattice description of classical Yang-Mills theory, starting from the CGC initial conditions. This was done in a static box framework as well as in an expanding geometry and both for N c = 2 and N c = 3 colours.
After the leading N c -dependence is factored out, deviations between the SU(2) and the SU(3) formulation are small and only visible in the details of the evolution during the early turbulent stage. This is not surprising in a classical treatment, since in the language of Feynman diagrams most of the subleading N c -behaviour is contained in loop, i.e. quantum, corrections.
Finite volume effects are related to the treatment of the boundary of the colliding nuclei and their embedding on the lattice. Given sizes of ∼ 10 fm, such effects are at a mild 5%level. Note, however, that this effect is larger than the finite size effects of the same box on the vacuum hadron spectrum, as expected for a many-particle problem.
The choice of the lattice spacing affects the number of modes available in the field theory and thus significantly influences the relation between the initial colour distribution and the total energy of the system. In the static box, all further evolution is naturally affected by this. Since the classical theory has no continuum limit, the lattice spacing would need to be fixed by some matching condition at the initial stage. By contrast, in the expanding system the energy density quickly diminishes and the effect of the lattice spacing is washed out.
A quantitatively much larger and significant role is played by the model parameters of the initial conditions, specifically additional IR and UV cutoffs affecting the distribution of modes and the amplitude of initial quantum fluctuations, whose presence is a necessary condition for isotropisation. For the static box we presented direct evidence for isotropisation to proceed through the emergence of chromo-Weibel instabilities, which are clearly visible as filamentation of the energy density. However, the hydrodynamisation time is unphysically large and gets increased further by additional IR-and UV-cutoffs in the initial condition. Without quantitative knowledge of these parameters, the hydrodynamisation time varies within a factor of five. We suggested a method to study the parameters' influence on the system at constant initial energy densities. This allows to establish relations between different parameter sets that should be useful to constrain their values.
Rather strikingly, no combination of model parameters leads to isotropisation in the expanding classical gluonic system.
International Center for FAIR within the framework of the LOEWE program launched by the State of Hesse. S.Z. acknowledges support by the DFG Collaborative Research Centre SFB 1225 (ISOQUANT).

Data Availability Statement
This manuscript has associated data in a data repository. [Author's comment: This is a theoretical work, no experimental data were used.] Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .