Conserved energy–momentum tensor for real-time lattice simulations

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 (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal {O}}(a^2)$$\end{document}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.


I. 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 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 space-time 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.
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, 22-24, 30, 31], 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 [32][33][34] 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 capture 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 T 00 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 T 0i are then obtained by requiring that a discretized continuity equation for energy conservation be satisfied.This procedure leads to T 0i given 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 Sec.III B 2. For the spatial components T ij , we start 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 T ij .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 heavy-ion collisions are real-time lattice simulations performed in the context of postinflatory cosmology using the classical approximation [35][36][37].
This paper is structured as follows.Sec.II describes our discretization framework.In Sec.III we construct the discrete energy-momentum tensor and introduce its time-symmetrized formulation.Sec.IV shows the numerical comparison among the naive and discretized approaches.We conclude in Sec.V.

II. 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 [38] of classical SU(N c ) Yang-Mills theory in temporal axial gauge (A t = 0), Here m denotes the lattice site on the spatial grid and a = 1, . . ., N 2 c − 1 is the color index.In this paper, we consider the Minkowski metric with g µν = (+1, −1, −1, −1) and √ −g = −det(g µν ) = 1.The coupling constant g enters in the form of the gauge coupling 2N c /g 2 .The discretization is performed on a three-dimensional lattice with size N x × N y × N z , lattice spacing a µ in the spatial direction μ, and the product a 3 s = a x a y a z .To guarantee gauge invariance, the theory is formulated in terms of gauge links U i,m = e iaigAi,m instead of gauge fields A i,m .Links enter Eq. ( 2) in terms of plaquettes where μ denotes a unit vector in the µ direction. 1In analogy to the links, plaquettes are related to the field-strength tensor F µν via U µν ≈ e igaµaν Fµν .Links and plaquettes are group elements U µ,m , U µν,m ∈ SU(N c ) in the fundamental representation, and correspondingly we will use the generators of the fundamental representation t a to go between the color components of the electric field and the fundamental representation matrix E µ m = E µ,a m t a .Using the (continuum) relation between the field-strength tensor and the chromomagnetic field B i = − 1 2 ϵ ijk F jk , the Hamiltonian (2) reduces to the correct expression d 3 x 1 2 j>0 (E j,a E j,a + B j,a B j,a ) in the continuum limit a µ → 0. The Hamiltonian (2) is associated with 2 the Hamiltonian equations of motion for the link matrices and the electric field variables, which are discrete in space but continuous in time 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 and introduced the backward and forward gauge covariant derivatives 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

III. 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 = m a 3 s √ −gL with respect to the metric tensor g µν , On the lattice, our starting point is a situation where the metric is purely diagonal.Indeed, both the Hamiltonian (1) and the corresponding action S only 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 Sec.III A 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 We will start in Sec.III B by allowing spatial redistribution for the energy density, i.e., the temporal component T 00 while requiring that it sums up to the Hamiltonian that corresponds to the total energy (2).The components T 0i are then constructed from (11) for ν = 0.The remaining components will then be constructed in Sec.III C in the same way by taking T 0i as the starting point and requiring that T ij satisfy the remaining continuity equations for ν > 0. This procedure determines T 0i up to a transformation T 0i → T 0i + ϕ i , where ∂ i ϕ i = 0.This transformation, however, leaves conserved quantities intact.

A. 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 E and B fields with their spatially averaged discretized counterparts where E 2 m,loc = E i,a m,loc E i,a m,loc and E i,a m,loc and B i,a m,loc are local electric and magnetic (cloverleaf) fields.Since the electric field labeled as E i,a m corresponds to the time derivative of the gauge field between m and m + î, it is rather centered at the point m + î/2.Similarly, a plaquette U ij,m is actually centered at m + î/2 + ĵ/2.Thus a natural  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,n but 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/2 to t − dt/2.
way to construct electric and magnetic fields at a site m is to take nearest neighbor averages to obtain symmetric expressions:

