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

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 over-occupied gluon fields, is typically modeled using classical Yang–Mills theory, and is often referred to as the Glasma [2,3,4]. Recently, the Glasma stage has received substantial attention [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22], especially in terms of its properties as a medium, described by transport coefficients [23,24,25,26,27,28,29,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,5,6, 10,11,12, 15,16,20, 23,24,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,34,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 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. (2) 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 \(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 \({\mathcal {O}}(a_s^2)\).

We believe that the idea and the algorithm presented here to derive a conserved \(T_{\mu \nu }\) (with high precision) can be utilized much more broadly. Nonconservation of \(T^{\mu \nu }\) 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 [36,37,38].

This paper is structured as follows. Section 2 describes our discretization framework. In Sect. 3 we construct the discrete energy–momentum tensor and introduce its time-symmetrized formulation. Section 4 shows 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(\(N_c\)) Yang–Mills theory in temporal axial gauge (\(A_t = 0\)),

$$\begin{aligned} {\mathcal {H}}&= \sum _{m}a_s^3\Bigg (-\frac{1}{2\sqrt{-\text {det}(g_{\mu \nu })}}\sum _{\mu >0}g_{\mu \mu }E^{\mu ,a}_mE^{\mu ,a}_m\nonumber \\&\quad + \frac{2N_c}{g^2} \sum _{0<\mu <\nu }\frac{1}{a_\mu ^2a_\nu ^2}g^{\mu \mu }g^{\nu \nu }\nonumber \\&\quad \times \sqrt{-\text {det}(g_{\mu \nu })}\Big (1-\frac{1}{N_c}\mathrm{{ReTr}}(U_{\mu \nu ,m})\Big )\Bigg ) \end{aligned}$$
(1)
$$\begin{aligned}&= \sum _{m}a_s^3\Bigg (\frac{1}{2}\sum _{\mu >0} E^{\mu ,a}_mE^{\mu ,a}_m \nonumber \\&\quad + \frac{2N_c}{g^2} \sum _{0<\mu <\nu }\frac{1}{a_\mu ^2a_\nu ^2}\Big (1-\frac{1}{N_c}\mathrm{{ReTr}}(U_{\mu \nu ,m})\Big )\Bigg ). \end{aligned}$$
(2)

Here m denotes the lattice site on the spatial grid and \(a=1, \ldots ,N_c^2-1\) is the color index. In this paper, we consider the Minkowski metric with \(g_{\mu \nu }=(+1,-1,-1,-1)\) and \(\sqrt{-\text {det}(g_{\mu \nu })} = 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 \times N_y \times N_z\), lattice spacing \(a_\mu \) in the spatial direction \({\hat{\mu }}\), and the product \(a_s^3 = a_x a_y a_z\). To guarantee gauge invariance, the theory is formulated in terms of gauge links \(U_{i,m}= e^{ia_i gA_{i,m}}\) instead of gauge fields \(A_{i,m}\). Links enter Eq. (2) in terms of plaquettes

$$\begin{aligned} U_{\mu \nu ,m}=U_{\mu ,m}U_{\nu ,m+{\hat{\mu }}}U^\dagger _{\mu ,m+{\hat{\nu }}}U^\dagger _{\nu ,m}, \end{aligned}$$
(3)

where \({{\hat{\mu }}}\) denotes a unit vector in the \(\mu \) direction.Footnote 1 In analogy to the links, plaquettes are related to the field-strength tensor \(F_{\mu \nu }\) via \(U_{\mu \nu } \approx e^{i g a_\mu a_\nu F_{\mu \nu }}\).

Links and plaquettes are group elements \(U_{\mu ,m}, U_{\mu \nu ,m} \in \) 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^{\mu }_m=E^{\mu ,a}_m t^a\). Using the (continuum) relation between the field-strength tensor and the chromomagnetic field \(B^i = -\frac{1}{2} \epsilon ^{ijk} F_{jk}\), the Hamiltonian (2) reduces to the correct expression \(\int d^3x \frac{1}{2}\sum _{j>0}(E^{j,a} E^{j,a} + B^{j,a} B^{j,a})\) in the continuum limit \(a_\mu \rightarrow 0\).

The Hamiltonian (2) is associated withFootnote 2 the Hamiltonian equations of motion for the link matrices and the electric field variables, which are discrete in space but continuous in time

$$\begin{aligned} \partial _0 U_{\mu ,m}&{}=iga_\mu E^{\mu ,a}_m t^a U_{\mu ,m} \end{aligned}$$
(4)
$$\begin{aligned} \partial _0 \,g E^{\mu ,a}_m&{}=\sum _{0<\nu \ne \mu }\frac{2 a_\mu }{a_\mu ^2a_\nu ^2}\textrm{ReTr}\big (it^aU_{\mu \nu ,m}+it^aU_{\mu -\nu ,m}\big ) \end{aligned}$$
(5)
$$\begin{aligned}&{}=\sum _{0<\nu \ne \mu }\frac{2}{a_\mu a_\nu }\textrm{ReTr}\big (it^a[D_\nu ^B,U_{\mu \nu ,m}]\big ). \end{aligned}$$
(6)

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_{\mu -\nu ,m}^\dagger = U_{\nu ,m-{{\hat{\nu }}}}^\dagger U_{\mu \nu ,m-{{\hat{\nu }}}}U_{\nu ,m-{{\hat{\nu }}}}\) and introduced the backward and forward gauge covariant derivatives

