Conserved energy–momentum tensor for real-time lattice simulations

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Conserved energy–momentum tensor for real-time lattice simulations © The Author(s) 2024 Published version Boguslavski, K.; Lappi, T.; Peuron, J.; Singh, P. Boguslavski, K., Lappi, T., Peuron, J., & Singh, P. (2024). Conserved energy–momentum tensor for real-time lattice simulations. European Physical Journal C, 84(4), Article 368. https://doi.org/10.1140/epjc/s10052-024-12725-6 2024 Eur. Phys. J. C (2024) 84:368 https://doi.org/10.1140/epjc/s10052-024-12725-6 Regular Article - Theoretical Physics Conserved energy–momentum tensor for real-time lattice simulations K. Boguslavski1, T. Lappi2,3, J. Peuron2,3,P.Singh 2,3,a 1Institute for Theoretical Physics, Technische Universität Wien, 1040 Vienna, Austria 2Department of Physics, University of Jyväskylä, P.O. Box 35, 40014 Jyväskylä, Finland 3Helsinki Institute of Physics, University of Helsinki, P.O. Box 64, 00014 Helsinki, Finland Received: 17 January 2024 / Accepted: 26 March 2024 © The Author(s) 2024 Abstract We derive an expression for the energy–momentum tensor in the discrete lattice formulation of pure glue QCD. The resulting expression satisfies the continuity equation for energy conservation up to numerical errors with a symmetric procedure for the time discretization. In the case of the momentum conservation equation, we obtain an expression that is of higher accuracy in lattice spacing (O(a2)) than the naive discretization where fields in the continuum expressions are replaced by discretized counterparts. The improvements are verified by performing numerical tests on the derived expressions using classical real-time lattice gauge theory simulations. We demonstrate substantial reductions in relative error of one to several orders of magnitude compared to a naive discretization for both energy and momentum conservation equations. We expect our formulation to have applications in the area of pre-equilibrium dynamics in ultrarelativistic heavy ion collisions, in particular for the extraction of transport coefficients such as shear viscosity. 1 Introduction Lattice gauge theory [1] is a powerful tool commonly used to address nonperturbative problems in Quantum Chromodynamics (QCD). This lattice framework encodes gauge invariance by construction and therefore preserves this crucial property of gauge theory. However, the lattice discretization modifies or breaks some of the underlying symmetries of the theory, such as translational or rotational invariance. The energy–momentum tensor (EMT) is an observable that encodes energy and momentum conservation in the form of a continuity equation. The conserved quantities in this case arise from invariance under space-time translations as ae-mail: pragya.phy[email protected] (corresponding author) described by Noether’s theorem. Furthermore, the energy– momentum tensor is traceless at the classical level due to conformal symmetry. In a lattice discretized formulation, these symmetries are modified, as invariance under spacetime translations becomes a discrete symmetry instead of a continuous one. The conformality of the theory is also broken by involving an explicit scale, the lattice spacing. In spite of these issues, one typically evaluates the classical energy– momentum tensor by taking its continuum expression and replacing the fields with their discretized counterparts. It is not clear a priori that the resulting expression satisfies the continuity equation. Our interest in this topic is motivated by the early-time dynamics in the context of ultrarelativistic heavy-ion collisions. Directly after the collision the system consists of overoccupied gluon fields, is typically modeled using classical Yang–Mills theory, and is often referred to as the Glasma [2–4]. Recently, the Glasma stage has received substantial attention [5–22], especially in terms of its properties as a medium, described by transport coefficients [23–30]. This paper aims to construct an improved expression for the lattice energy–momentum tensor of a nonabelian gauge field following classical equations of motion, such that the tensor satisfies the continuity equation exactly whenever possible and improves the accuracy of energy–momentum conservation when this is not the case. As our applications lie in the domain of real-time simulations [4–6,10–12,15– 20,23–25,31,32], where a classical-statistical approximation is employed, our approach will inevitably differ from that used in lattice QCD. This is due to the fact that the EMT in lattice QCD is constructed to satisfy Ward–Takahashi identities [33–35] at the quantum level after renormalization. These identities play the role of Noether’s theorem for quantum field theories. Hence, while our approach is similar aiming to cap- 0123456789().: V,-vol 123 368 Page 2 of 18 Eur. Phys. J. C (2024) 84:368 ture the same physical phenomena, the two approaches are nevertheless not identical. For an impatient reader, who is only interested in our main results and less in the technical details concerning their derivation, we summarize our results as follows. The energy density component T00 of the EMT is given by Eq. (18), which redistributes the total energy on the lattice into local energy densities that are spatially symmetric and centered at a lattice site. The Poynting vector components T0iare then obtained by requiring that a discretized continuity equation for energy conservation be satisfied. This procedure leads to T0igiven by Eq. (29) that satisfies the energy conservation continuity equation up to numerical errors. We will also discuss how to further improve this expression by making use of time symmetrization to synchronize the electric and magnetic fields in Sect. 3.2.2. For the spatial components Tij,westart from the Poynting vector and impose the continuity equation for momentum conservation. This leads to an expression for the discrete analog of the Maxwell stress tensor Tij.However, within this derivation, we will have to perform a few approximations due to additional complications. These arise from parallel transport on the lattice and from the fact that on the lattice we do not have a suitable analog of the continuum Bianchi identity. Our final result is given by Eq. (63) and, as we will show in a separate section, it satisfies the continuity equation to an accuracy of order O(a2 s). We believe that the idea and the algorithm presented here to derive a conserved Tμν (with high precision) can be utilized much more broadly. Nonconservation of Tμν is expected to occur in every theory discretized on a lattice due to the breaking of the space-time translational symmetry. Particularly interesting examples of potential applications beyond heavyion collisions are real-time lattice simulations performed in the context of postinflatory cosmology using the classical approximation [36–38]. This paper is structured as follows. Section 2describes our discretization framework. In Sect. 3we construct the discrete energy–momentum tensor and introduce its timesymmetrized formulation. Section 4shows the numerical comparison among the naive and discretized approaches. We conclude in Sect. 5. 2 General formalism and classical equations of motion In line with established practice within classical-statistical lattice gauge theory simulations, we begin with the discrete Kogut–Susskind Hamiltonian [39] of classical SU(Nc) Yang–Mills theory in temporal axial gauge (At=0), H= m a3 s−1 2−det(gμν) μ>0 gμμ Eμ,a mEμ,a m +2Nc g2 0<μ<ν 1 a2 μa2 ν gμμgνν ×−det(gμν)1−1 Nc ReTr(Uμν,m)(1) = m a3 s1 2 μ>0 Eμ,a mEμ,a m +2Nc g2 0<μ<ν 1 a2 μa2 ν1−1 Nc ReTr(Uμν,m).(2) Here mdenotes the lattice site on the spatial grid and a=1,...,N2 c−1 is the color index. In this paper, we consider the Minkowski metric with gμν =(+1,−1,−1,−1) and −det(gμν)=1. The coupling constant genters in the form of the gauge coupling 2Nc/g2. The discretization is performed on a three-dimensional lattice with size Nx×Ny×Nz, lattice spacing aμin the spatial direction ˆμ, and the product a3 s=axayaz. To guarantee gauge invariance, the theory is formulated in terms of gauge links Ui,m=eiaigAi,minstead of gauge fields Ai,m. Links enter Eq. (2) in terms of plaquettes Uμν,m=Uμ,mUν,m+ˆμU† μ,m+ˆνU† ν,m,(3) where ˆμdenotes a unit vector in the μdirection.1In analogy to the links, plaquettes are related to the field-strength tensor Fμν via Uμν ≈eigaμaνFμν . Links and plaquettes are group elements Uμ,m,Uμν,m∈ SU(Nc) in the fundamental representation, and correspondingly we will use the generators of the fundamental representation tato go between the color components of the electric field and the fundamental representation matrix Eμ m=Eμ,a mta. Using the (continuum) relation between the field-strength tensor and the chromomagnetic field Bi= −1 2ijkFjk, the Hamiltonian (2) reduces to the correct expression d3x1 2j>0(Ej,aEj,a+Bj,aBj,a)in the continuum limit aμ→0. The Hamiltonian (2) is associated with2the Hamiltonian equations of motion for the link matrices and the electric field variables, which are discrete in space but continuous in time ∂0Uμ,m=igaμEμ,a mtaUμ,m(4) ∂0gEμ,a m= 0<ν=μ 2aμ a2 μa2 ν ReTritaUμν,m+itaUμ−ν,m(5) = 0<ν=μ 2 aμaν ReTrita[DB ν,Uμν,m].(6) 1To simplify the notation, we will often write μinstead of ˆμ. 2We have rescaled the chromoelectric fields Eμ,a mto correspond to the continuum field, with dimension GeV2, thus they differ from the canonical momentum variables of the discrete theory by a multiplicative factor. 123 Eur. Phys. J. C (2024) 84:368 Page 3 of 18 368 These equations of motion are typically solved using a leapfrog algorithm, which is second-order accurate in the time step dt. In this discretization paradigm, the electric fields and links are located half a timestep apart. This is illustrated in the right panel of Fig. 1. In order to derive Eq. (6) from Eq. (5), we have used the identity U† μ−ν,m=U† ν,m−ˆνUμν,m−ˆνUν,m−ˆνand introduced the backward and forward gauge covariant derivatives DB μ,mX=Xm−U† μ,m−ˆμXm−ˆμUμ,m−ˆμ/aμ(7) DF μ,mX=Uμ,mXm+ˆμU† μ,m−Xm/aμ.(8) By utilizing the equation of motion for the link matrices Eq. (4), one can deduce the time derivative of the magnetic field part of the energy density ∂0ReTr Uμν,m=ReTrigaμaνDF μ,Eν mUμν,m −igaμaνDF ν,Eμ mUμν,m.(9) 3 Constructing the energy–momentum tensor In this section, we outline our methodology for acquiring a discretized representation of the energy–momentum tensor (EMT). In a broad sense, Noether’s theorem introduces the framework for establishing a connection between the temporal and spatial translation invariance inherent in Yang– Mills theory and the energy–momentum tensor. In the continuum, the canonical EMT can then be computed from Noether’s theorem and leads to a non-symmetric tensor that requires an additional gauge transformation to become symmetric. Instead, one can use the (Hilbert) EMT Tμν obtained by taking a functional derivative of the Yang–Mills action S=ma3 s−det(gμν)Lwith respect to the metric tensor gμν, Tμν =2 −det(gμν) δS δgμν =2∂L ∂gμν −gμνL.(10) On the lattice, our starting point is a situation where the metric is purely diagonal. Indeed, both the Hamiltonian (1) and the corresponding action Sonly include diagonal metric elements. Hence, there is no straightforward way to vary off-diagonal components of the metric in a continuous way around zero in order to calculate derivatives with respect to the metric as in (10). This approach is therefore not applicable to our purposes. Furthermore, off-diagonal terms play a crucial role when investigating transport coefficients, such as shear viscosity. We will first describe in Sect. 3.1 how the EMT could be constructed naively inspired by the continuum expressions. Then we will explain our improved procedure that is directly based on the continuity equation for energy–momentum conservation ∂μTμν =0.(11) We will start in Sect. 3.2 by allowing spatial redistribution for the energy density, i.e., the temporal component T00 while requiring that it sums up to the Hamiltonian that corresponds to the total energy (2). The components T0iare then constructed from (11)forν=0. The remaining components will then be constructed in Sect. 3.3 in the same way by taking T0ias the starting point and requiring that Tij satisfy the remaining continuity equations for ν>0. This procedure determines T0iup to a transformation T0i→T0i+φi, where ∂iφi=0. This transformation, however, leaves conserved quantities intact. 3.1 Energy momentum tensor in the continuum and naive discretization The naive discretization of the energy–momentum tensor proceeds by taking the continuum expression and replacing Eand Bfields with their spatially averaged discretized counterparts T00 m=1 2E2 m,loc +B2 m,loc)(12) T0i m=ijkEj m,loc Bk m,loc(13) Tij m=1 2δijE2 m,loc +B2 m,loc −Ei m,loc Ej m,loc −Bi m,loc Bj m,loc =δijT00 m−Ei m,loc Ej m,loc −Bi m,loc Bj m,loc (14) where E2 m,loc =Ei,a m,loc Ei,a m,loc and Ei,a m,loc and Bi,a m,loc are local electric and magnetic (cloverleaf) fields. Since the electric field labeled as Ei,a mcorresponds to the time derivative of the gauge field between mand m+ˆ i, it is rather centered at the point m+ˆ i/2. Similarly, a plaquette Uij,mis actually centered at m+ˆ i/2+ˆ j/2. Thus a natural way to construct electric and magnetic fields at a site mis to take nearest neighbor averages to obtain symmetric expressions: Ei,a m,loc =1 2Ei,a m+U† i,m−ˆ iEi,a m−ˆ iUi,m−ˆ i(15) Bi,a m,loc =− ijk 8gajak ReTr ×itaUjk,m+Uj−k,m+U−jk,m+U−j−k,m (16) 123 368 Page 4 of 18 Eur. Phys. J. C (2024) 84:368 3.2 Continuity equation for energy conservation Our starting point will be the scalar component of the continuity equation ∂0T00 =∂iT0i,(17) which will be used to obtain T0iafter T00 has been constructed. As the Hamiltonian density characterizes the energy density of the system, we utilize this insight to identify the μ=0, ν=0 component of the energy–momentum tensor Tμν . It is worth noting that in the Hamiltonian formulation, the electric field labeled Ei mis located at the position m+ˆ i 2+dt 2, taking also the finite timestep dt into account. The magnetic field strength, expressed by the plaquette Uij,m,is located at m+ˆ i 2+ˆ j 2. To determine the energy density at lattice site m, we compute an average over the “outgoing” and “incoming” electric fields Ej,a mand Ej,a m−jand various plaquette orientations. The left panel of Fig. 1illustrates the spatial averaging procedure for electric and magnetic contributions. This averaging procedure allows us to obtain a representative value for the energy density at a specific lattice site T00,m= j,k>0 1 g2a2 ja2 k ×Nc−1 4ReTrUjk,m+ReTrUj−k,m +ReTrU−jk,m+ReTrU−j−k,m +1 4 j>0Ej,a mEj,a m+Ej,a m−jEj,a m−j(18) Both the electric and magnetic field components of Eq. (18) lead to the same total energy as the Hamiltonian3 (2), while slightly redistributing the local energy density. This formulation preserves the interpretation of the temporal component of the energy–momentum tensor as energy density. In several works of the literature [15,18,21,22,24,25,40], the square of the symmetrized electric field Eq. (15) has been used for the energy density. However, this is not equivalent to the actual Hamiltonian (2): 3In Eq. (18), we employed a summation over all values of jand k instead of the constrained sum in the Hamiltonian (2). The equivalence between the two equations is established by the relation j,k>0= 20<j<k.  j,aEj,a m+U† j,m−jEj,a m−jUj,m−j2 = 2 j,aEj,a mEj,a m+Ej,a m−jEj,a m−j.(19) Furthermore, if one used the square of a symmetrized electric field as on the left-hand side of Eq. (19) to calculate the time derivative of T00, the equation of motion (4) would introduce cubic terms in the electric fields into T0i, which are not present in the continuum expression. Thus using the right-hand side of (19) is more suitable for the purposes of this paper. To write the magnetic field part of Eq. (18) in a form where the equivalence to the Hamiltonian (2) is more explicit, one needs to write the plaquettes that start from the base point m in Eq. (18) as plaquettes with a base point at the “lower left” corner, parallel transported to the site m. This can be done by making use of the identities U−kj,m=(Uj−k,m)†=U† k,m−kUjk,m−kUk,m−k(20) Uk−j,m=(U−jk,m)†=U† j,m−jUjk,m−jUj,m−j(21) U−j−k,m=U† j,m−jU† k,m−j−kUjk,m−j−kUk,m−j−kUj,m−j (22) and the fact that the parallel transporting links cancel in the trace Tr(U†MU)=TrM. This allows us to rewrite the energy density T00,mas T00,m= j,k>0 1 g2a2 ja2 kNc−1 4 ×ReTrUjk,m+Ujk,m−j+Ujk,m−k+Ujk,m−j−k +1 4 j>0 Ej,a mEj,a m+Ej,a m−jEj,a m−j.(23) 3.2.1 Constructing the Poynting vector We now want to use the fact that the continuity equation relates the time derivative of the energy density to the spatial derivative of the momentum density to deduce the components T0i. To employ this method, we apply the evolution equation for plaquettes (9) and electric fields (6) to calculate the time derivative of the energy density in Eq. (23): ∂0T00,m= j,k,j=k ajak 2ga2 ja2 kReTri[DF k,Ej m]Ujk,m +i[DF k,Ej m−j]Ujk,m−j+i[DF k,Ej m−k]Ujk,m−k +i[DF k,Ej m−j−k]Ujk,m−j−k+2ReTriEj m 123 Eur. Phys. J. C (2024) 84:368 Page 5 of 18 368 Fig. 1 (Left:) Illustration of discretized T00,mgiven by Eq. (18). The black lines with blue central points depict the plaquettes, while red circles represent electric fields. Additionally, the magnetic contribution is illustrated by green arrows indicating the orientation of the plaquettes. The blue and red points indicate the position of the plaquettes and electric fields respectively. The separate highlighted standard plaquette is shown in blue, denoted Uij,nbut in fact centered at n+i/2+j/2. (Right:) Presentation of the standard leapfrog algorithm used to update electric and magnetic fields: The electric field at t−dt/2 is used to update the magnetic field from t−dt to t, and similarly the magnetic field (links) at t−dt to update the electric field from t−3dt/2to t−dt/2 ×[DB k,Ujk,m]+iEj m−j[DB k,Ujk,m−j](24) = j,k,j=k aj 2ga2 ja2 kReTr ×iUk,mEj m+kU† k,mUjk,m+iEj mUjk,m +iUk,m−jEj m−j+kU† k,m−jUjk,m−j +iEj m−jUjk,m−j −iUk,m−kEj mU† k,m−kUjk,m−k −iEj m−kUjk,m−k −iUk,m−j−kEj m−jU† k,m−j−kUjk,m−j−k −iEj m−j−kUjk,m−j−k(25) Indeed, it is possible to write the right-hand side as a total spatial derivative. After some rearrangement, one achieves the following result: ∂0T(N) 00,m= j,k,k= j 1 2gajak ×ReTrDB k,iUk,mEj m+kU† k,mUjk,m+iEj mUjk,m ×+iUk,m−jEj m−j+kU† k,m−jUjk,m−j+iEj m−jUjk,m−j. (26) Using the definition of the covariant backward (or forward) derivative in Eqs. (7), (8), it is now easy to see that under the ReTr operation, this expression simplifies to an ordinary derivative of a scalar quantity ReTrDB μ,mX=1 aμ ReTrXm−U† μ,m−ˆμXm−ˆμUμ,m−ˆμ =1 aμ ReTrXm−Xm−ˆμ=∂B μ,mReTrX. (27) Thus we have arrived at the first important result of this paper, a discrete energy conservation law ∂0T00,m=∂B k,mT0k,m.(28) We emphasize that this equation is an exact relation even in the discrete case, although its correspondence with the continuum version is only realized in the limit of small lattice spacing. From Eq. (26) we can read off the momentum density along the k-th component as follows: T0k,m= j>0,j=k cjk (29) ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ReTriEj m+kU−kj,m+k+iEj mUjk,m   T1 0k +iEj m−j+kU−kj,m−j+k+iEj m−jUjk,m−j   same as T1 0kwith m→m−j⎫ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎭ 123 368 Page 6 of 18 Eur. Phys. J. C (2024) 84:368 where cjk =1 2gajak (30) It is important to note that what we label by T0k,mis actually centered at the position m+ˆ k 2,4which may initially seem unconventional. However, this arrangement finds justification in the discrete continuity equation (28), which has a backward derivative in the k-direction. Thus the derivative ∂B k,mT0k, is in fact centered around m, the same position where ∂0T00 is defined. We also note that Eq. (29) corresponds to a gauge covariant formulation of the Poynting vector in the continuum limit T0k→kij EiBj. This can be seen by expanding the plaquettes to quadratic order in as, thus Ujk ≈1+igajakFjk, and realizing that the term with the identity vanishes because the matrix Ejis traceless. 3.2.2 Symmetric time discretization Until now, we have treated tas a continuous variable, which is the limit dt/as1 of the numerical calculation. In practice, however, one wants to choose a larger timestep for numerical efficiency. In fact, we can make a further improvement to remove some of the finite timestep errors from the conservation law in the following way. To reach this goal we can further analyze the time dependence of ∂0T00 as described in Eq. (18). While the Hamiltonian approach we use considers time to be continuous, our numerical simulations do, in fact, have a discrete time step denoted by dt. Let us take a closer look at the time dependence of each of the factors and terms in the equation. To simplify the analysis, we will solely focus on the temporal derivative of one of the electric field components in the T00 expression ∂0TE 00,m=1 2E2 mt−dt 2−E2 mt−3dt 2 dt  =Ei mt−dt 2+Ei mt−3dt 2 2 ×Ei mt−dt 2−Ei mt−3dt 2 dt  =Ei,avg mt−dt∂0Ei mt−dt 2,(31) In the leapfrog scheme, the time difference labeled ∂0Ei mt− dt 2above, corresponding to stepping the electric fields from 4This can be seen by realizing that the first two terms correspond to the continuum gauge fields in the combination Ej(m+k+j/2)+Ej(m+j/2)Fjk(m+j/2+k/2), which is centered around m+j/2+k/2, and the second two terms that are shifted as m→m−jmake the full expression centered around m+k/2. t−3dt/2tot−dt/2, involves the plaquette at the time step t−dt, as illustrated in the right panel of Fig. 1. On the other side of the continuity equation, this term corresponds to the one with a spatial derivative of the plaquette, multiplied by the electric field. Equation (31) tells us that this term in the spatial derivative of T0kshould be evaluated with the following timestep assignment in terms of a time-symmetrized electric field: Ej m[DB k,Ujk,m]→1 2Ej mt−dt 2 +Ej mt−3dt 2DB k,Ujk,m(t−dt).(32) Similarly, the time derivative of the magnetic field part of the energy density corresponds to the term in ∂kT0kwhere one takes a spatial derivative of the electric field. A timesymmetric treatment of this term requires a symmetrization of the links in the evaluation of this term in ∂kT0kas [DF k,Ej m]Ujk,m→1 2DF k(t), Ej mt−dt 2 +DF k(t−dt), Ej mt−dt 2 ×1 2Ujk,m(t)+Ujk,m(t−dt),(33) where the notation DF k(t)refers to the link matrices in the covariant derivative being evaluated at the time t. Note that in the temporal gauge, such time-averagings are gauge invariant, since parallel transporters in the time direction are just identity matrices. We refer to the combination of the time averagings (32) and (33) as the time-symmetrized discrete formulation. 3.3 Continuity equation for momentum conservation Our goal is to construct Tjk in the same way as we constructed T0iabove, i.e., by utilizing the equations of motion and the continuity equation for i>0 ∂μTμi=0(34) to derive the Tij components from ∂0T0i. We first observe that T0kcomprises two kinds of terms, with one of them being merely shifted in the −ˆ jdirection. Henceforth, we will focus solely on the first one of the two, which we call T1 0kfor our subsequent calculations. The second term can then be restored in the end by shifting the site of T1 0k. Let us first take the time derivative of Eq. (29) and split it into terms where the time derivatives act either on the electric fields or the plaquettes. 123 Eur. Phys. J. C (2024) 84:368 Page 7 of 18 368 ∂0T1 0k,m=cjk Ej,a m+k∂0ReTr(itaU−kj,m+k)+Ej,a m∂0ReTr(itaUjk,m)   ∂0T1E 0k +(∂0Ej,a m)ReTr(itaUjk,m)+(∂0Ej,a m+k)ReTr(itaU−kj,m+k)   ∂0T1B 0k . (35) Since the time derivative of a plaquette (9) involves an electric field, the first term will end up being quadratic in E, and reduce in the continuum limit to the electric field part of Tij, which we refer to here as ∂0T1E 0k. Conversely, the time derivative of the electric field (6) is a (discrete) spatial derivative of plaquettes, and the second term will end up being quadratic in the magnetic field in the continuum limit, thus denoted by ∂0T1B 0k.5 3.3.1 Chromoelectric field contribution Let us start with the chromoelectric field contribution ∂0T1E 0k. We begin by employing the equation of motion for the plaquette to express the electric field part ∂0T1E 0kas ∂0T1E 0k,m= j=k gcjkReTr ×aj−Ej m+kU† k,mEj mUk,mU−kj,m+k +Ej m+kEj m+kU−kj,m+k +akEj m+kU† k,mEk mUk,mU−kj,m+k −Ej m+kU† k,mUj,mEk m+jU† j,mUk,mU−kj,m+k +aj−Ej mEj mUjk,m+Uk,mEj m+kU† k,mEj mUjk,m +ak−Ej mUj,mEk m+jU† j,mUjk,m+Ek mEj mUjk,m. (36) Here, we have utilized the evolution equation of the plaquette U−kj,m, which can be derived in analogy to Eq. (9), ∂0ReTr(itaU−kj,m) =igReTritaajU† k,m−kEj m−kUk,m−kU−kj,m−U−kj,mEj m −itaakU† k,m−kEk m−kUk,m−kU−kj,m −U† k,m−kUj,m−kEk m+j−kU† j,m−kUk,m−kU−kj,m.(37) 5These notations should be understood as abbreviated forms of (∂0T1 0k)Eand (∂0T1 0k)B, i.e., the E-field and B-field parts of the time derivative, rather than time derivatives of E-field and B-field parts. By employing the expression in Eq. (20), we can rewrite Eq. (36)asfollows: ∂0T1E 0k,m= j=k gcjkReTr ×ajEj m+kEj m+kU† k,mUjk,mUk,m−Ej mEj mUjk,m +akEj m+kU† k,mEk mUjk,mUk,m −Ej m+kU† k,mUj,mEk m+jU† j,mUjk,mUk,m −Ej mUj,mEk m+jU† j,mUjk,m+Ek mEj mUjk,m. (38) Note that in the continuum limit, only the term where the plaquette on the r.h.s is replaced by an identity matrix survives. On the lattice, however, it is not possible to rewrite derivatives of the plaquette in such a way that they would not themselves involve plaquettes. In contrast to the previous subsection, where obtaining T0ifrom ∂0T00 was straightforward, this part is more challenging. In Fig. 2, we illustrate the issue by examining the time derivative of T1 0k(see Eq. 35) at lattice point m.The components involving two plaquettes i.e the first line on the right-hand side depicts the magnetic field component that we will discuss in a later section, while the second line refers to the electric field part ∂0T1E 0k. Notably, each term entering the latter includes an extra plaquette, as explicitly stated in Eq. (36). This differs from the continuum expression TE jk =EjEk+1 2δjkE2, where the electric terms are solely quadratic in nature and not entangled with magnetic field components that are encoded in the plaquettes. We have not found an exact way of writing ∂0T1E 0kas a total spatial derivative of a quantity that could be identified as a contribution to Tij. We, therefore, use the approximation of replacing the plaquettes in ∂0T1E 0kwith the identity matrix. This introduces a relative error O(a2)and agrees with the original expression in the continuum limit. After rearranging certain terms and adding the contribution for ∂0T1E 0k|m→m−j (as in Eq. 29) the resulting expression takes the following form: ∂0TE 0k,m= j=k gcjkReTrajEj m+kEj m+k −Ej mEj m+Ej m+k−jEj m+k−j−Ej m−jEj m−j −akajEj mDF j,Ek m +Uk,mEj m+kU† k,mDF j,Ek m]+Ej m−jDF j,Ek m−j +Uk,m−jEj m+k−jU† k,m−jDF j,Ek m−j]+O(a2). (39) 123 368 Page 8 of 18 Eur. Phys. J. C (2024) 84:368 Fig. 2 Illustration of the time derivative of T1 0kin Eq. (35). The filled yellow circle symbolizes the lattice point, while the electric field at a specific point is indicated by red lines. The initial four terms on the right-hand side represent the magnetic field component of ∂0T1B 0kas defined in Eq. (53), while the subsequent eight terms on the following line correspond to the electric field component of ∂0T1E 0kas defined in Eq. (36) It is important to highlight that when jis equal to k,different terms in the first and second lines of ∂0TE 0kcancel each other. This permits us to interchange the summations j=k and j,kin the above equation without altering the final outcome. Having established this, we will apply the product rule, allowing us to extract the spatial derivative along the jdirection. This step is essential to obtain a discretized second (momentum conservation) continuity equation, ultimately enabling us to identify TE jk. The discretized form of the product rule is given as DF k,Ei mEj m =DF k,Ei mUk,mEj m+kUk,m+Ei mDF k,Ej m,(40) where we note the shift from site mto site m+kin the second electric field factor in the first term. In the case of a parallel transported field, this can be written as DF k,U† k,m−kEi m−kUk,m−kEj m =[DF k,U† m−kEi m−kUk,m−k]Ej m+Ei m[DF k,Ej m].(41) We utilize either of these relations to rewrite the terms in Eq. (39), which enables us to employ Gauss’ law and eliminate certain contributions from the equation [DB j,U† j,mEj mUj,mEk m+j]=[DF j,U† j,m−jEj m−jUj,m−jEk m] =Ej mUj,mEk m+jU† j,m −Ej mEk m−U† j,m−jEj m−jUj,m−jEk m+Ej mEk m =Ej m[DF j,Ek m]+Ek m[DB j,Ej m]   =0GaussLaw (42) [DB j,U† j,mUk,mEj m+kU† k,mUj,mEk m+j] =DF j,U† j,m−jUk,m−jEj m−j+kU† k,m−jUj,m−jEk m =Uk,mEj m+kU† k,mUj,mEk m+jU† j,m−Uk,mEj m+kU† k,mEk m −Uk,m−jEj m+k−jU† k,m−jUj,m−jEk mU† j,m−j +Uk,mEj m+kU† k,mEk m =Uk,mEj m+kU† k,mDF j,Ek m+DB j,Uk,mEj m+kU† k,mEk m (43) DB j,Ej mEk m=DF j,Ej m−jEk m−j =Ej m−j[DF j,Ek m−j]+Ek m[DB j,Ej m]   =0GaussLaw (44) DB j,Uk,mEj m+kU† k,mEk m =DF j,Uk,m−jEj m+k−jU† k,m−jEk m−j =Uk,m−jEj m+k−jU† k,m−jDF j,Ek m−j +DB j,Uk,mEj m+kU† k,mEk m(45) We can make slight modifications to the underlined terms in Eqs. (43) and (45) to make use of Gauss’ law. DB j,Uk,mEj m+kU† k,mEk m =Uk,mEj m+kU† k,mEk m−U† j,m−j ×Uk,m−jEj m+k−jU† k,m−jUj,m−jEk m =Uk,mEj m+kU† k,mEk m−U−jk,mUk,mUj m−j+kEj m−j+k ×U−kj,m−j+kUj,m−j+kU† k,mEk m ≃Uk,mEj m+kU† k,mEk m−1Uk,mUj m−j+kEj m−j+k ×1Uj,m−j+kU† k,mEk m =Uk,m[DB j,Ej m+k]U† k,mEk m =0[Gauss Law](46) As we transition from the first to the second equality, we introduce additional gauge links in the second term on the righthand side, thereby forming plaquettes denoted as U−jk,mand U−kj,m−j+k. Subsequently, in the third equality, we approx- 123 Eur. Phys. J. C (2024) 84:368 Page 15 of 18 368 Fig. 5 Relative error (72) in the momentum conservation equation as a function of time for Ns=32 on the left and Ns=256 on the right. Both are shown for the electric and magnetic field contributions in the continuum (14) and discretized (63) formulations Fig. 6 Measuring the violation (72) in the second continuity equation (34) by utilizing TB ij derived from the energy density T00 (63)andthe conjectured Bianchi identity (60), considering two distinct lattice sizes: Ns=32 and Ns=128. Fig. 7 Relative error (72) in the second continuity equation (34)asa function of lattice spacing a2 sfor discretized (63) and naive (continuum limit) electric and magnetic fields (14). Black dashed lines correspond to the a2 spower-law approach. This justifies our choice of the approach chosen in (63). In future research, it would be beneficial to develop a more refined expression for the Bianchi identity to ensure that the derivation of a discretized Tij aligns seamlessly with the continuum case. The denominator of the relative error (72) can be understood as the space-integrated squared rate of change of energy flux in the i-direction. Thus it corresponds to a physical quantity, which takes a nonzero value at the continuum limit. Hence, the relative error should tend to zero in the continuum limit, and the asscaling of the ratio is that of the numerator (discretization errors of the denominator represent a subleading correction to this behavior). Figure7illustrates this quantity for the naive and improved discretizations as functions of the quadratic lattice spacing a2 s. The black dashed line corresponds to a power law ∼a2 s. We observe that in the limit of small lattice spacing, the relative error of the naive expression goes to zero slower than a2 s. In contrast to this, the improved discrete expressions follow the a2 spower law very closely, thus indicating that they are O(a2 s). The discretized version of the second continuity equation (34) would also benefit from the time-symmetrizing procedure that was performed on T00 and T0iabove. However, while deriving Eq. (62), we have already performed an error O(a2 s)when replacing plaquettes with identity operators. Similarly, the approximations in the magnetic sector are at most of the same accuracy. Artifacts arising from the time discretization are typically subleading to the spatial discretization effects since we employ dt/as1 to guarantee numerical stability. Hence, we do not consider corrections to time discretization for the momentum conservation equation. 5 Conclusions The use of numerical simulations to gain insights into nonperturbative aspects of classical and quantum field theories requires discretizing space-time on a lattice and demands a 123 368 Page 16 of 18 Eur. Phys. J. C (2024) 84:368 systematic approach to study the energy–momentum tensor for nonabelian gauge fields. We present an improved expression for the this purpose, given by Eqs. (18), (29), and (63), that are derived using classical field equations of motion in conjunction with the energy–momentum conservation law for Tμν. In comparison to a naive discretization method, where chromoelectric and -magnetic fields are replaced by lattice counterparts, our formulation improves the relative violation of the conservation laws by several orders of magnitude. The energy conservation equation (17) offers a means of obtaining the T0kcomponents of the energy–momentum tensor that satisfies an exact energy conservation relation on the lattice. Challenges arise when deriving Tjk using the momentum conservation equation (34). The terms involving electric fields introduce spurious contributions like EjEjBk that are suppressed by the lattice spacing and are not present in the continuum expression. Furthermore, the terms involving magnetic fields cannot be written as a spatial derivative of Tjk due to the lack of a suitable lattice Bianchi identity. We have avoided this by replacing elements of TB jk that are proportional to δjk with parts of the discrete energy density, which has led to a significant reduction of the violation. In the future, our focus will be on addressing these issues to obtain the subleading terms O(a2 s)in the expressions of Tjk. For the energy-conserving continuity expression, we also see that the relative error due to finite timesteps can be further improved by orders of magnitude when using a timesymmetrized discretization of the equation. As illustrated by Eqs. (32) and (33), the electric and magnetic fields then lie on the same time slice in the standard leapfrog algorithm. We expect our work to have several interesting applications, especially with an extension to account for small perturbations on top of a nonequilibrium plasma [42–44]. Examples of these are transport coefficients, e.g., shear viscosity in an over-occupied gluon plasma. In this context, the energy conservation given by the continuity equation can hopefully be used to prevent the activation of other modes, like sound modes in the case of shear viscosity, ensuring a more accurate depiction of the system’s behavior. We have made a first step toward a direct measurement of such transport coefficients in App. B, where we have derived an expression for the perturbed energy–momentum tensor after introducing small fluctuations. Another intriguing future research direction is to expand our work to a wider range of metric tensors, such as the Friedmann–Lemaître–Robertson–Walker (FLRW) metric and a longitudinally expanding (Bjorken) metric in the contexts of cosmology and heavy-ion collisions, respectively. Indeed, including expansion in the framework would permit a more realistic treatment of the Glasma at the initial stages of heavy-ion collisions. Furthermore, it would extend the applicability of our framework to cosmological applications as well. Acknowledgements We would like to thank D. I. Müller, H. Matsuda, and S. Schlichting for valuable discussions. This work is supported by the European Research Council, ERC-2018-ADG-835105 YoctoLHC. This work was also supported by the European Union’s Horizon 2020 research and innovation by the STRONG-2020 project (grant agreement No. 824093). TL, JP, and PS have been supported by the Academy of Finland, by the Centre of Excellence in QuarkMatter (project 346324), and project 321840. KB would like to thank the Austrian Science Fund (FWF) for support under project P 34455. The authors wish to acknowledge the Vienna Scientific Cluster (VSC) under project 71444 and the CSC-IT Center for Science Finland, for computational resources on the supercomputer Puhti. Data availability This manuscript has no associated data or the data will not be deposited. [Authors’ comment: Data sharing not applicable to this article as no datasets were generated or analysed during the current study.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP3. Appendix A: Derivation of a relation equivalent to the Bianchi Identity The general proof of Eq. (51) can be given as δij 2[Dj,Fmn]Fmn =δijδmm 2[Dj,Fmn]Fmn =1 2kmjkmi+δmiδmj[Dj,Fmn]Fmn =1 2kmjkmi[Dj,Fmn]Fmn+1 2[Dj,Fin]Fjn (A1) where the first term can be further modified as 1 2kmjkmi[Dj,Fmn]Fmn =− 1 2kmjkmimnl[Dj,Bl]Fmn =1 2mkjmnl[Dj,Bl]Fmnkmi =1 2δknδjl −δklδjn[Dj,Bl]Fmnkmi 123 Eur. Phys. J. C (2024) 84:368 Page 17 of 18 368 =1 2[Dj,Bj]Fmnnmi−1 2[Dj,Bk]Fmjkmi =1 2[Dj,Fmi]Fmj =1 2[Dj,Fin]Fjn (A2) On combining the above two expressions, one recovers Eq. (51). Appendix B: Energy momentum tensor for fluctuations Over the years, significant effort has been devoted to investigating the linear response of non-Abelian plasma, focusing on the analysis of fluctuations superimposed on background fields [42–44]. Calculating the linear response involves decomposing the gauge field and electric field as follows: Ai(x)→Ai(x)+ai(x)(B1) Ei(x)→Ei(x)+ei(x)(B2) The link matrix representing the combination of background and fluctuating fields is expressed as UBG+fluct i,m=1+igai,mas iUi,m(B3) resulting in the following form of the plaquette: UBG+fluct jk,m=Ujk,m+δUjk,m =Ujk,m+igaj,mas jUjk,m+Uj,mak,m+jas kU† j,mUjk,m −Ujk,mUk,maj,m+kas jU† k,m−Ujk,mak,mas k(B4) where as iis the lattice spacing along the idirection. Following this decomposition, one can formulate expressions for the energy–momentum tensor, such as: δT00,m= j,k>0 −1 g2a2 ja2 k 1 4ReTrδUjk,m +δUjk,m−j+δUjk,m−k+δUjk,m−j−k +1 4 j>0 2ej,a mEj,a m+2ej,a m−jEj,a m−j(B5) Similarly, the expression for the density along the k-th component of linear momentum can be provided as: δT0k,m= j>0 1 2g ajak a2 ja2 kReTriej m+kU−kj,m+k +iEj m+kδU−kj,m+k+iej mUjk,m+iEj mδUjk,m +iej m−j+kU−kj,m−j+k+iEj m−j+kδU−kj,m−j+k +iej m−jUjk,m−j+iEj m−jδUjk,m−j(B6) References 1. K.G. Wilson, Confinement of quarks. Phys. Rev. D 10, 2445 (1974). https://doi.org/10.1103/PhysRevD.10.2445 2. F. Gelis, E. Iancu, J. Jalilian-Marian, R. Venugopalan, The color glass condensate. Annu. Rev. Nucl. Part. Sci. 60, 463 (2010). https://doi.org/10.1146/annurev.nucl.010909.083629. arXiv:1002.0333 [hep-ph] 3. F. Gelis, Color glass condensate and glasma. Int. J. Mod. Phys. A 28, 1330001 (2013). https://doi.org/10.1142/ S0217751X13300019.arXiv:1211.3327 4. T. Lappi, L. McLerran, Some features of the glasma. Nucl. Phys. A 772, 200 (2006). https://doi.org/10.1016/j.nuclphysa.2006.04.001. arXiv:hep-ph/0602189 5. B. Schenke, P. Tribedy, R. Venugopalan, Fluctuating Glasma initial conditions and flow in heavy ion collisions. Phys. Rev. Lett. 108, 252301 (2012). https://doi.org/10.1103/PhysRevLett.108.252301. arXiv:1202.6646 [nucl-th] 6. B. Schenke, P. Tribedy, R. Venugopalan, Event-by-event gluon multiplicity, energy density, and eccentricities in ultrarelativistic heavy-ion collisions. Phys. Rev. C 86, 034908 (2012). https://doi. org/10.1103/PhysRevC.86.034908.arXiv:1206.6805 [hep-ph] 7. T. Lappi, S. Schlichting, Linearly polarized gluons and axial charge fluctuations in the Glasma. Phys. Rev. D 97(3), 034034 (2018). https://doi.org/10.1103/PhysRevD.97.034034.arXiv:1708.08625 [hep-ph] 8. J.L. Albacete, P. Guerrero-Rodríguez, C. Marquet, Initial correlations of the Glasma energy–momentum tensor. JHEP 01, 073 (2019). https://doi.org/10.1007/JHEP01(2019)073. arXiv:1808.00795 [hep-ph] 9. H. Fujii, K. Fukushima, Y. Hidaka, Initial energy density and gluon distribution from the Glasma in heavy-ion collisions. Phys. Rev. C 79, 024909 (2009). https://doi.org/10.1103/PhysRevC.79.024909. arXiv:0811.0437 [hep-ph] 10. J. Berges, K. Boguslavski, S. Schlichting, R. Venugopalan, Turbulent thermalization process in heavy-ion collisions at ultrarelativistic energies. Phys. Rev. D 89(7), 074011 (2014). https://doi.org/10. 1103/PhysRevD.89.074011.arXiv:1303.5650 [hep-ph] 11. J. Berges, K. Boguslavski, S. Schlichting, R. Venugopalan, Universal attractor in a highly occupied non-Abelian plasma. Phys. Rev. D89(11), 114007 (2014). https://doi.org/10.1103/PhysRevD.89. 114007.arXiv:1311.3005 [hep-ph] 12. K. Boguslavski, A. Kurkela, T. Lappi, J. Peuron, Highly occupied gauge theories in 2+1 dimensions: a self-similar attractor. Phys. Rev. D 100(9), 094022 (2019). https://doi.org/10.1103/PhysRevD. 100.094022.arXiv:1907.05892 [hep-ph] 13. M.E. Carrington, A. Czajka, S. Mrówczy´nski, Physical characteristics of glasma from the earliest stage of relativistic heavy ion collisions. Phys. Rev. C 106(3), 034904 (2022). https://doi.org/10. 1103/PhysRevC.106.034904.arXiv:2105.05327 [hep-ph] 14. M.E. Carrington, A. Czajka, S. Mrowczynski, The energymomentum tensor at the earliest stage of relativistic heavy-ion collisions. Eur. Phys. J. A 58(1), 5 (2022). https://doi.org/10.1140/ epja/s10050-021-00600-x.arXiv:2012.03042 [hep-ph] 15. A. Ipp, D. Müller, Broken boost invariance in the Glasma via finite nuclei thickness. Phys. Lett. B 771, 74 (2017). https://doi.org/10. 1016/j.physletb.2017.05.032.arXiv:1703.00017 [hep-ph] 16. S. McDonald, S. Jeon, C. Gale, The 3+1D initialization and evolution of the Glasma. arXiv:2306.04896 [hep-ph] 123 368 Page 18 of 18 Eur. Phys. J. C (2024) 84:368 17. B. Schenke, S. Schlichting, 3D glasma initial state for relativistic heavy ion collisions. Phys. Rev. C 94(4), 044907 (2016). https://doi.org/10.1103/PhysRevC.94.044907.arXiv:1605.07158 [hep-ph] 18. S. Schlichting, P. Singh, 3-D structure of the Glasma initial state–breaking boost-invariance by collisions of extended shock waves in classical Yang-Mills theory. Phys. Rev. D 103(1), 014003 (2021). https://doi.org/10.1103/PhysRevD.103.014003. arXiv:2010.11172 [hep-ph] 19. A. Ipp, D.I. Müller, S. Schlichting, P. Singh, Spacetime structure of (3+1)D color fields in high energy nuclear collisions. Phys. Rev. D 104(11), 114040 (2021). https://doi.org/10.1103/PhysRevD.104. 114040.arXiv:2109.05028 [hep-ph] 20. B. Schenke, S. Schlichting, P. Singh, Rapidity dependence of initial state geometry and momentum correlations in p+Pb collisions. Phys. Rev. D 105(9), 094023 (2022). https://doi.org/10. 1103/PhysRevD.105.094023.arXiv:2201.08864 [nucl-th] 21. H. Matsuda, X.-G. Huang, Simulation of 3+1D glasma in Milne coordinates I: development of the framework. arXiv:2308.15269 [hep-ph] 22. A. Ipp, M. Leuthner, D.I. Müller, S. Schlichting, K. Schmidt, P. Singh, Energy-momentum tensor of the dilute (3+1)D Glasma. arXiv:2401.10320 [hep-ph] 23. K. Boguslavski, A. Kurkela, T. Lappi, J. Peuron, Heavy quark diffusion in an overoccupied gluon plasma. JHEP 09, 077 (2020). https:// doi.org/10.1007/JHEP09(2020)077.arXiv:2005.02418 [hep-ph] 24. D. Avramescu, V. B˘aran, V. Greco, A. Ipp, D.I. Müller, M. Ruggieri, Simulating jets and heavy quarks in the glasma using the colored particle-in-cell method. Phys. Rev. D 107(11), 114021 (2023). https://doi.org/10.1103/PhysRevD.107.114021. arXiv:2303.05599 [hep-ph] 25. A. Ipp, D.I. Müller, D. Schuh, Jet momentum broadening in the preequilibrium Glasma. Phys. Lett. B 810, 135810 (2020). https://doi. org/10.1016/j.physletb.2020.135810.arXiv:2009.14206 [hep-ph] 26. M.E. Carrington, A. Czajka, S. Mrowczynski, Heavy quarks embedded in glasma. Nucl. Phys. A 1001, 121914 (2020). https:// doi.org/10.1016/j.nuclphysa.2020.121914.arXiv:2001.05074 [nucl-th] 27. Pooja, M. Ruggieri, S.K. Das, Heavy quark diffusion in glasma and gluonic plasma. DAE Symp. Nucl. Phys. 66, 924 (2023) 28. S.K. Pooja, V. Das, M. Ruggieri, Anisotropic fluctuations of angular momentum of heavy quarks in the Glasma. Eur. Phys. J. Plus 138(4), 313 (2023). https://doi.org/10.1140/epjp/ s13360-023-03913-6.arXiv:2212.09725 [hep-ph] 29. M.E. Carrington, A. Czajka, S. Mrowczynski, Transport of hard probes through glasma. Phys. Rev. C 105(6), 064910 (2022). https://doi.org/10.1103/PhysRevC.105.064910. arXiv:2202.00357 [nucl-th] 30. M.E. Carrington, A. Czajka, S. Mrowczynski, Jet quenching in glasma. Phys. Lett. B 834, 137464 (2022). https://doi.org/10.1016/ j.physletb.2022.137464.arXiv:2112.06812 [hep-ph] 31. T. Lappi, Production of gluons in the classical field model for heavy ion collisions. Phys. Rev. C 67, 054903 (2003). https://doi.org/10. 1103/PhysRevC.67.054903.arXiv:hep-ph/0303076 32. C. Shen, B. Schenke, Longitudinal dynamics and particle production in relativistic nuclear collisions. Phys. Rev. C 105(6), 064905 (2022). https://doi.org/10.1103/PhysRevC.105.064905. arXiv:2203.04685 [nucl-th] 33. S. Caracciolo, G. Curci, P. Menotti, A. Pelissetto, The energy momentum tensor for lattice gauge theories. Ann. Phys. 197, 119 (1990). https://doi.org/10.1016/0003-4916(90)90203-Z 34. S. Caracciolo, P. Menotti, A. Pelissetto, One loop analytic computation of the energy momentum tensor for lattice gauge theories. Nucl. Phys. B 375, 195 (1992). https://doi.org/10.1016/ 0550-3213(92)90339-D 35. M. Bochicchio, L. Maiani, G. Martinelli, G.C. Rossi, M. Testa, Chiral symmetry on the lattice with Wilson fermions. Nucl. Phys. B 262, 331 (1985). https://doi.org/10.1016/0550-3213(85)90290-1 36. A. Tranberg, G. Ungersbäck, Bubble nucleation and quantum initial conditions in classical statistical simulations. JHEP 09, 206 (2022). https://doi.org/10.1007/JHEP09(2022)206. arXiv:2206.08691 [hep-lat] 37. A. Arrizabalaga, J. Smit, A. Tranberg, Tachyonic preheating using 2PI-1/N dynamics and the classical approximation. JHEP 10, 017 (2004). https://doi.org/10.1088/1126-6708/2004/10/017. arXiv:hep-ph/0409177 38. M. D’Onofrio, K. Rummukainen, A. Tranberg, Sphaleron rate in the minimal standard model. Phys. Rev. Lett. 113(14), 141602 (2014). https://doi.org/10.1103/PhysRevLett.113.141602. arXiv:1404.3565 [hep-ph] 39. J.B. Kogut, L. Susskind, Hamiltonian formulation of Wilson’s lattice gauge theories. Phys. Rev. D 11, 395 (1975). https://doi.org/ 10.1103/PhysRevD.11.395 40. H. Matsuda, T. Kunihiro, B. Müller, A. Ohnishi, T.T. Takahashi, Shear viscosity of classical Yang–Mills field. Phys. Rev. D102, 114503 (2020). https://doi.org/10.1103/PhysRevD.102. 114503.arXiv:2007.06886 [hep-ph] 41. J.E. Kiskis, The Bianchi identity for nonabelian lattice gauge fields. Phys. Rev. D 26, 429 (1982). https://doi.org/10.1103/PhysRevD. 26.429 42. K. Boguslavski, A. Kurkela, T. Lappi, J. Peuron, Spectral function for overoccupied gluodynamics from real-time lattice simulations. Phys. Rev. D 98(1), 014006 (2018). https://doi.org/10.1103/ PhysRevD.98.014006.arXiv:1804.01966 [hep-ph] 43. A. Kurkela, T. Lappi, J. Peuron, Time evolution of linearized gauge field fluctuations on a real-time lattice. Eur. Phys. J. C 76(12), 688 (2016). https://doi.org/10.1140/epjc/ s10052-016-4523-9.arXiv:1610.01355 [hep-lat] 44. K. Boguslavski, A. Kurkela, T. Lappi, J. Peuron, Broad excitations in a 2+1D overoccupied gluon plasma. JHEP 05, 225 (2021). https://doi.org/10.1007/JHEP05(2021)225.arXiv:2101.02715 [hep-ph] 123