B. Continuity equation for energy conservation
Our starting point will be the scalar component of the continuity equation which will be used to obtain T 0i after T 00 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 E i m is located at the position m + î 2 + dt 2 , taking also the finite timestep dt into account.The magnetic field strength, expressed by the plaquette U ij,m , is located at m + î 2 + ĵ 2 .To determine the energy density at lattice site m, we compute an average over the "outgoing" and "incoming" electric fields E j,a m and E j,a m−j and various plaquette orientations.The left panel of Fig. 1 illustrates 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 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.
Taking a square of the symmetrized electric field from Eq. ( 15) as is done in [15,18,21,23,24,39] on the other hand, is not equivalent to the actual Hamiltonian (2): 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 T 00 , the equation of motion (4) would introduce cubic terms in the electric fields into T 0i , 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 and the fact that the parallel transporting links cancel in the trace Tr(U † M U ) = TrM .This allows us to rewrite the energy density T 00,m as

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 T 0i .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): Indeed, it is possible to write the right-hand side as a total spatial derivative.After some rearrangement, one achieves the following result: 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 Thus we have arrived at the first important result of this paper, a discrete energy conservation law 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: where It is important to note that what we label by T 0k,m is actually centered at the position m + k 2 ,4 which 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,m T 0k , is in fact centered around m, the same position where ∂ 0 T 00 is defined.
We also note that Eq. ( 29) corresponds to a gauge covariant formulation of the Poynting vector in the continuum limit T 0k → ϵ kij E i B j .This can be seen by expanding the plaquettes to quadratic order in a s , thus U jk ≈ 1 + iga j a k F jk , and realizing that the term with the identity vanishes because the matrix E j is traceless.

Symmetric time discretization
Until now, we have treated t as a continuous variable, which is the limit dt/a s ≪ 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 ∂ 0 T 00 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 T 00 expression In the leapfrog scheme, the time difference labeled above, corresponding to stepping the electric fields from t − 3dt/2 to t − 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 T 0k should be evaluated with the following timestep assignment in terms of a time-symmetrized electric field: FIG. 2. Illustration of the time derivative of T 1 0k in Eq. ( 29).The filled yellow circle symbolizes the lattice point, while the electric field at a specific point is indicated by red lines.
Similarly, the time derivative of the magnetic field part of the energy density corresponds to the term in ∂ k T 0k where one takes a spatial derivative of the electric field.A time-symmetric treatment of this term requires a symmetrization of the links in the evaluation of this term in where the notation D F 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.

C. Continuity equation for momentum conservation
Our goal is to construct T jk in the same way as we constructed T 0i above, i.e., by utilizing the equations of motion and the continuity equation for i > 0 to derive the T ij components from ∂ 0 T 0i .We first observe that T 0k comprises two kinds of terms, with one of them being merely shifted in the − ĵ direction.Henceforth, we will focus solely on the first one of the two, which we call T 1 0k for our subsequent calculations.The second term can then be restored in the end by shifting the site of T 1 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.
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 T ij , which we refer to here as ∂ 0 T 1E 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 ∂ 0 T 1B 0k . 5

