Time evolution of linearized gauge field fluctuations on a real-time lattice

Classical real-time lattice simulations play an important role in understanding non-equilibrium phenomena in gauge theories and are used in particular to model the prethermal evolution of heavy-ion collisions. Due to instabilities, small quantum fluctuations on top of the classical background may significantly affect the dynamics of the system. In this paper we argue for the need for a numerical calculation of a system of classical gauge fields and small linearized fluctuations in a way that keeps the separation between the two manifest. We derive and test an explicit algorithm to solve these equations on the lattice, maintaining gauge invariance and Gauss's law.


I. INTRODUCTION
Particle production at central rapidities in collisions of high energy hadrons or nuclei is dominated by the clouds of small-x gluons surrounding the projectiles. The high density of these gluons has been argued to lead to "gluon saturation", i.e., the emergence of a dominant semihard transverse momentum scale Q s Λ QCD where the physics become nonperturbative due to the nonlinear interactions of the gluons even at weak coupling [1]. The saturation picture of a weak coupling and a nonperturbatively large phase space density of gluons f ∼ 1/g 2 leads to a description of the initial stages of a heavy ion collision in terms of "glasma" fields [2], strong boost invariant color fields with transverse coherence length ∼ Q −1 s . How these maximally anisotropic far-from-equilibrium gauge fields hydrodynamize, isotropize, and reach local thermal equilibrium to form quark-gluon plasma has been a central open question in understanding the spacetime evolution of the matter produced in a heavy-ion collisions.
The large phase space occupancy, or equivalently the strength of the gauge fields, at the early stages of the collision admits a classical description of the glasma fields accurate to leading order in g. The classical description, however, poses a problem phenomenologically as the boost invariance of the fields is not broken and the system remains anisotropic at all times, never thermalizing or reaching hydrodynamical flow.
For the process of isotropization to proceed, it is necessary (but not sufficient) that the boost invariance is bro- * Electronic address: a.k@cern.ch † Electronic address: tuomas.v.v.lappi@jyu.fi ‡ Electronic address: jarkko.t.peuron@student.jyu.fi ken by small rapidity-dependent fluctuations. The origin of the fluctuations may be quantum [3][4][5][6][7] or arise from the longitudinal structure of the colliding nuclei [8,9]. It is then expected that in the presence of the anisotropic background, some of these fluctuations are unstable and experience a period of exponential growth, playing an important role in the isotropization process [10][11][12][13]. Assuming a parametric scale separation between the dominant scale Q s and the inverse wavelength of the unstable modes gf 1/2 Q s , the growth and saturation of the plasma instabilities can be studied in a "hard loop" (HL) framework in which the modes at the scale Q s are treated as quasiparticles and the unstable modes as classical fields. Many calculations have been performed in this framework both analytically [13][14][15][16][17][18][19][20][21] and numerically [22][23][24][25][26]. This is indeed a valid approach when the isotropization process is already under way and the system is only moderately anisotropic and the occupation numbers f of gluonic states with p T ∼ Q s have decreased from their initial value ∼ g −2 . The method however fails at the earliest time scale after the collision, τ ∼ 1/Q s , when the role of the instabilities are expected to be the most important.
The contribution of plasma instabilities to isotropization has also been studied using purely classical field simulations [27][28][29][30] without performing the Hard Loop approximation. These calculations typically proceed using the so called "classical statistical approximation" (CSA). This consists of identifying the initial field fluctuations of the fields, adding these to the classical background field, and then solving the time evolution of the system using the full classical equations of motion on a discrete lattice. Some of these calculations have pointed towards the possibility of a very rapid isotropization caused by the plasma instabilities seeded by the quantum fluctuations of the gauge fields [31].
The treatment of quantum fluctuations in CSA how-ever is problematic due the backreaction of the fluctuations on the background field. Including the quantum fluctuations in the equations of motion of the background is justified only for the modes that grow large and become effectively classical [32,33]; for the other modes, the time evolution of the fluctuations is mistreated. The problem is severe in the case of quantum fluctuations, which have a highly UV-divergent spectrum and occupy modes with f ∼ 1/2 at all scales supported by the lattice. In the CSA these fluctuations are superimposed on top of the background with f ∼ 1/g 2 [31]. Even though the occupancy of the mistreated fluctuations is parametrically smaller than that of the background field, the phase space opens up like 1/a d 3 p, with lattice spacing a. Therefore, on a fine enough lattice the UV tail of the fluctuation spectrum dominates the energy density, particle number, and eventually the dynamics of the system 1 and the time evolution of the combined system cannot be reliably followed in a classical simulation [34]. No continuum limit may be taken (see also [35]).
To avoid this problem, we propose to study the evolution of the fluctuations on a mode-by-mode basis in a setup where the evolution of the fluctuation is explicitly linearized. In this case one can treat the fluctuations to one-loop order, explicitly excluding interactions between the fluctuations and any backreaction to the classical field. One loses the ability to resum late-time "secular divergences" that was one of the motivations for adopting the CSA [4,6]. However, the later-time behavior and eventual hydrodynamization in the context of a heavy ion collision is in any case better described in terms of kinetic theory [36][37][38]. Instead, one keeps the analytical control given by a well defined weak coupling expansion, where different orders in g remain separate. The growth and evolution of the unstable modes can be followed in a clean numerical setup, and one may choose to include only the unstable modes in the simulation. One can also formulate the calculation of gluon production in a dense-dense collision system to NLO accuracy [39,40] analogously to the way quark pair production from the classical field is calculated by solving the Dirac equation in the classical background [41][42][43][44][45].
We will write down the equations of motion for the system of a classical gauge field and linearized fluctuations in Sec. II, noting in particular that maintaining Gauss's law in a calculation with discretized time requires some care. In Sec. III we will present results from simple numerical tests of our algorithm, before pointing in Sec. IV towards some of its potential future applications.

