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) and SU(3), 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. The application of hydrodynamic models requires the thermalisation time scale from the initial non-equilibrium stage of the collision to the QGP to be very fast [1,2] compared to the lifetime of the QGP. 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 [3]. The resulting strong gauge field dynamics constitutes the first stage of the evolution. The following second stage is then governed by hydrodynamic equations until the medium becomes too dilute for this long wavelength description. The precise duration of stage I is not yet known for realistic values of the coupling. Hydrodynamic models of stage II constrain it to be around or less than 1 fm [4].
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.
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.
In the following section we summarise the theoretical framework of our approach and give the CGC initial conditions this work is based on. In Section 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 [25].
2 Classical Yang-Mills theory on the lattice

Hamiltonian formulation
Our starting point is the Yang-Mills action in general coordinates, For a treatment on an anisotropic, hypercubic lattice in Minkowski spacetime we employ Wilson's formulation 1 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, 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 (2.9) or (2.12), respectively. For a system out of equilibrium there is no partition function. Rather, initial fields satisfying the Gauss constraint have to be specified 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 [3,26]. 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 [27]. 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 [28] the distribution is taken to be Gaussian, with charge den- 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 [29], with Q s ≈ Q . . = g 2 µ.
For P b − P b 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 [30]) 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 gauge-covariance 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 [29,31]. 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 y-direction 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 [29] 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 [29]. 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 Rp , 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 cutoff Λ is used in the literature [18], while solving Poisson's equation (2.19a). It can be interpreted as an additional model parameter, which restricts the colour sources in Fourier space to modes further in the IR. Again, we shall investigate how results depend on the presence and size of this model 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, For the case of N c = 3 we do this numerically using multidimensional root finding methods of the GSL library [32]. For the case of N c = 2, one can find a closed-form expression and circumvent this procedure, i.e. (2.21) reduces to The remaining field components are U ζ (x) = 1, E a k (x) = 0 and with the index convention introduced in Section 2.1.
To make the initial conditions more realistic, fluctuations can be added on top of this background [14,33], which are supposed to represent quantum corrections to the purely classical fields. They are low momentum modes constructed to satisfy the Gauss constraints (2.10) and (2.13), respectively, where χ k (x ⊥ ) are standard Gaussian distributed random variables on the transverse plane. The amplitude of the fluctuations is parametrised by ∆. Since there is no theoretical prediction for its value, it is yet another model parameter which we shall vary in order to study its effect on the physical results.

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 (2.16) and we follow again [18] 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, 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 [34].

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 [35], and reflects the force of the colliding nuclei. In complete equilibrium both pressures are equal.

Validity of the classical approximation
One requirement for a 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. Thus, in order to study the population of the different momentum modes, it is customary in the literature to compute the Fourier components of the chromo-electric field and their contribution to the energy. However, chromo-electric fields and their Fourier modes are gauge-dependent. Besides the ambiguities this causes in the interpretation of the momentum distribution, it also introduces a significant computational overhead for the process of gauge fixing. For this reason we consider the spectral decomposition for a manifestly gauge-invariant quantity, the 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 the 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. This definition is motivated by considering a gas of free gluons [36][37][38], where the system's energy E (in d dimensions) is related to the gluon mode density n via 34) and the eigenfrequency ω corresponds to the free massless dispersion relation ω(p) ≈ p for With a gauge-invariant definition, the occupation number per mode can be determined fairly uniquely. A much harder 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 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 [39], which allows for data-parallel programming on high performance clusters. Unless stated differently, we will use QL = 120 throughout this section. Furthermore, as introduced in Section 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 density. For the saturation scale we have Q s ∼ √ N c Q [29] and for the initial energy density g 2 (t|τ = 0) ∼ N c N g [31]. 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 Figure 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 Section 3.4. 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 Figure 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.

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. (2.25). 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 Figure 3 (left). 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.
Note however, 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 Figure 3 (left), 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 Figure 3 (right). 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 [31], 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 [29] 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 initialisation 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 [40]. While m = 0.1 Q ≈ 1 Rp , 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 Figure 5 (left). With a higher cutoff, less modes are populated to contribute to the energy density. As studied in [29], 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 Section 2.2, one can apply a UV cutoff Λ while solving Poisson's equation (2.19a), 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 (right) 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 Λ [41]. 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 [18] in some of our later investigations. As a 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. Figure 6 and the previous Figure 3 (left).

The mode spectrum
To our knowledge, the occupation number of the field modes in Fourier space is currently the only criterion applied to judge the validity of the classical approximation during the time evolution of the system. It is well-established that, starting from CGC 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 [18,30,42,43]. We confirm these earlier findings by plotting the occupation number as a function of the momentum modes defined via (2.33). Figure 7 (top) shows the mode spectra for different model parameter values at initial time. In order to study the full range of the additional UV cutoff, we deliberately chose Λ = Q as its smallest value, cf. (2.35). One clearly sees that the additional UV cutoff causes a strong suppression of higher modes, thus strengthening the validity of the classi- cal approximation. Another observation is that the distribution is rather independent of the IR cutoff value. In Figure 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. 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 Figure 7. At this stage of the evolution a better description might be obtained by an effective kinetic theory [44][45][46], where quantum effects are already included.

Isotropisation
In this section we add small quantum fluctuations on the initial conditions, as described by eq. (2.24). 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 Figure 8. After a peak at around Qt ≈ 0.6 follows an oscillating stage until the system isotropises. The oscillating stage originates from turbulent pattern formation and diffusion [17,18] and precludes a hydrodynamical description. We see a strong finite size effect in N z , Figure 8 (left), 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 Figure 8 (right). 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. [18]. The pressure ratio is highly sensitive both to the additional IR and to the UV cutoff introduced in the initial condition, cf. Figure 9 (left). 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 hydrodynamisation 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.
In accord with Section 3.1, we see no significant change in the isotropisation time when using N c = 3 colours instead of 2, cf. Figure 9 (right). 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 [47], 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 [6,7]. Indeed, the final increase of the pressure ratio towards isotropisation, as observed in Figure 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 modes during the evolution in time, which is clearly realised in our system, as shown by the occupation number in Figure 7. The most striking illustration that this indeed corresponds to 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 10 shows the amplitude of the x-component of the chromo-magnetic energy density in the yz-plane while averaging  Figure 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 Figure 8 (right): 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, as shown in Figure 11. This is in accord with the findings in [48] and robust under variation of all model parameters. In particular, it also holds for the largest fluctuation seed considered, cf. Figure 11 (right). Correspondingly, in the expanding system no dynamic filamentation Figure 12: SU(2) total energy density on a 400 2 lattice with N l = 30 as a function of QL and ΛL at initial time and at Qτ = 0.3. The dashed lines represent constant energy density levels at integer multiples of 10 8 . 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. 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.

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 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 Figure 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. Figure 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.