Chromoelectric field contribution
Let us start with the chromoelectric field contribution ∂ 0 T 1E 0k .We begin by employing the equation of motion for the plaquette to express the electric field part ∂ 0 T 1E 0k as Here, we have utilized the evolution equation of the plaquette U −kj,m , which can be derived in analogy to Eq. ( 9), 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.By employing this expression in Eq. ( 20), we can rewrite Eq. ( 36) as follows: In contrast to the previous subsection, where obtaining T 0i from ∂ 0 T 00 was straightforward, this part is more challenging.In Fig. 2, we illustrate the issue by examining the time derivative of T 1 0k (see Eq. 29) at lattice point m.The components involving two plaquettes depict the magnetic field component that we will discuss in a later section, while the second line refers to the electric field part ∂ 0 T 1E 0k .Notably, each term entering the latter includes an extra plaquette, as explicitly stated in Equation (38).This differs from the continuum expression T E jk = E j E k + 1 2 δ jk E 2 , 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 ∂ 0 T 1E 0k as a total spatial derivative of a quantity that could be identified as a contribution to T ij .We, therefore, use the approximation of replacing the plaquettes in ∂ 0 T 1E 0k with the identity matrix.This introduces a relative error O(a 2 ) and agrees with the original expression in the continuum limit.After rearranging certain terms and adding the contribution for ∂ 0 T 1E 0k | m→m−j (as in Eq. 29) the resulting expression takes the following form: It is important to highlight that when j is equal to k, different terms in the first and second lines of ∂ 0 T E 0k cancel each other.This permits us to interchange the summations j̸ =k and j,k in 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 j direction.This step is essential to obtain a discretized second (momentum conservation) continuity equation, ultimately enabling us to identify T E jk .The discretized form of the product rule is given as where we note the shift from site m to site m + k in the second electric field factor in the first term.In the case of a parallel transported field, this can be written as 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 We can make slight modifications to the underlined terms in Eqns.( 43) and (45) to make use of Gauss' law.
As we transition from the first to the second equality, we introduce additional gauge links in the second term on the right-hand side, thereby forming plaquettes denoted as U −jk,m and U −kj,m−j+k .Subsequently, in the third equality, we approximate the plaquettes with the identity, similarly as discussed earlier, and at the same relative O(a 2 ) error, which enables us to apply Gauss' law.By employing Eqns.( 42) -(45), we can rephrase Eq. ( 39) as follows: This allows us to use the second continuity equation (34), to identify the electric field component of the stress tensor as follows: The contribution T E jk is located at the position6 m + j 2 + k 2 , precisely as expected from ∂ 0 T 0k calculations.

Chromomagnetic field contribution
Having found an expression for the electric contribution in (48), our next objective is to derive an expression for the magnetic field component of T jk .However, before delving into that calculation, it is beneficial to examine the continuum expression to illustrate the issues we will encounter along the way.For this discussion to be closer to the discretized version, we will formulate it in terms of the field strength tensor instead of the magnetic field.We begin by computing the magnetic field term in the time derivative of Comparing this expression to the spatial derivative of the spacelike part of the energy-momentum tensor T B ij = δij 4 F pq F pq − F iq F jq , it is clear that additional terms need to be added to and subtracted from (49) to write it as a spatial derivative.In fact, to see how this happens in the continuum, it is easier to start from T B ij = δij 4 F pq F pq −F iq F jq and differentiate it: In obtaining the last equality, we have utilized the relation B r = − 1 2 ϵ rpq F pq and the Bianchi identity [D j , B j ] = 0.The essential part of this manipulation was the cancellation of the first two terms on the right-hand side of the first line of Eq. ( 50): Appendix A presents an alternative derivation of this same relation.An important part of the derivation was the Bianchi identity.A version of the Bianchi identity exists on the lattice [40].It involves the product of six plaquettes, which, in the continuum, reduces to the continuum Bianchi identity.However, the identity we would need here would involve discretized derivatives of plaquettes, i.e., their nearest-neighbor differences.We have not been able to formulate a suitable exact discrete identity that could play the role of the Bianchi identity in the derivation of the energy-momentum tensor.Thus again, as in the case of the electric field component, we have to resort to an approximation that reduces to (51) in the continuum limit, Within this approximation, we are set to calculate the magnetic field contribution of T jk by employing the dynamical equations for the field on the ∂ 0 T 1B 0k component of Eq. (35).These contributions are visualized in the first line of Fig. 2, which shows the terms that appear after plugging in the temporal derivatives of the electric fields , where Transitioning from the second equality to the third, we took into account ReTr(it a 1) = 0 under the conditions j = k and j = l, and subsequently exchanged the indices j and l.Here, the notation T B jk,m,CB is employed to represent the magnetic-field component of the stress tensor derived from the conjectured Bianchi (CB) identity.Note that in order to write the δ jk -part of T B jk explicitly in this form, we have used the continuum form of the product rule, rather than the exact one (40), and the approximate Bianchi identities (52).
We have checked the quality of the conjectured approximate Bianchi identity (52) approximation numerically, using the setup that will be discussed in more detail in Sec.IV.The result is demonstrated in Fig. 3.Here we show the relative error in the lattice counterpart of the conjectured Bianchi identity, calculated as at distinct lattice sites labeled as A, B and C within the transverse plane.It is evident that the relative error is large and exhibits local variations, thereby suggesting the potential for devising an improved representation of the lattice Bianchi identity in future studies.
For now, we need an alternative, more accurate expression.Since our conjectured Bianchi identity is not satisfied well for the field configurations that we would like to apply it to, we will use another approach that circumvents the main issues.Here, one first notices that in our derivation leading to T B jk , the Bianchi identity has been used only to manipulate the total energy contribution that is proportional to δ jk .In the continuum, the magnetic field part of the δ jk -term is just the same magnetic field squared that appears in the energy density.Thus, we will follow another approach that does not rely on Eq. ( 52), which is to just take the magnetic field part of the total energy density T B 00 and use it to get the magnetic δ jk part of T jk .This leads to This is achieved by subtracting T B 00,m+j based on the first term of the continuum expression T ij = −1/4δ ij F mn F mn + F in F jn .Additionally, this contribution from the energy density needs to be taken at the point m + j, as T jk for j = k needs to be situated at the point m + j in order to yield ∂ 0 T 0k centered at m + k/2 as a backward derivative in the j = k direction.Equation (62) thus does not use the approximate Bianchi identity and turns out to introduce only a small relative error to the continuity equation, as we will show below.At the beginning, we take the fields from initial conditions of the form where g is the coupling constant, Q s a dimensionful initial momentum scale, and a 3 s N 3 s the lattice volume.The expectation value ⟨•⟩ implies an average over classical configurations of the lattice fields.In Eq. ( 64) we only initialize transverse modes, p i A i (p, t=0) = 0.With these initial conditions, all of the energy of the system resides initially in the chromomagnetic fields, and electric fields are generated over the time evolution of the system.These initial conditions are the magnetic field part of the ones used, e.g., in [22,41].Choosing the initial condition to have a zero electric field allows Gauss' law to be exactly satisfied at the initial condition, while the leapfrog algorithm for the time evolution preserves it to machine precision. 8After the characteristic timescale t ∼ 1/Q s the energy becomes roughly evenly distributed between the electric and magnetic fields.
Within this setup, we will present our results for lattices with a constant physical volume N s a, where we change the lattice spacing a s = {0.25,0.125, 0.0625, 0.03125} and the lattice size N s = {32, 64, 128, 256} simultaneously.We will also vary the time step dt.Due to computational complexity, we will conduct all our simulations using the SU (2) gauge group instead of the physical SU (3) gauge group of QCD.