$$\begin{aligned} D_{\mu ,m}^BX={}&\big (X_{m}-U^\dagger _{\mu ,m-{\hat{\mu }}}X_{m-{\hat{\mu }}}~U_{\mu ,m-{\hat{\mu }}}\big )/a_\mu \end{aligned}$$
(7)
$$\begin{aligned} D_{\mu ,m}^FX={}&\big (U_{\mu ,m}X_{m+{\hat{\mu }}}U^\dagger _{\mu ,m}-X_{m}\big )/a_\mu . \end{aligned}$$
(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

$$\begin{aligned} \partial _0\mathrm{{ReTr}}\,U_{\mu \nu ,m}&{}=\mathrm{{ReTr}}\Big (ig a_\mu a_\nu \big [D_\mu ^F,E^\nu _m\big ]U_{\mu \nu ,m}\nonumber \\&\quad -iga_\mu a_\nu \big [D_\nu ^F,E^\mu _m\big ]U_{\mu \nu ,m} \Big ). \end{aligned}$$
(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^{\mu \nu }\) obtained by taking a functional derivative of the Yang–Mills action \(S = \sum _m a_s^3 \sqrt{-\text {det}(g_{\mu \nu })} {\mathcal {L}}\) with respect to the metric tensor \(g_{\mu \nu }\),

$$\begin{aligned} T^{\mu \nu } = \frac{2}{\sqrt{-\text {det}(g_{\mu \nu })}} \frac{\delta S}{\delta g_{\mu \nu }} = 2\frac{\partial {\mathcal {L}}}{\partial g_{\mu \nu }}-g^{\mu \nu }{\mathcal {L}}. \end{aligned}$$
(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 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 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

$$\begin{aligned} \partial _\mu T^{\mu \nu } =0. \end{aligned}$$
(11)

We will start in Sect. 3.2 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 \(\nu = 0\). The remaining components will then be constructed in Sect. 3.3 in the same way by taking \(T^{0i}\) as the starting point and requiring that \(T^{ij}\) satisfy the remaining continuity equations for \(\nu > 0\). This procedure determines \(T^{0i}\) up to a transformation \(T^{0i} \rightarrow T^{0i} + \phi ^i\), where \(\partial _i \phi ^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 E and B fields with their spatially averaged discretized counterparts

$$\begin{aligned} T^{00}_{m}&=\frac{1}{2}\big (E^2_{m,\textrm{loc}}+B^2_{m,\textrm{loc}}) \end{aligned}$$
(12)
$$\begin{aligned} T^{0i}_m&=\epsilon ^{ijk}\big (E^j_{m,\textrm{loc}} B^k_{m,\textrm{loc}}\big ) \end{aligned}$$
(13)
$$\begin{aligned} T^{ij}_m&=\frac{1}{2}\delta ^{ij}\big (E^2_{m,\textrm{loc}} + B^2_{m,\textrm{loc}}\big )\nonumber \\&\quad - E^i_{m,\textrm{loc}} E^j_{m,\textrm{loc}} - B^i_{m,\textrm{loc}} B^j_{m,\textrm{loc}} \nonumber \\&= \delta ^{ij}T^{00}_m - E^i_{m,\textrm{loc}} E^j_{m,\textrm{loc}} - B^i_{m,\textrm{loc}} B^j_{m,\textrm{loc}} \end{aligned}$$
(14)

where \(E^2_{m,\textrm{loc}}=E^{i,a}_{m,\textrm{loc}}E^{i,a}_{m,\textrm{loc}}\) and \(E^{i,a}_{m,\textrm{loc}}\) and \(B^{i,a}_{m,\textrm{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+{\hat{i}}\), it is rather centered at the point \(m+{\hat{i}}/2\). Similarly, a plaquette \(U_{ij,m}\) is actually centered at \(m+{\hat{i}}/2+{\hat{j}}/2\). Thus a natural way to construct electric and magnetic fields at a site m is to take nearest neighbor averages to obtain symmetric expressions:

$$\begin{aligned} E^{i,a}_{m,\textrm{loc}}&= \frac{1}{2}\big (E^{i,a}_m + U^\dagger _{i,m-\hat{i}}E^{i,a}_{m-\hat{i}}U_{i,m-\hat{i}}\big ) \end{aligned}$$
(15)
$$\begin{aligned} B^{i,a}_{m,\textrm{loc}}&= -\frac{\epsilon ^{ijk}}{8ga_ja_k} \textrm{ReTr}\nonumber \\&\quad \times \Big (it^a \big (U_{jk,m} + U_{j-k,m} + U_{-jk,m} + U_{-j-k,m}\big ) \Big ) \end{aligned}$$
(16)

3.2 Continuity equation for energy conservation

Our starting point will be the scalar component of the continuity equation

$$\begin{aligned} \partial _0 T_{00} = \partial _i T_{0i}, \end{aligned}$$
(17)

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 \(\mu =0\), \(\nu =0\) component of the energy–momentum tensor \(T_{\mu \nu }\). It is worth noting that in the Hamiltonian formulation, the electric field labeled \(E_{m}^i\) is located at the position \(m + \frac{\hat{i}}{2} + \frac{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 + \frac{\hat{i}}{2} + \frac{\hat{j}}{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

figure a
Fig. 1
figure 1

(Left:) Illustration of discretized \(T_{00,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 \(U_{ij,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\)

Full size image

Both the electric and magnetic field components of Eq. (18) lead to the same total energy as the HamiltonianFootnote 3 (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):

$$\begin{aligned}&\sum _{j,a}\left( E^{j,a}_m + U^\dagger _{j,m-j}E^{j,a}_{m-j}U_{j,m-j}\right) ^2 \nonumber \\&\quad \ne 2\sum _{j,a} \left[ E^{j,a}_mE^{j,a}_m+ E^{j,a}_{m-j}E^{j,a}_{m-j} \right] . \end{aligned}$$
(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 \(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

$$\begin{aligned} U_{-kj,m}&{} =(U_{j-k,m})^\dagger =U^\dagger _{k,m-k}U_{jk,m-k}U_{k,m-k} \end{aligned}$$
(20)
$$\begin{aligned} U_{k-j,m}&{} = (U_{-jk,m})^\dagger = U^\dagger _{j,m-j}U_{jk,m-j}U_{j,m-j} \end{aligned}$$
(21)
$$\begin{aligned} U_{-j-k,m}&{}= U^\dagger _{j,m-j}U^\dagger _{k,m-j-k}U_{jk,m-j-k}U_{k,m-j-k}U_{j,m-j} \end{aligned}$$
(22)

and the fact that the parallel transporting links cancel in the trace \(\textrm{Tr}(U^\dag M U) = \textrm{Tr}M\). This allows us to rewrite the energy density \(T_{00,m}\) as

$$\begin{aligned}&T_{00,m}= \sum _{j,k>0}\frac{1}{g^2a_j^2a_k^2}\Big [N_c-\frac{1}{4}\nonumber \\&\quad \times \Big (\textrm{ReTr}\big (U_{jk,m} + U_{jk,m-j} + U_{jk,m-k} + U_{jk,m-j-k}\big )\Big )\Big ]\nonumber \\&\quad +\frac{1}{4}\sum _{j>0}E^{j,a}_mE^{j,a}_m+ E^{j,a}_{m-j}E^{j,a}_{m-j}. \end{aligned}$$
(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 \(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):

$$\begin{aligned} \partial _0T_{00,m}&{}=\sum _{j,k, j\ne k}\frac{a_ja_k}{2ga_j^2a_k^2} \Big \{\textrm{ReTr}\big (i[D_k^F,E^j_m]U_{jk,m} \nonumber \\&\quad + i[D_k^F,E^j_{m-j}]U_{jk,m-j} + i[D_k^F,E^j_{m-k}]U_{jk,m-k}\nonumber \\&\quad {}+i[D_k^F,E^j_{m-j-k}]U_{jk,m-j-k}\big ) +2\textrm{ReTr}\big (iE^j_m\nonumber \\&\quad \times [D_k^B,U_{jk,m}]+iE^j_{m-j}[D_k^B,U_{jk,m-j}]\big )\Big \} \end{aligned}$$
(24)
$$\begin{aligned}&= \sum _{j,k, j\ne k}\frac{a_j}{2ga_j^2a_k^2}\Big \{\textrm{ReTr}\nonumber \\&\quad \times \big ({iU_{k,m}E^j_{m+k}U^\dagger _{k,m}U_{jk,m}} + {iE^j_mU_{jk,m}}\nonumber \\&\quad + {iU_{k,m-j}E^j_{m-j+k}U^\dagger _{k,m-j}U_{jk,m-j}} \nonumber \\&\quad + iE^j_{m-j}U_{jk,m-j}\nonumber \\&\quad - {iU_{k,m-k}E^j_mU^\dagger _{k,m-k}U_{jk,m-k}} \nonumber \\&\quad - {iE^j_{m-k}U_{jk,m-k}} \nonumber \\&\quad - {iU_{k,m-j-k}E^j_{m-j}U^\dagger _{k,m-j-k}U_{jk,m-j-k}} \nonumber \\&\quad - iE^j_{m-j-k}U_{jk,m-j-k}\big )\Big \} \end{aligned}$$
(25)

Indeed, it is possible to write the right-hand side as a total spatial derivative. After some rearrangement, one achieves the following result:

$$\begin{aligned}&\partial _0T_{00,m}^{(N)}=\sum _{j,k, k\ne j}\frac{1}{2ga_j a_k}\nonumber \\&\quad \times \Big \{\textrm{ReTr}\big [D_k^B,iU_{k,m}E^j_{m+k}U^\dagger _{k,m}U_{jk,m} + i E^j_mU_{jk,m} \nonumber \\&\quad \times +iU_{k,m-j}E^j_{m-j+k}U^\dagger _{k,m-j}U_{jk,m-j}+iE^j_{m-j}U_{jk,m-j}\big ]\Big \}. \end{aligned}$$
(26)

Using the definition of the covariant backward (or forward) derivative in Eqs. (7), (8), it is now easy to see that under the \(\textrm{ReTr}\) operation, this expression simplifies to an ordinary derivative of a scalar quantity

$$\begin{aligned} \textrm{ReTr} D_{\mu , m}^BX= & {} \frac{1}{a_\mu }\,\textrm{ReTr} \big (X_{m}-U^\dagger _{\mu , m-{\hat{\mu }}}X_{m-{\hat{\mu }}}~U_{\mu , m-{\hat{\mu }}}\big ) \nonumber \\= & {} \frac{1}{a_\mu }\,\textrm{ReTr} \big (X_{m}-X_{m-{\hat{\mu }}}\big ) = \partial _{\mu , m}^B \textrm{ReTr} X.\nonumber \\ \end{aligned}$$
(27)

Thus we have arrived at the first important result of this paper, a discrete energy conservation law

$$\begin{aligned} \partial _0T_{00,m} = \partial _{k, m}^B T_{0k,m}. \end{aligned}$$
(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:

figure b

where

$$\begin{aligned} c_{jk} = \frac{1}{2 g a_j a_k} \end{aligned}$$
(30)

It is important to note that what we label by \(T_{0k,m}\) is actually centered at the position \(m+\frac{\hat{k}}{2}\),Footnote 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 \(\partial ^B_{k, m} T_{0k}\), is in fact centered around m, the same position where \(\partial _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} \rightarrow \epsilon ^{kij} E^i B^j\). This can be seen by expanding the plaquettes to quadratic order in \(a_s\), thus \(U_{jk} \approx 1+ i g a_j a_k F_{jk}\), and realizing that the term with the identity vanishes because the matrix \(E^j\) is traceless.

3.2.2 Symmetric time discretization

Until now, we have treated t as a continuous variable, which is the limit \(dt/a_s \ll 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 \(\partial _0T_{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

$$\begin{aligned} \partial _0T_{00,m}^E&{}= \frac{1}{2}\bigg ( \frac{E^2_{m}\big (t-\frac{dt}{2}\big ) - E^2_{m}\big (t-\frac{3dt}{2}\big )}{dt}\bigg ) \nonumber \\&{}=\bigg (\frac{E^i_m\big (t-\frac{dt}{2}\big ) + E^i_m\big (t-\frac{3dt}{2}\big )}{2}\bigg )\nonumber \\&\quad \times \bigg (\frac{E^i_m\big (t-\frac{dt}{2}\big ) - E^i_m\big (t-\frac{3dt}{2}\big )}{dt}\bigg )\nonumber \\&{}=E^{i,\textrm{avg}}_m\big (t-dt\big )\partial _0E^{i}_m\left( t-\frac{dt}{2}\right) , \end{aligned}$$
(31)

In the leapfrog scheme, the time difference labeled \(\partial _0E^{i}_m\big (t-\frac{dt}{2}\big )\) 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:

$$\begin{aligned}&E^j_m[D_k^B,U_{jk,m}] \rightarrow \frac{1}{2}\left( E^j_m \left( t-\frac{dt}{2}\right) \right. \nonumber \\&\quad \left. + E^j_m \left( t-\frac{3dt}{2}\right) \right) \big [D_k^B,U_{jk,m}(t-dt)\big ]. \end{aligned}$$
(32)

Similarly, the time derivative of the magnetic field part of the energy density corresponds to the term in \(\partial _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 \(\partial _k T_{0k}\) as

$$\begin{aligned}&[D_k^F,E^j_m]U_{jk,m}\rightarrow \frac{1}{2}\left( \left[ D_k^{F}(t),E^j_m\left( t-\frac{dt}{2}\right) \right] \right. \nonumber \\&\quad \left. + \left[ D_k^{F}(t-dt),E^j_m\left( t-\frac{dt}{2}\right) \right] \right) \nonumber \\&\quad \times \frac{1}{2}\Big (U_{jk,m}(t) + U_{jk,m} (t-dt) \Big ), \end{aligned}$$
(33)

where the notation \(D_k^{F}(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 \(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\)

$$\begin{aligned} \partial _\mu T^{\mu i} = 0 \end{aligned}$$
(34)

to derive the \(T_{ij}\) components from \(\partial _0 T_{0i}\). We first observe that \(T_{0k}\) comprises two kinds of terms, with one of them being merely shifted in the \(-\hat{j}\) direction. Henceforth, we will focus solely on the first one of the two, which we call \(T_{0k}^1\) for our subsequent calculations. The second term can then be restored in the end by shifting the site of \(T_{0k}^1\). Let us first take the time derivative of Eq. (2) and split it into terms where the time derivatives act either on the electric fields or the plaquettes.

$$\begin{aligned}&\partial _0 T_{0k,m}^{1} = c_{jk} \nonumber \\&\quad \Bigg \{ \underbrace{E^{j,a}_{m+k}\partial _0\textrm{ReTr}(it^a U_{-kj,m+k})+ E^{j,a}_m \partial _0\textrm{ReTr} (it^aU_{jk,m})}_{\partial _0T_{0k}^{1E}} \nonumber \\&\quad {+} \underbrace{(\partial _0 E^{j,a}_m) \textrm{ReTr} (it^aU_{jk,m}) {+} (\partial _0E^{j,a}_{m{+}k})\textrm{ReTr}(it^a U_{-kj,m+k})}_{{\partial _0T_{0k}^{1B}}}\Bigg \}. \end{aligned}$$
(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 \(T_{ij}\), which we refer to here as \(\partial _0T_{0k}^{1E}\). 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 \(\partial _0T_{0k}^{1B}\).Footnote 5

Fig. 2
figure 2

Illustration of the time derivative of \(T_{0k}^{1}\) in 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 \(\partial _0T_{0k}^{1B}\) as defined in Eq. (53), while the subsequent eight terms on the following line correspond to the electric field component of \(\partial _0T_{0k}^{1E}\) as defined in Eq. (36)

Full size image

3.3.1 Chromoelectric field contribution

Let us start with the chromoelectric field contribution \(\partial _0T_{0k}^{1E}\). We begin by employing the equation of motion for the plaquette to express the electric field part \({\partial _0T_{0k}^{1E}}\) as

$$\begin{aligned} \partial _0T_{0k,m}^{1E}&= \sum _{j\ne k}gc_{jk}\textrm{ReTr}\nonumber \\&\quad \times \Big \{a_j\Big (-E^j_{m+k}U^\dagger _{k,m}E^j_mU_{k,m}U_{-kj,m+k} \nonumber \\&\quad + E^j_{m+k}E^j_{m+k}U_{-kj,m+k}\Big )\nonumber \\&\quad + a_k\Big (E^j_{m+k}U^\dagger _{k,m}E^k_mU_{k,m}U_{-kj,m+k} \nonumber \\&\quad - E^j_{m+k}U^\dagger _{k,m}U_{j,m}E^k_{m+j}U^\dagger _{j,m}U_{k,m}U_{-kj,m+k} \Big )\nonumber \\&\quad + a_j\Big (-E^j_mE^j_mU_{jk,m} + U_{k,m}E^j_{m+k}U^\dagger _{k,m}E^j_{m}U_{jk,m}\Big ) \nonumber \\&\quad {+} a_k\Big (\!-E^j_{m}U_{j,m}E^k_{m+j}U^\dagger _{j,m}U_{jk,m} {+} E^k_mE^j_mU_{jk,m}\Big ) \Big \}. \end{aligned}$$
(36)

Here, we have utilized the evolution equation of the plaquette \(U_{-kj,m}\), which can be derived in analogy to Eq. (9),

$$\begin{aligned}&\partial _0\textrm{ReTr}(it^aU_{-kj,m})\nonumber \\&\quad =ig\textrm{ReTr}\Big (it^a a_j\big (U^\dagger _{k,m-k}E^j_{m-k}U_{k,m-k}U_{-kj,m} - U_{-kj,m}E^j_m\big )\nonumber \\&\qquad - it^a a_k\big (U^\dagger _{k,m-k}E^k_{m-k}U_{k,m-k}U_{-kj,m} \nonumber \\&\qquad - U^\dagger _{k,m-k}U_{j,m-k}E^k_{m+j-k}U^\dagger _{j,m-k}U_{k,m-k}U_{-kj,m}\big ) \Big ). \end{aligned}$$
(37)

By employing the expression in Eq. (20), we can rewrite Eq. (36) as follows:

$$\begin{aligned} \partial _0T_{0k,m}^{1E}&= \sum _{j\ne k} gc_{jk}\textrm{ReTr}\nonumber \\&\quad \times \Big \{ a_j\Big ( E^j_{m+k}E^j_{m+k}U^\dagger _{k,m}U_{jk,m}U_{k,m} - E^j_mE^j_mU_{jk,m}\Big ) \nonumber \\&\quad +a_k\Big (E^j_{m+k}U^\dagger _{k,m}E^k_mU_{jk,m}U_{k,m} \nonumber \\&\quad - E^j_{m+k}U^\dagger _{k,m}U_{j,m}E^k_{m+j}U^\dagger _{j,m}U_{jk,m}U_{k,m}\nonumber \\&\quad - E^j_mU_{j,m}E^k_{m+j}U^\dagger _{j,m}U_{jk,m} + E^k_mE^j_mU_{jk,m} \Big )\Big \}. \nonumber \\ \end{aligned}$$
(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 \(T_{0i}\) from \(\partial _0T_{00}\) was straightforward, this part is more challenging. In Fig. 2, we illustrate the issue by examining the time derivative of \(T_{0k}^1\) (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 \(\partial _0T_{0k}^{1E}\). Notably, each term entering the latter includes an extra plaquette, as explicitly stated in Eq. (36). This differs from the continuum expression \(T_{jk}^{E} = E_jE_k + \frac{1}{2} \delta _{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 \(\partial _0T_{0k}^{1E}\) 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 \(\partial _0T_{0k}^{1E}\) with the identity matrix. This introduces a relative error \({\mathcal {O}}(a^2)\) and agrees with the original expression in the continuum limit. After rearranging certain terms and adding the contribution for \(\partial _0T_{0k}^{1E}|_{m\rightarrow {m-j}}\) (as in Eq. 29) the resulting expression takes the following form:

$$\begin{aligned} \partial _0T_{0k,m}^{E}&= \sum _{j\ne k} gc_{jk} \textrm{ReTr}\Big \{ a_j\Big (E^j_{m+k}E^j_{m+k} \nonumber \\&\quad - E^j_mE^j_m + E^j_{m+k-j}E^j_{m+k-j} - E^j_{m-j}E^j_{m-j}\Big )\nonumber \\&\quad -a_ka_j\Big (E^j_m\big [D_j^F,E^k_m\big ] \nonumber \\ {}&\quad + U_{k,m}E^j_{m+k}U^\dagger _{k,m}\big [D_j^F,E^k_m] + E^j_{m-j}\big [D_j^F,E^k_{m-j}\big ]\nonumber \\&\quad + U_{k,m-j}E^j_{m+k-j}U^\dagger _{k,m-j}\big [D_j^F,E^k_{m-j}] \Big ) \Big \} + {\mathcal {O}}(a^2) . \end{aligned}$$
(39)

It is important to highlight that when j is equal to k, different terms in the first and second lines of \(\partial _0T_{0k}^{E}\) cancel each other. This permits us to interchange the summations \(\sum _{j\ne k}\) and \(\sum _{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_{jk}^E\). The discretized form of the product rule is given as

$$\begin{aligned}&\big [D_k^{F}, E^i_mE^j_m\big ] \nonumber \\&\quad = \big [D_k^{F}, E^i_m\big ]U_{k,m}E^j_{m+k}U_{k,m} + E^i_m\big [D_k^{F}, E^j_m\big ], \end{aligned}$$
(40)

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

$$\begin{aligned}{} & {} \big [D_k^{F}, U_{k,m-k}^\dagger E^i_{m-k} U_{k,m-k}E^j_m\big ]\nonumber \\{} & {} \quad = [D_k^F, U^\dagger _{m-k}E^i_{m-k}U_{k,m-k}]E^j_m + E^i_m[D_k^F, E^j_m]. \end{aligned}$$
(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

$$\begin{aligned}&[D_j^B, U^\dagger _{j,m}E^j_mU_{j,m}E^k_{m+j}] = [D_j^F, U^\dagger _{j,m-j}E^j_{m-j}U_{j,m-j}E^k_m] \nonumber \\&\quad = E^j_mU_{j,m}E^k_{m+j}U^\dagger _{j,m} \nonumber \\&\qquad - E^j_mE^k_m - U^\dagger _{j,m-j}E^j_{m-j}U_{j,m-j}E^k_m + E^j_mE^k_m\nonumber \\&\qquad =E^j_{m}[D_j^F,E^k_m] + E^k_m\hspace{-0.2cm}\underbrace{[D_j^B,E^j_m]}_{=~0~ \mathrm{Gauss~Law}} \end{aligned}$$
(42)
$$\begin{aligned}&[D_j^B,U^\dagger _{j,m}U_{k,m}E^j_{m+k}U^\dagger _{k,m}U_{j,m}E^k_{m+j}] \nonumber \\&\quad = \big [D_j^F,U^\dagger _{j,m-j}U_{k,m-j}E^j_{m-j+k}U^\dagger _{k,m-j}U_{j,m-j}E^k_m\big ] \nonumber \\&\quad =U_{k,m}E^j_{m+k}U^\dagger _{k,m}U_{j,m}E^k_{m+j}U^\dagger _{j,m} - U_{k,m}E^j_{m+k}U^\dagger _{k,m}E^k_m\nonumber \\&\qquad - U_{k,m-j}E^j_{m+k-j}U^\dagger _{k,m-j}U_{j,m-j}E^k_mU^\dagger _{j,m-j} \nonumber \\&\qquad + U_{k,m}E^j_{m+k}U^\dagger _{k,m}E^k_m\nonumber \\&\quad = U_{k,m}E^j_{m+k}U^\dagger _{k,m}\big [D_j^F,E^k_m\big ] + \underline{\big [D_j^B,U_{k,m}E^j_{m+k}U^\dagger _{k,m}\big ]E^k_m} \end{aligned}$$
(43)
$$\begin{aligned}&\big [D_j^B,E^j_{m}E^k_{m}\big ] =\big [D_j^F,E^j_{m-j}E^k_{m-j}\big ]\nonumber \\&\quad = E^j_{m-j}[D_j^F,E^k_{m-j}] + E^k_m\hspace{-0.2cm}\underbrace{[D_j^B,E^j_m]}_{= ~0~ \mathrm{Gauss~Law}} \end{aligned}$$
(44)
$$\begin{aligned}&\big [D_j^B, U_{k,m}E^j_{m+k}U^\dagger _{k,m}E^k_{m}\big ]\nonumber \\&\quad =\big [D_j^F, U_{k,m-j}E^j_{m+k-j}U^\dagger _{k,m-j}E^k_{m-j}\big ]\nonumber \\&\quad = U_{k,m-j}E^j_{m+k-j}U^\dagger _{k,m-j}\big [D_j^F,E^k_{m-j}\big ] \nonumber \\&\qquad + \underline{\big [D_j^B,U_{k,m}E^j_{m+k}U^\dagger _{k,m}\big ]E^k_m} \end{aligned}$$
(45)

We can make slight modifications to the underlined terms in Eqs. (43) and (45) to make use of Gauss’ law.

$$\begin{aligned}&\underline{\big [D_j^B,U_{k,m}E^j_{m+k}U^\dagger _{k,m}\big ]E^k_m}\nonumber \\&\quad = U_{k,m}E^j_{m+k}U^\dagger _{k,m}E^k_m - U^\dagger _{j,m-j}\nonumber \\&\qquad \times U_{k,m-j}E^j_{m+k-j}U^\dagger _{k,m-j}U_{j,m-j}E^k_m\nonumber \\&\quad {}= U_{k,m}E^j_{m+k}U^\dagger _{k,m}E^k_m - U_{-jk,m}U_{k,m}U^j_{m-j+k}E^j_{m-j+k}\nonumber \\&\qquad \times U_{-kj,m-j+k}U_{j,m-j+k}U^\dagger _{k,m}E^k_m \nonumber \\&\quad {}\simeq U_{k,m}E^j_{m+k}U^\dagger _{k,m}E^k_m - \mathbb {1}U_{k,m}U^j_{m-j+k}E^j_{m-j+k}\nonumber \\&\qquad \times \mathbb {1}U_{j,m-j+k}U^\dagger _{k,m}E^k_m \nonumber \\&\quad {}=U_{k,m}[D^B_j,E^j_{m+k}]U^\dagger _{k,m}E^k_m \nonumber \\&\quad {}=0 \hspace{1cm}\mathrm{[Gauss~Law]} \end{aligned}$$
(46)

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 \({\mathcal {O}}(a^2)\) error, which enables us to apply Gauss’ law. By employing Eqs. (42)–(45), we can rephrase Eq. (39) as follows:

$$\begin{aligned} \partial _0T_{0k,m}^{E}&= g \sum _j\sum _lc_{lj}\delta _{jk}a_la_j \textrm{ReTr}\nonumber \\&\quad \times \Big \{\big [D_j^B, E^l_{m+j}E^l_{m+j} + E^l_{m+j-l}E^l_{m+j-l}\big ]\Big \}\nonumber \\&\quad - g \sum _jc_{jk}a_k a_j \textrm{ReTr}\Big \{ \big [D_j^B, E^j_{m}E^k_{m} \nonumber \\&\quad + U^\dagger _{j,m}U_{k,m}E^j_{m+k }U^\dagger _{k,m}U_{j,m}E^k_{m+j} \nonumber \\&\quad + U^\dagger _{j,m}E^j_{m}U_{j,m}E^k_{m+j} + U_{k,m}E^j_{m+k}U^\dagger _{k,m}E^k_{m}\big ]\Big \}\nonumber \\&\quad + {\mathcal {O}}(a^2). \end{aligned}$$
(47)

This allows us to use the second continuity equation (34), to identify the electric field component of the stress tensor as follows:

$$\begin{aligned} T_{jk,m}^{E}&= g\bigg \{\sum _l c_{lj}a_j a_l\delta _{jk} \textrm{ReTr}\nonumber \\&\quad \times \Big [E^l_{m+j}E^l_{m+j} + E^l_{m+j-l}E^l_{m+j-l}\Big ]\nonumber \\&\quad {}-a_ja_kc_{jk}\textrm{ReTr}\Big [U^\dagger _{j,m}U_{k,m}E^j_{m+k }U^\dagger _{k,m}U_{j,m}E^k_{m+j}\nonumber \\&\quad + E^j_{m}E^k_{m}+U^\dagger _{j,m}E^j_{m}U_{j,m}E^k_{m+j}\nonumber \\&\quad + U_{k,m}E^j_{m+k}U^\dagger _{k,m}E^k_{m}\Big ]\bigg \}. \end{aligned}$$
(48)

The contribution \(T_{jk}^{E}\) is located at the positionFootnote 6\(m+\frac{j}{2}+\frac{k}{2}\), precisely as expected from \(\partial _{0}T_{0k}\) calculations.

3.3.2 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 \(T_{0i} = E^jF_{ij}\)Footnote 7

$$\begin{aligned} \partial _0 T_{0i}^{B}&{}=[D_j,F_{jq}]F_{iq}. \end{aligned}$$
(49)

Comparing this expression to the spatial derivative of the spacelike part of the energy–momentum tensor \(T_{ij}^B = \frac{\delta _{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_{ij}^B = \frac{\delta _{ij}}{4}F_{pq}F_{pq}-F_{iq}F_{jq}\) and differentiate it:

$$\begin{aligned} \partial _j T_{ij}^{B}&= [D_j,T_{ij}^{B}] =-\frac{\delta _{ij}}{2}\nonumber \\&\qquad \times [D_j,F_{pq}]F_{pq}+[D_j,F_{iq}]F_{jq}+F_{iq}[D_j,F_{jq}]\nonumber \\&\quad =-\delta _{ij}[D_j,-\frac{1}{2}\epsilon ^{pqr}F_{pq}]B^r\nonumber \\&\qquad +\delta _{ij}[D_j,B^r]B^r-[D_j,B^j]B^i+F_{iq}[D_j,F_{jq}]\nonumber \\&\quad =[D_j,F_{jq}]F_{iq} \end{aligned}$$
(50)

In obtaining the last equality, we have utilized the relation \(B^r = -\frac{1}{2}\epsilon ^{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):

$$\begin{aligned} \frac{\delta _{ij}}{2}[D_j,F_{pq}]F_{pq} = [D_j,F_{iq}]F_{jq} \end{aligned}$$
(51)

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 [41]. 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,

$$\begin{aligned}&\frac{\delta _{ij}}{2}\textrm{ReTr}\Big (it^a\big [D_j^B,U_{pq,m}\big ]\Big )\textrm{ReTr}\Big (it^aU_{pq,m}\Big )\nonumber \\&\quad \approx \textrm{ReTr}\Big (it^a\big [D_j^B,U_{iq,m}\big ]\Big )\textrm{ReTr}\Big (it^aU_{jq,m}\Big ) . \end{aligned}$$
(52)

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 \(\partial _0 T_{0k}^{1B}\) 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

$$\begin{aligned} \partial _0 T_{0k,m}^{1B}&= \sum _{j\ne k} c_{jk}\Big \{\partial _0E^{j,a}_{m+k}\textrm{ReTr}(it^a U_{-kj,m+k})\nonumber \\&\quad + \partial _0 E^{j,a}_m \textrm{ReTr} (it^aU_{jk,m}) \Big \}\nonumber \\&= \sum _{\begin{array}{c} j\ne k\\ j\ne l \end{array}}c_{jk}d_{jl}\Big \{\textrm{ReTr}\big (it^a (U_{jl,m+k} + U_{j-l,m+k})\big )\textrm{ReTr}\nonumber \\&\quad \times \big (it^aU_{-kj,m+k}\big ) \nonumber \\&\quad + \textrm{ReTr}\big (it^a(U_{jl,m}+U_{j-l,m})\big )\textrm{ReTr}\big (it^aU_{jk,m}\big ) \Big \}\nonumber \\&=\sum _{l,j}c_{lk}d_{lj}a_j \nonumber \\&\quad \times \Big \{\underbrace{\textrm{ReTr}\big (it^a \big [D_j^B,U_{lj,m+k}\big ] \big ) \textrm{ReTr} \big (it^aU_{-kl,m+k}\big )}_{(A)} \nonumber \\&\quad + \underbrace{\textrm{ReTr}\big (it^a\big [D_j^B,U_{lj,m}\big ]\big ) \textrm{ReTr}\big (it^aU_{lk,m}\big )}_{(B)}\Big \}, \end{aligned}$$
(53)

where

$$\begin{aligned} d_{jl}=\frac{2a_j}{ga_j^2a_l^2}. \end{aligned}$$
(54)

Transitioning from the second equality to the third, we took into account \(\textrm{ReTr}(it^a\mathbb {1}) = 0\) under the conditions \(j =k\) and \(j=l\), and subsequently exchanged the indices j and l.

By employing the Fierz identity \(t^a_{ij}t^a_{kl} = \frac{1}{2}\big (\delta _{il}\delta _{jk}- \frac{1}{N_c}\delta _{ij}\delta _{kl}\big )\) for the generators in the fundamental representation, the multiplication of two arbitrary plaquettes can be streamlined. Let us at this point simplify our notations and leave out contributions arising from the \(N_c\)-suppressed term in the Fierz identity. For SU(2) this term is in fact zero because the trace of a unitary matrix is purely real. For SU(3) the trace of a plaquette can have an imaginary part, but it is suppressed by powers of the lattice spacing. Since we are in any case not obtaining an expression that is exact at all orders in the lattice spacing, we will not write these trace terms here. With these approximations, we obtain

$$\begin{aligned}&\textrm{ReTr}(it^aU_{\mu \nu ,m})\textrm{ReTr}(it^aU_{\alpha \beta ,m}) \nonumber \\&{} \approx \frac{i^2}{4}\textrm{ReTr}\big (U_{\mu \nu ,m}U_{\alpha \beta ,m} -U_{\mu \nu ,m}U_{\beta \alpha ,m}\big ), \end{aligned}$$
(55)

where we have utilized the fact that the anti-Hermitian part of the SU(2) matrix is traceless to remove the \(1/N_c\) term from the equation. Now, with the approximation (55), part (A) of Eq. (53) can be rewritten as

$$\begin{aligned} (A)&{} =\frac{i^2}{4}\textrm{ReTr}\big (U_{-kl,m+k} \big [D_j^B, U_{lj,m+k}-U_{jl,m+k}\big ]\big )\nonumber \\&{} = \frac{i^2}{4}\Big \{ \textrm{ReTr}\big (U_{-kl,m+k}\big [D_j^B,U_{lj,m+k}\big ] \nonumber \\&\quad + \big [D_j^B, U_{j,m+k}U_{-kl,m+k+j}U^\dagger _{j,m+k}\big ]U_{lj,m+k} \nonumber \\&\quad {}-\big [D_j^B, U_{j,m+k}U_{-kl,m+k+j}U^\dagger _{j,m+k}\big ]U_{lj,m+k}\nonumber \\&\quad - U_{-kl,m+k}\big [D_j^B,U_{jl,m+k}\big ] \nonumber \\&\quad {}- \big [D_j^B, U_{j,m+k}U_{-kl,m+k+j}U^\dagger _{j,m+k}\big ]U_{jl,m+k} \nonumber \\&\quad + \big [D_j^B, U_{j,m+k}U_{-kl,m+k+j}U^\dagger _{j,m+k}\big ]U_{jl,m+k}\Big \}. \end{aligned}$$
(56)

Here, we employed insights from our prior discussion on the continuum case to add and subtract terms strategically. By utilizing the discretized product rule, certain terms can be combined

$$\begin{aligned} (A)&{} =\frac{i^2}{4}\Big \{\textrm{ReTr}\Big (\big [D_j^B, U_{j,m+k}U_{-kl,m+k+j}U^\dagger _{j,m+k}U_{lj,m+k}\big ]\nonumber \\&\quad -\big [D_j^B, U_{j,m+k}U_{-kl,m+j+k}U^\dagger _{j,m+k}U_{jl,m+k}\big ]\nonumber \\&\quad {}-\big [D_j^B,U_{j,m+k}U_{-kl,m+k+j}U^\dagger _{j,m+k}\big ] U_{lj,m+k} \nonumber \\&\quad +\big [D_j^b,U_{j,m+k}U_{-kl,m+k+j}U^\dagger _{j,m+k}\big ]U_{jl,m+k} \Big )\Big \}\nonumber \\&{}=\frac{i^2}{4}\textrm{ReTr}\Big (\big [D_j^B,U_{j,m+k}U_{-kl,m+j+k}U^\dagger _{j,m+k}\nonumber \\&\quad \times \big (U_{lj,m+k}-U_{jl,m+k}\big )\big ]\nonumber \\&\quad {}-\big [D_j^B, U_{j,m+k}U^\dagger _{k,m+j}U_{lk,m+j}U_{k,m+j}U^\dagger _{j,m+k}\big ]\nonumber \\&\quad \times \big (U_{lj,m+k}-U_{jl,m+k}\big )\Big ), \end{aligned}$$
(57)

while the remaining terms can be substituted with the approximate lattice Bianchi Identity in Eq. (52) to give the final expression as

$$\begin{aligned}&(A) \approx \frac{i^2}{4}\textrm{ReTr}\Big (\big [D_j^B, U_{j,m+k}U_{-kl,m+k+j}U^\dagger _{j,m+k}\nonumber \\&\quad \times \big (U_{lj,m+k}-U_{jl,m+k}\big )\big ] \nonumber \\&\quad -\frac{\delta _{jk}}{4}\big [D_j^B, U_{pl,m+k}\big (U_{pl,m+k}-U_{lp,m+k}\big )\big ]\Big ). \end{aligned}$$
(58)

With this, we employ a similar approach for the term labeled as (B) in Eq. (53)

$$\begin{aligned} (B)&= \frac{i^2}{4}\textrm{ReTr}\Big (U_{lk,m}\big [D_j^B, U_{lj,m}-U_{jl,m}\big ]\Big )\nonumber \\&{}=\frac{i^2}{4}\textrm{ReTr}\Big (U_{lk,m}\big [D_j^B, U_{lj, m}\big ] + U^\dagger _{j,m-j}U_{lj,m-j}U_{j,m-j}\nonumber \\&\quad \times \big [D_j^B, U_{lk,m}\big ] - U^\dagger _{j,m-j}U_{lj,m-j}U_{j,m-j}\big [D_j^B, U_{lk,m}\big ] \nonumber \\&\quad {}- U_{lk,m}\big [D_j^B, U_{jl,m}\big ] \nonumber \\&\quad - U^\dagger _{j,m-j}U_{jl,m-j}U_{j,m-j}\big [D_j^B,U_{lk,m}\big ]\nonumber \\&\quad + U^\dagger _{j,m-j}U_{jl,m-j}U_{j,m-j}\big [D_j^B,U_{lk,m}\big ]\Big )\nonumber \\&{}=\frac{i^2}{4}\textrm{ReTr}\Big (\big [D_j^B, U_{lk,m}U_{lj,m} - U_{lk,m}U_{jl,m}\big ] \nonumber \\&\quad + \big [D_j^B,U_{lk,m}\big ] \times \big (U_{l-j,m} - U_{-jl,m}\big )\Big ) \nonumber \\&{}\approx \frac{i^2}{4}\textrm{ReTr}\Big (\big [D_j^B, U_{lk,m}\big (U_{lj,m} - U_{jl,m}\big )\big ]\nonumber \\&\quad -\frac{\delta _{jk}}{4} \big [D_j^B, U_{pl,m}\big (U_{pl,m}-U_{lp,m}\big )\big ]\Big ) \Big ). \end{aligned}$$
(59)

By combining the aforementioned equations (58) and (59), along with their shifted counterparts \(\partial _0 T_{0k}^{1B}|_{m\rightarrow {m-j}}\) in Eq. (2), we derive the approximate magnetic contribution of the stress tensor as

$$\begin{aligned}&T_{jk,m,CB}^{B} \approx \sum _{l}c_{lk}d_{lj}\frac{a_j}{4}\textrm{ReTr}\nonumber \\&\quad \times \bigg \{-U_{j,m+k}U^\dagger _{k,m+j}U_{lk,m+j}U_{k,m+j}U_{j,m+k}^\dagger \big (U_{lj,m+k}\nonumber \\&\quad - U_{jl,m+k}\big ) - U_{lk,m}\big (U_{lj,m}-U_{jl,m}\big ) +\frac{\delta _{jk}}{4}\nonumber \\&\quad \Big (U_{pl,m+k}\big (U_{pl,m+k} - U_{lp,m+k}\big ) \nonumber \\&\quad + U_{pl,m}\big (U_{pl,m}-U_{lp,m}\big )\Big )\nonumber \\&\quad -U_{j,m+k-l}U^\dagger _{k,m+j-l}U_{lk,m+j-l}U_{k,m+j-l}\nonumber \\&\quad \times \big (U_{lj,m+k-l}- U_{jl,m+k-l}\big ) - U_{lk,m-l}\big (U_{lj,m-l}-U_{jl,m-l}\big )\nonumber \\&\quad +\frac{\delta _{jk}}{4}\Big (U_{pl,m+k-l}\big (U_{pl,m+k-l} - U_{lp,m+k-l}\big ) \nonumber \\&\quad + U_{pl,m-l}\big (U_{pl,m-l}-U_{lp,m-l}\big )\Big ) \bigg \}. \end{aligned}$$
(60)

Here, the notation \(T_{jk, m, CB}^B\) 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 \(\delta _{jk}\)-part of \(T_{jk}^B\) 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).

Fig. 3
figure 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

Full size image

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 Sect. 4. The result is demonstrated in Fig. 3. Here we show the relative error in the lattice counterpart of the conjectured Bianchi identity, calculated as

$$\begin{aligned} \frac{\textrm{ReTr}\big (it^a\big [D_j^B,U_{lk,m}\big ]\big )\textrm{ReTr}\big (it^aU_{lj,m}\big ) - \frac{\delta _{jk}}{2}\textrm{ReTr}\big (it^a\big [D_j^B,U_{pl,m}\big ]\Big )\textrm{ReTr}\big (it^aU_{pl,m}\big )}{\frac{1}{V}\sum _{x}\textrm{ReTr}\big (it^a\big [D_j^B,U_{lk,m}\big ]\big )\textrm{ReTr}\big (it^aU_{lj,m}\big )} \end{aligned}$$
(61)

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_{jk}^B\), the Bianchi identity has been used only to manipulate the total energy contribution that is proportional to \(\delta _{jk}\). In the continuum, the magnetic field part of the \(\delta _{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_{00}^B\) and use it to get the magnetic \(\delta _{jk}\) part of \(T_{jk}\). This leads to

$$\begin{aligned} T_{jk,m}^{B}&{}= \sum _{l}c_{lk}d_{lj}\frac{a_j}{4}\textrm{ReTr}\nonumber \\&\quad \times \bigg \{-U_{j,m+k}U^\dagger _{k,m+j}U_{lk,m+j}U_{k,m+j}\nonumber \\&\quad U_{j,m+k}^\dagger \big (U_{lj,m+k}- U_{jl,m+k}\big )\nonumber \\&\quad - U_{lk,m}\big (U_{lj,m}-U_{jl,m}\big )\nonumber \\&\quad - U_{j,m+k-l}U^\dagger _{k,m+j-l}U_{lk,m+j-l}U_{k,m+j-l}\big (U_{lj,m+k-l}\nonumber \\&\quad - U_{jl,m+k-l}\big ) - U_{lk,m-l}\big (U_{lj,m-l}-U_{jl,m-l}\big )\bigg \} \nonumber \\&\quad - \delta _{jk}T_{00, m+j}^B. \end{aligned}$$
(62)

This is achieved by subtracting \(T_{00,m+j}^B\) based on the first term of the continuum expression \(T_{ij} = -1/4 \delta _{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 \(\partial _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.

3.3.3 Final formula

With this, we write our final expression for \(T_{jk}\) by unifying the electric and magnetic components (48) and (62)

figure c
Fig. 4
figure 4

The relative error (71) in the energy conservation equation is shown as a function of time on the left for \(N_s=32\) and on the right for \(N_s=128\). The lines correspond to different formulations of the Poynting vector using the ‘naive’ continuum (13), ‘discrete’ (2) and time-symmetrized (obtained using the procedure described by (31) and (33)) expressions. The upper panels represent \(dt/a_s\) values of 0.01 and 0.1 for (\(N_s=32\) and \(Q_sa_s = 0.25\) ) and for (\(N_s=128\) and \(Q_s a_s = 0.0625\)), respectively, while the lower panels employ 0.001 and 0.01, respectively

Full size image

3.4 Final expression

For readers seeking a concise overview, we summarize the approximation and provide the ultimate expression of the energy–momentum tensor here. The determination of \(T_{00}\) at a specific lattice site involves averaging the incoming and outgoing electric fields and considering the various orientations of plaquettes within the Hamiltonian density:

$$\begin{aligned} T_{00,m}&{}= \sum _{j,k>0}\frac{1}{g^2a_j^2a_k^2}\Big [N_c-\frac{1}{4}\Big (\textrm{ReTr} U_{jk,m} \nonumber \\&\quad + \textrm{ReTr} U_{j-k,m} + \textrm{ReTr} U_{-jk,m}+ \textrm{ReTr} U_{-j-k,m}\Big )\Big ]\nonumber \\&\quad +\frac{1}{4}\sum _{j>0} \Big (E^{j,a}_mE^{j,a}_m+ E^{j,a}_{m-j}E^{j,a}_{m-j}\Big ) \end{aligned}$$
(64)

Subsequently, by examining the equation of motion of the plaquette and electric fields, along with the energy conservation equation, we derive the Poynting vectors \(T_{0i}\):

$$\begin{aligned} T_{0k,m}&=\sum _{j>0,j\ne k}c_{jk}\Big \{\textrm{ReTr}\big (iE^j_{m+k}U_{-kj,m+k}\nonumber \\&\quad + i E^j_mU_{jk,m}+iE^j_{m-j+k}U_{-kj,m-j+k}\nonumber \\&\quad +iE^j_{m-j}U_{jk,m-j}\big )\Big \} \end{aligned}$$
(65)

Following a similar methodology, we utilize the momentum conservation equation to construct \(T_{jk}\), which presents certain challenges. Notably, terms such as \(E^{k} E^{j} U_{jk}\) emerge in the electric field contribution of \(\partial _0 T_{0k}\), which are absent in the continuum. Since these terms are suppressed by the lattice spacing, we substitute the plaquette with the identity matrix. Additionally, for the magnetic field contribution, the lack of a suitable Bianchi identity poses difficulties in obtaining the spatial derivative of \(T_{jk}^{B}\) from \(\partial _0 T_{0k}\). To address this issue, we replace the term proportional to \(\delta _{jk}\) in \(T_{jk}^{B}\) with magnetic components of the discrete energy density, resulting in a notable reduction in violation. The final expression for \(T_{jk}\) is as follows:

$$\begin{aligned} T_{jk,m}&= \bigg \{\sum _l gc_{lj}a_j a_l\delta _{jk} \textrm{ReTr}\nonumber \\&\quad \Big [E^l_{m+j}E^l_{m+j} + E^l_{m+j-l}E^l_{m+j-l}\Big ]\nonumber \\&\quad -ga_ja_kc_{jk}\textrm{ReTr}\Big [U^\dagger _{j,m}U_{k,m}E^j_{m+k }U^\dagger _{k,m}U_{j,m}E^k_{m+j}\nonumber \\&\quad + E^j_{m}E^k_{m}+U^\dagger _{j,m}E^j_{m}U_{j,m}E^k_{m+j} + U_{k,m}E^j_{m+k}U^\dagger _{k,m}E^k_{m}\Big ]\bigg \}\nonumber \\&\quad +\sum _{l}c_{lk}d_{lj}\frac{a_j}{4}\textrm{ReTr}\bigg \{-U_{-kl,m+j+k}\nonumber \\&\quad \big (U_{-jl,m+j+k}- U_{l-j,m+j+k}\big )\nonumber \\&\quad - U_{lk,m}\big (U_{lj,m}-U_{jl,m}\big )\nonumber \\&\quad - U_{-kl,m+j+k-l}\big (U_{-jl,m+j+k-l}- U_{l-j,m+j+k-l}\big ) \nonumber \\&\quad - U_{lk,m-l}\big (U_{lj,m-l}-U_{jl,m-l}\big )\bigg \} - \delta _{jk}T_{00, m+j}^B . \end{aligned}$$
(66)

where

$$\begin{aligned} c_{jk} = \frac{1}{2 g a_j a_k} \quad \quad \textrm{and} \quad \quad d_{jl}=\frac{2a_j}{ga_j^2a_l^2}. \end{aligned}$$
(67)

4 Numerical results

4.1 Initial conditions and numerical setup

To quantitatively assess our expressions for the energy–momentum tensor, we test the conservation laws in a classical Yang–Mills simulation. In our numerical simulations, we solve the evolution equations (4) and (6) for the lattice fields using the standard leapfrog algorithm, and then calculate components of the energy–momentum tensor. At the beginning, we take the fields from initial conditions of the form

$$\begin{aligned} \langle A_i^a({\textbf{p}},t{=}0) A_i^{*b}({\textbf{p}},t{=}0) \rangle&= 2 a_s^3 N_s^3 \delta _{ab}\,\frac{0.2}{g^2}\,\frac{Q_s}{p^2}\,e^{-p^2/2Q_s^2} \end{aligned}$$
(68)
$$\begin{aligned} \langle A_i^a(t=0, p) \rangle&= 0 \end{aligned}$$
(69)
$$\begin{aligned} E_i^a(t=0)&= 0, \end{aligned}$$
(70)

where g is the coupling constant, \(Q_s\) a dimensionful initial momentum scale, and \(a_s^3 N_s^3\) the lattice volume. The expectation value \(\langle \cdot \rangle \) implies an average over classical configurations of the lattice fields. In Eq. (68) we only initialize transverse modes, \(p^i A_i({\textbf{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 [23, 42]. 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.Footnote 8 After the characteristic timescale \(t\sim 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(N_c=2)\) gauge group instead of the physical SU(3) gauge group of QCD.

Fig. 5
figure 5

Relative error (72) in the momentum conservation equation as a function of time for \(N_s=32\) on the left and \(N_s=256\) on the right. Both are shown for the electric and magnetic field contributions in the continuum (14) and discretized (63) formulations

Full size image

4.2 Results

Based on the above expressions (18), (2) and (63) of the energy–momentum tensor, we will now study the violation of the continuity equations \(\partial _\mu T^{\mu \nu } = 0\). We start our analysis with the equation for energy conservation, where we quantify the deviation in the form of a relative error

$$\begin{aligned} \sqrt{\frac{\sum \nolimits _{m} \big (\partial _0 T_{00,m} -\partial _{i}T_{0i,m}\big )^2}{\sum \nolimits _{m} \big (\partial _0 T_{00,m}\big )^2}}\,. \end{aligned}$$
(71)

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. (2) or the time-symmetrized discretization outlined in Sect.  3.2.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

$$\begin{aligned} \sqrt{\frac{\sum \nolimits _{m} \big (\partial _0 T_{0i,m}^{E/B} -\partial _{j}T_{ij,m}^{E/B}\big )^2}{\sum \nolimits _{m} \big (\partial _0 T_{0i,m}^{E/B}\big )^2}} \end{aligned}$$
(72)

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.  (2) 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.

Fig. 6
figure 6

Measuring the violation (72) in the second continuity equation (34) by utilizing \(T_{ij}^B\) derived from the energy density \(T_{00}\) (63) and the conjectured Bianchi identity (60), considering two distinct lattice sizes: \(N_s=32\) and \(N_s=128.\)

Full size image
Fig. 7
figure 7

Relative error (72) in the second continuity equation (34) as a function of lattice spacing \(a_s^2\) for discretized (63) and naive (continuum limit) electric and magnetic fields (14). Black dashed lines correspond to the \(a_s^2\) power-law

Full size image

As discussed in Sect. 3.3.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 \(\delta _{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 (72). 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.

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 \(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_s^2\). The black dashed line corresponds to a power law \(\sim a_s^2\). We observe that in the limit of small lattice spacing, the relative error of the naive expression goes to zero slower than \(a_s^2\). In contrast to this, the improved discrete expressions follow the \(a_s^2\) power law very closely, thus indicating that they are \({\mathcal {O}}(a_s^2)\).

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 \({\mathcal {O}}(a_s^2)\) 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 \ll 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 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), (2), and (63), that are derived using classical field equations of motion in conjunction with the energy–momentum conservation law for \(T^{\mu \nu }\).

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^jE^jB^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_{jk}^B\) that are proportional to \(\delta _{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 \({\mathcal {O}}(a_s^2)\) 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 [42,43,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.