II. EQUATIONS OF MOTION FOR FLUCTUATIONS
In this Section we construct the equations of motion for the linearized fluctuations of the gauge and the chromoelectric field {a i , e i } on top of the background field. On the lattice we will use the Kogut-Susskind Hamiltonian [46] for the background field and in discretizing the equations of motion for the background field we will take special care to make sure that the discretized and linenarized equations of motion exactly conserve the Gauss's law constraint.
In this paper, for simplicity, we will constrain the discussion to a system not undergoing longitudinal expansion (fixed box), however, the extension to a expanding coordinate system is trivial.

A. Small fluctuations in the continuum
In the continuum the Hamiltonian of a pure gauge theory can be written, in temporal gauge A 0 = 0 as with field strength tensor Here we write the gauge and chromoelectric fields in matrix form A i = A a i t a , with the fundamental representation generators t a normalized as Tr t a t b = 1 2 δ ab . From this Hamiltonian one derives the equations of motionȦ In order to project to the physical charge sector, also the Gauss's law constraint must be fulfilled which is conserved exactly by the equations of motion ∂ t C(x, t) = 0. Dividing the field into a background field and linearized fluctuations the equations of motion and Gauss's law for the fluctuations becomė where in terms of an adjoint representation scalar field equation for a i supplemented with a gluon chromomagnetic moment term (see, e.g., [47]) . Similarly, the Gauss's law constraint for the fluctuation reads

B. Discretized equations for background
In order to conserve the gauge symmetry exactly, it is convenient to trade the gauge fields A i belonging to the Lie algebra of the group to link matrices U i which are members of the group. The Kogut-Susskind Hamiltonian [46] in terms of the link matrices reads where the spatial coordinate x takes discrete values on a cartesian lattice x = a(n i , n j , n k ), with integers n i , n j , n k and lattice spacing a. Here the plaquette i,j (x) is written in terms of the link matrices whereî, are unit vectors in the i, j directions; see Fig. 1 for an illustration. 2 The lattice fields are related to continuum quantities by U i (x) ≈ e iagA i (x) and E i lat ≈ agE i cont . The Kogut-Susskind Hamiltonian gives us equations of motion that are discrete in space but continuous in time: where the plaquette in the negative j direction is . Here 2 Note that in the discrete formulation from now on we abandon the summation convention for spatial indices i, j, . . . (but not for color indices).
the notation [] ah denotes the antihermitian traceless part of a matrix: where N c is the number of colors.
In order to perform a practical simulations, also the time direction must be discretized. To guarantee time reversal invariance and second order accuracy in the time step dt, the time direction is commonly discretized with the leaprog algorithm, where the electric fields and the links live on alternate timesteps where we have dropped the explicit position arguments for brevity. It is a straightforward exercise to show that both the link and electric field timesteps (15) and (16) separately conserve the discretized version of Gauss's law constraint Finally, let us recall that under a lattice gauge transformation V (x) (which must be time-independent in order to conserve the temporal gauge condition) the links and electric fields transform as It is easy to see that the Hamiltonian (10) is gauge invariant and the equations of motion (15), (16) and Gauss's law (17) gauge covariant under these transformations.