B. Results
Based on the above expressions ( 18), ( 29) and (63) of the energy-momentum tensor, we will now study the violation of the continuity equations ∂ µ T µν = 0. We start our analysis with the equation for energy conservation, where we quantify the deviation in the form of a relative error Here, the time derivative is calculated explicitly as a time difference between two discrete timesteps.Figure 4 illustrates the evolution of the relative error in three scenarios: one employing the straightforward "naive" discretization method outlined in Eqs.(12) and (13), and the others employing the discretized approach specified in Eq. ( 18) along with the Poynting vector formulated either by Eq. ( 29) or the time-symmetrized discretization outlined in Sec.III B 2. The results are presented for a spatially three-dimensional lattice with N s = 32 in the left and N s = 128 in the right panels.The top panel corresponds to dt/a s values of 0.01 and 0.1 for N s = 32 and 128, respectively.The lower panels illustrate the relative error for one-tenth of these values for dt/a s .First of all, it is evident that the violation remains relatively consistent over time, even after the energy becomes distributed over both electric and magnetic fields.Secondly, the discrete approach with time symmetrization exhibits exceptional performance for different lattice spacings and timesteps.Thirdly, for smaller lattices, the discrete method outperforms the naive description.However, as we approach the continuum limit by increasing N s and decreasing the lattice spacing accordingly, the naive description improves and approaches the discrete one as expected.
When examining the upper panel of Fig. 4 in comparison to the lower one with a smaller dt/a s ratio by a factor of ten, one observes that the relative error of the discrete approach reduces roughly by the same factor.This is again consistent with the expectation that the discrete formulation in this regime is dominated by timestep effects.However, we note that the time-symmetric discrete method is still several orders of magnitude more precise even with a rather small timestep.The error in the time-symmetrized formulation does not decrease anymore when the timestep is decreased from the top row to the bottom, which we take as an indication that it has already reached a limit where it is dominated by machine precision effects.
In Fig. 5, we examine the relative error of the second continuity equation (34), quantified by as a function of time for i = 1.The violation is measured separately for the electric and magnetic field components for N s = 32 on the left and N s = 256 on the right-hand side as before.The naive approach is based on Eqs. ( 13) and ( 14), while the discrete approach utilizes Eqs. ( 29) and (63).In the figure, one observes that the discrete formulation surpasses the naive approach.This effect is particularly large for the electric field part but is also present for the magnetic field contribution.Indeed, the violation of the electric contribution in the discrete case is roughly an order of magnitude smaller than for the magnetic one, which indicates that our approximations for the electric sector of T ij are more accurate than for the magnetic one.As we approach the continuum limit by considering the right panel, the violation further diminishes for both the naive and the discretized expressions, while other trends remain consistent.Thus, while we have not been able to achieve a description with an exact conservation law on a discrete lattice, we have managed to find a formulation where the lattice discretization effects on the conservation law are significantly smaller than for a naive discretization.As discussed in Sec.III C 2, we had two formulations for the magnetic part of T ij , one based on a conjectured approximation that reduces to the Bianchi identity in the continuum, the other one using the discrete energy density T 00 to reconstruct the part of T ij that is proportional to δ ij .Figure 6 offers a comprehensive analysis of the difference between these two approaches for the violation of the momentum conservation equation, as parametrized by the relative error (68).We investigate the evolution of the relative error for lattice sizes of N s = 32 and N s = 128 for those two distinct stress tensor formulations.We note that the expression based on T 00 (63) performs significantly better than the one derived utilizing the lattice analog of the Bianchi identity (60), demonstrating an improvement of at least one order of magnitude.As we approach the continuum limit with N s = 128, we observe a decrease in violation for both cases, as anticipated.Notably, the reduction occurs at a more favorable rate for the T 00 -based 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 T ij aligns seamlessly with the continuum case.34) as a function of lattice spacing a 2 s for discretized (63) and naive (continuum limit) electric and magnetic fields (14).Black dashed lines correspond to the a 2 s power-law.
The denominator of the relative error (68) 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 a s scaling of the ratio is that of the numerator (discretization errors of the denominator represent a subleading correction to this behavior).Figure 7 illustrates this quantity for the naive and improved discretizations as functions of the quadratic lattice spacing a 2 s .The black dashed line corresponds to a power law ∼ a 2 s .We observe that in the limit of small lattice spacing, the relative error of the naive expression goes to zero slower than a 2 s .In contrast to this, the improved discrete expressions follow the a 2 s power law very closely, thus indicating that they are O(a 2 s ).The discretized version of the second continuity equation ( 34) would also benefit from the time-symmetrizing procedure that was performed on T 00 and T 0i above.However, while deriving Eq. (62), we have already performed an error O(a 2 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/a s ≪ 1 to guarantee numerical stability.Hence, we do not consider corrections to time discretization for the momentum conservation equation.

V. 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 systematic approach to study the energymomentum 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 energymomentum 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 T 0k components of the energy-momentum tensor that satisfies an exact energy conservation relation on the lattice.Challenges arise when deriving T jk using the momentum conservation equation (34).The terms involving electric fields introduce spurious contributions like E j E j B k 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 T jk due to the lack of a suitable lattice Bianchi identity.We have avoided this by replacing elements of T B 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(a 2 s ) in the expressions of T jk .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 time-symmetrized 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 [41][42][43].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.

FIG. 1 .
FIG.1.(Left:) Illustration of discretized T00,m given 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,n but 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/2 to t − dt/2.

FIG. 3 .
FIG. 3. Relative error (61) in the conjectured Bianchi identity at different lattice sites in the transverse plane (A, B and C) as a function of the longitudinal direction

FIG. 5 .
FIG. 5. Relative error (68) 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. 7 .
FIG.7.Relative error (68) in the second continuity equation(34) as a function of lattice spacing a 2 s for discretized (63) and naive (continuum limit) electric and magnetic fields(14).Black dashed lines correspond to the a 2 s power-law.