C. Discretized equations for fluctuations
After these preliminaries, let us move to the lattice equations of motion for the small fluctuations. Naturally, there is a certain freedom in writing down the discretized equations; here, we choose to construct the discretized equations so that they satisfy the following requirements: 1. Reduction to the continuum equations of motion (6), (8) in the limit a → 0, dt → 0.
3. Linearity in a i and e i . 4. An exact conservation of a lattice version of a Gauss's law that reduces to (9) in the limit a → 0, dt → 0 at every time step.
We choose here to start from condition 2 by defining the required gauge transformation properties as those of an adjoint representation scalar field: From these it follows that a i must correspond to a variation of the link matrix U i (x) on the left: In the continuum limit, the fluctuation field on the lattice is related to the continuum equivalent though a lat i = aga cont i . We then choose to discretize the perturbation of the electric field by linearizing the r.h.s. of (16), so that which is easily seen to be gauge covariant. Figure 2 illustrates the ordering of the plaquettes and the field fluctuations in Eq. (23). Here we denote by the fluctuation parallel transported from site x +î to site x by and similarly for the fields parallel transported over two links 3 and so on. The links and gauge field fluctuations a i in (23) are evaluated according to the leapfrog scheme at 3 Note that, in our notation there are two identical ways of writing the most complicated terms involving parallel transports over two links time t + dt/2. We emphasize that the choice (23) is not unique, but one could add terms proportional to (dt) 2 or higher powers.
Similarly to the timestep of E i , Gauss's law (17) is linear in the chromoelectric field. The natural choice is then to derive Gauss's law for the fluctuations by replacing E i with E i + e i , U i (x) with U i (x) + ia i (x)U i (x), and taking the linear terms in the fluctuation fields. This yields We now have equations for the timestep of the electric field fluctuation e i and Gauss's law for the fluctuations.
To complete the set of equations, we need also to specify the timestep for a i (x). The first guess would be a straightforward discretization of the continuumȧ i = e i . However, this naive discretization is inadmissible, since it does not conserve the linearized Gauss's law (27). Physically this would mean an unphysical creation of "charges" in the lattice. This can be traced to the fact that Gauss's law involves a covariant derivative using links U i that advance in time simultaneously as their fluctuations a i , and the timestep must reflect this change. Another hint of the subtlety of the step for a i is to see that a linearization of the timestep for the link U i (x) in Eq. (15) would involve developing the exponential e i(E i +e i )dt to linear order in e i , which is a rather complicated expression when E i and e i do not commute. We may, however, construct a valid update for the gauge fields by demanding the Gauss's law constraint to be conserved, It is straightforward to see that this condition holds if the update satisfies 4 where we use a shorthand for the "timelike plaquette" 0i = e iE i dt . Imposing this condition on the gauge field update leads by construction to a time step that conserves Gauss's law.
It is convenient here to separate the parts of a i , e i that are parallel and perpendicular to E i in color space. Denoting Eq. (29) can be solved for a ⊥ i (t + dt) in terms of a ⊥ i (t) and e i⊥ giving the equation of motion for the perpendicular component. Because of the commutator, Eq. (29) gives no condition for the parallel component, and we may complete the equations of motion with the naive discretization that already satisfies Eq. (29). For a practical algorithm, it still remains a technical problem to solve the perpendicular components of the gauge field fluctuations from Eq. (29). For a general gauge group we can write the solution more compactly in the adjoint representation: in terms of the unitary matrix 0i ab = 2 Tr t a 0i t b † 0i , the hermitian matrix c and the N 2 c − 1 component vectors a i and e i with components (a i ) a and (e i ) a . 4 We drop the explicit time argument for the electric field from now on; this will always be dictated by the leapfrog scheme.
In this notation Eq. (29) becomes The parallel components are the null space of the matrix E i , and in this subspace the timelike plaquette acts like the identity: † 0i f = f . Thus parallel components of a i (t) and e i only generate parallel components of a i (t + dt). In the perpendicular color directions, on the other hand, the matrix E i is invertible, and we can write the gauge field timestep as where the notation () −1 ⊥ means a projection to the subspace where the matrix E i is invertible followed by an inversion in that subspace. This equation is our general result for the timestep of the gauge field fluctuation.
In the small dt limit 0i ≈ 1 + i E i dt and we see that Eq. (34) reduces to a i (t + dt) = e i dt + a i (t) as desired. It may seem like a disproportionate amount of trouble to formulate the equation in this way, when the result reduces to the naive discretization in the limit dt → 0 which one wants to take in the end. However, we have found that in practical computations it is essential for a good precision to conserve Gauss's law also in discrete time and not only in the continuous time limit. At this point it is also straightforward to check that the equation is time reversal invariant, ensuring second order accuracy in dt.
Note that the form (34) of the timestep results from a choice made in writing the timestep for e i and Gauss's law in the form Eqs. (23), (27). We could have resolved the ambiguity in linearizing the fluctuations of the timelike plaquette in another way by defining a different electric field fluctuation e.g. by This would make the timestep for a i simpler, but the timestep and Gauss's law for e i mod would have a more complicated form, with the appearence of terms proportional to E i dt.
The general result (34) requires the solution of a system of N c 2 − 1 linear equations. For the special case of SU(2) we can invert the matrix E i analytically using the fact that in the absence of symmetric structure constants the Fierz identity for f abc f ade is particularly simple ijk ilm = δ jl δ km − δ jm δ kl . Thus if, for the perpendicular part, E i a a ⊥ i,a = 0, we have  3: Test of the decomposition of the field in the background field and fluctuation after some finite time. All runs have been evolved to same physical time, which is dtN t = 2 here with dt = 0.01 and N t = 200. The upper points correspond to δĖ and the lower δ A (see Eqs (38) and (39)). The straight lines are fits of the form ax 4 , showing that the observables decrease with the correct power law. and we can write (34) as (36) or in the fundamental representation as (37) We stress that these final versions (36) and (37) are valid for SU(2) only, and e.g. for SU(3) one must use Eq. (34).

III. NUMERICAL TESTS
We now have the equations of motion for the linearized fluctuation: the timestep for e i (Eq. (23)), for a i (general equation in Eq. (34) and SU(2)-specific ones in Eqs. (36) and (37)) and Gauss's law (27). We present here some simple test results from an implementaion of these equations for the SU(2) gauge group.
We construct initial conditions for a background field configuration by setting the gauge fields A i to random values uniformly distributed in the interval [0, 0.9]. The electric fields are set to zero initially in order to satisfy the Gauss's law C(x, 0) = 0. We then construct the link  (41) and (42).
matrices by exponentiating the gauge fields U i = e iagA i . We similarly construct initial fluctuation fields. We then choose a small parameter ε ranging from 0.5 to 0.0001 and multiply the fluctuations by ε, effectively setting ε as the scale of these fluctuations, i.e. e i , a i ∼ ε. We can now evolve separately in time: 1. The system of the background field and linearized fluctuations E i , A i , e i , a i and 2. A different pure background field configuration initialized asÊ i (t = 0) = E i (t = 0) + e i (t = 0) and A i (t = 0) = A i (t = 0) + a i (t = 0).
If we have now successfully linearized the classical equations of motion, the squared differences and should scale as ε 4 with the magnitude of the fluctuation. For numerical convenience it is easier for us to plot the corresponding differences for the time derivatives E i (t + dt) − E i (t) etc. as the expression involving time derivatives is easily obtained during the time-evolution. In Fig. 3 we show that indeed these differences scale in the correct way for a large range of ε. We can note here that for a naive a i timestep a i (t + dt) = a i (t) + e i the correct ε-scaling for small fluctuations is only obtained for prohibitively expensive small values of dt because the correct scaling of δ E and δ A is violated by terms of the order ε 2 dt 4 . We have verified this numerically by observing that for larger dt δ E and δ A begin to scale as ε 2 . This means that the naive timestep for a i does not correctly capture the differenceÂ i − A i even to leading order in ε.
We also show, in Fig. 4, the violation of Gauss's law constraint as a function of time. To quantify the violation we consider for the background field and similarly for the fluctuations These two quantities measure how well the components in different spatial directions cancel each other. Due to numerical roundoff error, Gauss's law is never satisfied exactly. However, our algorithm preserves it equally well for the fluctuations as for the background field. Also, the fact that Gauss's law remains satisfied orders of magnitude more precisely in double than single precision shows that that the remaining values are purely due to the limited machine precision.

IV. CONCLUSIONS AND OUTLOOK
Unstable fluctuations seeded by quantum effects around a boost invariant classical background field play an important part in the pre-equilibrium evolution of heavy-ion collisions. Until now, the Classical Statistical Approximation has been common tool to study these phenomena. However, the very UV dominated spectrum of vacuum fluctuations in field theory makes attaining the continuum limit in CSA calculations very difficult if not impossible.
We have argued in this paper that it would be desirable to address these issues by real time lattice calculations with an explicitly linearized fluctuation around the classical field. We have here explicitly derived and tested equations of motion for these fluctuations, showing that satisfying Gauss's law for the fluctuations requires a careful treatment in the discretization of the timestep. By giving up the attempt to resum asymptotical long time "secular divergences," which are not a problem with a matching to kinetic theory, one stands to gain better control of the UV dynamics in the classical gauge field calculation. We expect this formalism to have several interesting applications, which we plan to return to in future work.