Implicit schemes for real-time lattice gauge theory

We develop new gauge-covariant implicit numerical schemes for classical real-time lattice gauge theory. A new semi-implicit scheme is used to cure a numerical instability encountered in three-dimensional classical Yang-Mills simulations of heavy-ion collisions by allowing for wave propagation along one lattice direction free of numerical dispersion. We show that the scheme is gauge covariant and that the Gauss constraint is conserved even for large time steps.


Introduction
Color Glass Condensate (CGC) effective theory [1] applies classical Yang-Mills theory to the area of high energy heavyion collisions. In the CGC description high energy nuclei can be treated as thin sheets of fast moving color charge which generate a classical gluon field. The collision of two such sheets produces the Glasma [2], which behaves classically at the earliest stages of the collision. Due to classical Yang-Mills theory being non-linear and the non-perturbative nature of the CGC, computer simulations are commonly used to investigate the time evolution of such systems [3][4][5][6][7]. Here, real-time lattice gauge theory provides a numerical treatment of classical Yang-Mills theory while retaining exact gauge invariance on the lattice. To name a few applications besides CGC and Glasma simulations, real-time lattice gauge theory is widely used for non-Abelian plasma simulations [8,9] and hard thermal loop (HTL) simulations [10][11][12], classical statistical simulations of fermion production [13,14], in studying sphalerons (in electroweak theory) [15,16], for determining the plasmon mass scale in non-Abelian gauge theory [17,18] or when studying perturbations on top of a non-Abelian background field [19].
The real-time lattice gauge theory approach is based on the discretization of Yang-Mills fields on a lattice in terms of a e-mail: ipp@hep.itp.tuwien.ac.at b e-mail: david.mueller@tuwien.ac.at so-called gauge links, i.e. Wilson lines connecting adjacent lattice sites. Using gauge link variables one can discretize the continuum Yang-Mills action in various ways, the simplest of which is the Wilson gauge action [20]. Varying this action with respect to the link variables one obtains discretized classical field equations, which are of the explicit leapfrog type. Moreover, in addition to the equations of motion, one also obtains the Gauss constraint, which is exactly conserved by the leapfrog scheme, even for finite time steps. Going further, the accuracy of the numerical approximation can be systematically improved by adding higher order terms to the standard Wilson gauge action [21].
In previous publications [22,23] we developed lattice Yang-Mills simulations for genuinely three-dimensional heavy-ion collisions in the CGC framework. Unlike the usual boost-invariant approach, we consider collisions of nuclei with thin, but non-vanishing support along the longitudinal direction (the beam axis) and simulate them in the laboratory frame. This has enabled us to study the effects of finite nuclear longitudinal extent (which is inversely proportional to the Lorentz gamma factor γ ) on the rapidity profile of the produced Glasma after the collision. The numerical scheme in these simulations is based on the standard Wilson gauge action with the fields coupled to external color currents. The treatment of these color charges is closely related to the colored particle-in-cell (CPIC) method [24,25], which is a non-Abelian extension of the particle-in-cell (PIC) method [26] commonly used in (Abelian) plasma simulations.
Unfortunately, these simulations suffer from a numerical instability that leads to an artificial increase of total energy if the lattice resolution is too coarse: even a single nucleus propagating along the beam axis, which should remain static and stable, eventually becomes unstable. Improving the resolution simply postpones the problem at the cost of much higher computational resources. We realize that this instability is due to numerical dispersion on the lattice inherent to the leapfrog scheme, which renders the dispersion relation of plane waves non-linear. As a consequence of numerical dis-persion, high frequency plane waves exhibit a phase velocity that is noticeably less than the speed of light on the lattice. The shape of the pulse of color fields is lost over time. At the same time, the color current "driving" the nucleus forward will not disperse by construction: the point-like color charges making up the current are simply moved from one cell to the next as the simulation progresses. Thus the shape of the current is always kept intact. This mismatch and the resulting instability is therefore related to the numerical Cherenkov instability [27], which can occur in (Abelian) particle-in-cell simulations. Notably, simulations of laser wakefield acceleration [28], where electric charges moving at relativistic speeds are coupled to discretized electromagnetic fields, suffer from the same type of instability and many numerical schemes have been devised to cure it [29][30][31][32]. A particularly simple solution to the problem is the use of semi-implicit schemes to repair the dispersion relation (i.e. making it linear) for one direction of propagation [33], which is the approach we take in this work.
In this paper we derive an implicit and a semi-implicit scheme for real-time lattice gauge theory by modifying the standard Wilson gauge action. We obtain two new actions that are gauge invariant (in the lattice sense), are of the same order of accuracy as the original action, but yield an implicit or semi-implicit scheme upon variation. In the case of the semi-implicit scheme setting the lattice spacing and the time step to specific values can fix the dispersion relation along the longitudinal direction and thus suppress the numerical Cherenkov instability. We also obtain a modified version of the Gauss constraint that is conserved up to (in principle) arbitrary numerical precision under the discrete equations of motion.
We start with a discussion of the main ideas behind the semi-implicit scheme for the two-dimensional wave equation in Sect. 2 and for Abelian gauge fields on the lattice in Sect. 3. The concepts are then generalized to non-Abelian lattice gauge theory in Sect. 4, where we derive both a fully implicit and the semi-implicit scheme. Finally, we verify numerically that the Cherenkov instability can be suppressed using the new scheme and that the Gauss constraint is conserved in Sect. 5.

A toy model: the 2D wave equation
The basic ideas behind the numerical scheme we are after can be most easily explained using a simple toy model, namely the two-dimensional wave equation. We start by giving a few definitions and then derive three different numerical schemes by discretizing the action of the system in different ways and using a discrete variational principle. The schemes obtained through this procedure are known as variational integrators, which exhibit useful numerical properties such as conserv-ing symplectic structure and retaining symmetries of the discrete action [34,35]. We will see how the exact discretization of the action affects the properties of the numerical scheme and in particular how numerical dispersion can be eliminated.
We consider a real-valued scalar field φ(x) in 2+1 with mostly minuses metric signature (+1, −1, −1) and set the speed of light to c = 1. The action is given by which upon demanding that the variation of the action vanishes yields the equations of motion (EOM) We use Latin indices i, j, k, . . . to denote the spatial components, ∂ 2 i is a shorthand for ∂ i ∂ i (no sum implied) and x = dx 0 dx 1 dx 2 . Inserting a plane-wave ansatz with φ 0 ∈ R and k μ = (ω, k 1 , k 2 ) μ into the EOM gives the dispersion relation Obviously, the phase velocity v = ω/|k| = 1 is constant, i.e. there is no dispersion. Now let us consider a discretized version of this system by approximating space time as an infinite rectangular lattice with grid spacings a μ . We refer to a 0 as the time step and a i as the spatial lattice spacings. The field φ(x) is replaced with field values φ x defined at the lattice sites x and derivatives are replaced with finite difference expressions. We define the forward difference and the backward difference where we introduced another shorthand notation: φ x±μ denotes the field at a neighboring lattice site x ± a μê μ (no implicit sum over μ) with the unit vector in the μ direction e μ . We also define the second order central difference The forward and backward differences are linear approximations to the first order derivative, while the second order difference is accurate up to second order in the lattice spacing a. Equipped with these definitions we could directly discretize the EOM (3), but this is not the approach we will take. The strategy behind variational integrators is to first discretize the action (1) and then demand that the discrete variation vanishes.

Leapfrog scheme
One possible way of discretizing the action is where x is the sum over all lattice sites and V = a 0 a 1 a 2 is the space-time volume of a unit cell. Introducing small variations δφ x of the discrete field at each point, the discrete variation of this action reads which upon setting it to zero yields the discretized EOM where ∂ 2 0 and ∂ 2 i are second order finite differences. Here, we made use of summation by parts, i.e.
which is the discrete analogue of integration by parts. If the field is known in two consecutive time slices we can explicitly solve the EOM (11) for the field values in the next time slice: In fact, this scheme is identical to the explicit leapfrog scheme, 1 which is accurate up to second order in the time 1 The connection to the leapfrog scheme becomes more apparent if we introduce an approximation of the conjugate momentum step a 0 and spatial lattice spacings a i . Using the plane-wave ansatz (4) we find the dispersion relation which is in general non-linear and only yields real-valued (stable) frequencies ω for all wave vectors k if the Courant-Friedrichs-Lewy (CFL) condition holds The discretization errors of this finite difference scheme result in a non-linear dispersion relation, which is usually referred to as numerical dispersion, since this kind of artificial dispersive behavior of plane waves does not show up in the continuum. If it were possible to set a 0 = a 1 = a 2 (the socalled "magic time-step") the leapfrog scheme would actually be non-dispersive along the lattice axes, but this choice of the parameters is forbidden by the CFL condition in higher dimensions than 1 + 1 and would lead to unstable modes.

Implicit scheme
Let us consider a different discretization: we define a new action where φ x is the temporally averaged field Note that only one of the spatial finite differences in the squared term is temporally averaged. Since this action differs from the leapfrog action (9) only up to an error term quadratic in a 0 , the numerical scheme derived from this action will have the same accuracy as the leapfrog scheme.
Repeating the steps as before we obtain the discretized EOM Footnote 1 continued which is defined naturally between time slices x 0 and x 0 + a 0 (hence the index "x + 0 2 " in our notation). The EOM can then be written as This is an implicit scheme, which is more complicated to solve compared to the explicit leapfrog scheme given by Eq. (13). Here we have to find the solution to a system of linear equations, which can be accomplished using (for instance) iterative methods. The dispersion relation for this scheme reads This relation can always be solved for real-valued frequencies ω and therefore the implicit scheme is unconditionally stable. Unfortunately this does not solve the problem of numerical dispersion either, because there is no choice of lattice parameters that results in a linear dispersion relation. We quickly summarize: the first action we considered given by Eq. (9) gave us the explicit leapfrog scheme, which is rendered non-dispersive but unstable using the "magic timestep". The second action, Eq. (19), which we obtained by replacing one of the spatial finite differences with a temporally averaged expression, yields an implicit scheme. This scheme is unconditionally stable, but always dispersive. This suggests that a mixture of both discretizations might solve our problem.

Semi-implicit scheme
Finally, we consider the action where the derivatives w.r.t. x 1 are treated like in the leapfrog scheme and the derivatives w.r.t. x 2 involve a temporally averaged expression as in the implicit scheme. The EOM now read We call this numerical scheme semi-implicit, because the finite difference equation contains both explicitly and implicitly treated spatial derivatives. The dispersion relation associated with Eq. (24) is given by which is stable if The CFL condition (26) now allows us to set a 0 = a 1 . Looking at the dispersion relation (25) we notice that for k 1 = 0, but k 2 = 0 the propagation becomes non-dispersive, i.e. ω = k 1 . For k 2 = 0 and k 1 = 0 the propagation still exhibits numerical dispersion. The scheme defined by the action (23) therefore allows for non-dispersive, stable wave propagation along one particular direction on the lattice. This principle also extends to systems with more spatial dimensions, where one treats a preferred direction explicitly and all other spatial directions implicitly.

Solution method and numerical tests
To solve the EOM of the implicit or the semi-implicit scheme one has to solve a linear system of equations. This can be accomplished for instance by inverting a band matrix. Alternatively the equations can also be solved in an iterative manner. Taking the latter approach will be readily applicable to lattice gauge theory. One example for an iterative method is damped (or relaxed) fixed point iteration: the idea is to first rewrite the EOM (24) as a fixed point equation and then, starting with an initial guess φ (0) x+0 from e.g. the explicit leapfrog evolution, use the iteration to obtain a new guess φ (n+1) x+0 . Here the real-valued parameter α acts as a damping coefficient. Using fixed point iteration might induce other numerical instabilities not covered by the CFL condition (26). To analyze this we make the ansatz where φ (∞) x+0 is the true solution to Eq. (27) and ϕ x represents a time-independent error term. The growth of the error is determined by the modulus of λ. Employing a Fourier ansatz which is independent of the k 1 component and the corresponding lattice spacing a 1 . Requiring convergence for highk modes, i.e. |λ| < 1 for k 2 = ±π/a 2 , we find where Note that for δ < 1/2 damping might not even be necessary. A similar stability condition can be derived also for the implicit scheme, which by itself (just from the plane wave analysis) is unconditionally stable. It is important to keep in mind that the use of fixed point iteration can introduce new instabilities depending on the lattice spacing and the time step. Finally, we perform a crucial numerical test: we compare the propagation of a Gaussian pulse using the three different schemes to show the effects of numerical dispersion and in particular that the semi-implicit scheme is dispersion-free. For simulations using the implicit or semi-implicit method we solve the equations using damped fixed point iteration. The results are shown in Fig. 1.
The main insight of this section is that the specific discretization of the action completely fixes the numerical scheme of the discrete equations of motion (and the associated stability and dispersion properties) which one obtains from a discrete variational principle. The use of temporally averaged quantities in the action leads to implicit schemes. If we treat some derivatives explicitly and some implicitly we can end up with a semi-implicit scheme that can be nondispersive and still stable for propagation along a single direc-tion on the lattice. As it turns out, this is just what we need to suppress the numerical Cherenkov instability we encountered in our heavy-ion collision simulations. In the next section we will see how we can use the same "trick" for Abelian gauge fields on the lattice.

Abelian gauge fields on the lattice
Before tackling the problem of non-Abelian gauge fields on the lattice, it is instructive to see how we can derive a dispersion-free semi-implicit scheme for discretized Abelian gauge fields. We will approach the problem as before: starting with a discretization of the action, we apply a discrete variational principle to derive discrete equations of motion and constraints. Then we will see what modifications to the action are required to obtain implicit and semi-implicit numerical schemes. Since we are dealing with gauge theory we will take care to retain gauge invariance also for the discretized system.
In the continuum the action for Abelian gauge fields reads with the field strength tensor given by The field strength and the action are invariant under gauge transformations defined by where α(x) is an at least twice differentiable function which defines the gauge transformation.
Varying the action with respect to the gauge field leads to The term proportional to δ A 0 (x) leads to the Gauss constraint while the term proportional to the variation of spatial components gives the EOM It is trivial to see that the EOM imply the conservation of the Gauss constraint. As we will see next, it is also possible to formulate a discretization of the system where this conservation of the constraint holds exactly.

Leapfrog scheme
We consider a discretized gauge field A x,μ at the lattice sites x. The field strength tensor F x,μν at x is defined using forward differences which is antisymmetric in the Lorentz index pair μ, ν like its continuum analogue. Furthermore, the lattice field strength is invariant under lattice gauge transformations given by where α x defines the local transformation at each lattice site x. A straightforward discretization of the gauge field action is given by where V = μ a μ is the space-time volume of a unit cell. Due to the invariance of F x,μν under lattice gauge transformations the action is also invariant. Performing the variation of (42) with respect to spatial components A x,i yields the equations of motion. Using the variation of the magnetic part of the action 1 4 and the variation of the electric part w.r.t. spatial components (denoted by δ s ) we find the discrete EOM which are of the explicit leapfrog type. We also obtain a discretized version of the Gauss constraint by considering the variation w.r.t. temporal components A x,0 (denoted by δ t ). With we get the constraint Since both the discrete EOM and the constraint follow from the same discretized, gauge-invariant action (42), the Gauss constraint is guaranteed to be automatically conserved under the EOM. We can show this explicitly via which is equivalent to This means that if the Gauss constraint is satisfied in one time slice then the EOM will ensure that it remains satisfied in the next time slice. We can also give a more general proof: consider an infinitesimal gauge transformation and expand the action S[A ] for small α. We then find Since S[A] is invariant for any α it must hold that or written slightly differently If the equations of motion are satisfied in every time slice, i.e.

∂ S[A]
∂ A x,i = 0, then the Gauss constraint ∂ S[A] ∂ A x,0 must be conserved from one slice to the next: This holds regardless of the exact form of the gauge invariant action S [A]. Consequently, it does not matter what kind of discretization of the action we use. As long as S[A] retains lattice gauge invariance in the sense of Eq. (41), we are guaranteed to find that the discrete Gauss constraint is conserved under the discrete equations of motion. The EOM (45) alone are not enough to uniquely determine the time evolution of the field A x,μ : we must specify a gauge condition. Here we use temporal gauge which we will also use in the case of non-Abelian lattice gauge fields. The EOM (45) then read Using a plane wave ansatz with amplitude A i , we find the same non-trivial dispersion relation and CFL stability condition as in the case of the leapfrog scheme for the wave equation in 2+1 [see Eq. (17) and Eq. (18)]. To show this explicitly we first introduce some notation. Taking either a forward (backward) finite difference of the plane wave ansatz yields which suggests the definition of the forward (backward) lattice momentum We define the squared lattice momentum as Furthermore, we can render the lattice momenta dimensionless by multiplying with a 0 . We define the dimensionless lattice momentum as and The factor of 1/2 in (63) is introduced for convenience. For differences with respect to the time coordinate we find χ 2 0 = sin 2 ωa 0 /2 . Using these definitions the Gauss constraint (47) for the plane wave ansatz can be reduced to The discrete EOM (56) can be written as One can eliminate the mixed terms χ B j χ F i A j using the Gauss constraint and find the dispersion relation which is equivalent to the dispersion relation of the wave equation Eq. (17), i.e.

Implicit scheme
An implicit scheme analogous to the one we derived for the wave equation [see Eq. (19)] can be found using the action where we introduce the temporally averaged field-strength Note that M x,i j differs from F x,i j only by an error term proportional to a 0 2 . Replacing one of the field strengths F x,i j in the quadratic term in the action with its temporally averaged expression M x,i j is analogous to replacing the wave (19). The averaged field strength M x,i j is also invariant under lattice gauge transformations: Varying the action as we did for the leapfrog scheme we find the EOM and employing temporal gauge we have where the fields A x,i have been replaced with temporally averaged expressions A x,i = 1 2 A x+0,i + A x−0,i on the right-hand side. The Gauss constraint that arises from varying w.r.t. temporal components is simply Eq. (47), because we did not modify the term involving F x,0i or introduce new dependencies on A x,0 . The discrete EOM (72) still conserve the Gauss constraint due to lattice gauge invariance. Performing the same steps as in the previous section, we can show that this implicit scheme is unconditionally stable and exhibits the same non-trivial dispersion relation as the implicit scheme for the wave equation, see Eq. (22).

Semi-implicit scheme
In this section we want to develop a semi-implicit scheme for Abelian gauge fields. Specifically we need an action that allows for dispersion-free propagation of waves in the direction of the x 1 coordinate (which we refer to as the longitudinal direction). We call x 2 and x 3 the transverse coordinates. From now on Latin indices i, j, k, . . . refer to transverse indices and the longitudinal index will always be explicit. Our goal is to define the action in such a way that we end up with equations of motion that include explicit differences in the x 1 direction, but temporally averaged finite differences in the x 2 and x 3 direction. This means that we have to modify the F 2 x,i1 term of the leapfrog action (42) such that it results in terms like A first guess for a semi-averaged version of field strength F x,i1 that could accomplish this is with A x,1 = 1 2 A x+0,1 + A x−0,1 . However, it turns out that such a term is not invariant under lattice gauge transformations (41). The problem is that A x,1 transforms differently than A x,i . We have which yields where the last term To fix this we introduce the "properly" averaged field strength A x,1 given bỹ where the last term transforms as Here we made use of the exact relation Therefore, the transformation property of the properly averaged gauge fieldÃ x,1 is the same as A x,1 : It still holds that in the continuum limit the properly averaged gauge fieldÃ x,1 is the same as A x,1 up to an error term quadratic in a 0 . While the definition (78) seems a bit arbitrary at first sight, this way of averaging is more natural using the language of lattice gauge theory: in Sect. 4.3 we will find an intuitive picture in terms of Wilson lines that reduces to (78) in the Abelian limit for small lattice spacings. The properly semi-averaged gauge-invariant field strength is then given by In order to keep the field strength explicitly antisymmetric we define W x,1i = −W x,i1 . Using these definitions we can guess the action where the indices i, j denote transverse components. Since we have built the new action from gauge invariant expressions it is also invariant under lattice gauge transformations. Note that the use ofÃ x,i in W x,1i introduces new terms in the action (83) dependent on the temporal component of the gauge field. Although these terms disappear in temporal gauge (our preferred choice), they still have an effect on the scheme since we have to perform the variation before choosing a gauge. Therefore we will obtain a modified Gauss constraint compatible with the equations of motion derived from the action (83) even after setting A x,0 = 0. At this point one might ask if the "proper" averaging procedure has any effect on the implicit scheme of the previous section as well. It turns out that if one defines the averaged field-strength M x,i j of (70) using the properly averaged gauge fieldÃ x,i , the action of the implicit scheme (and by extension the EOM and the constraint) remains unchanged. This is due to the fact that the terms proportional to A x,0 in (78) cancel: Therefore, no such modification is required in the implicit scheme.
Varying this action w.r.t temporal components A x,0 yields the modified Gauss constraint Since W x,1i explicitly depends on A x,0 , we obtain a correction term to the standard leapfrog Gauss constraint (47). The discrete EOM read We have separate EOM for the longitudinal and transverse components of the gauge field. Note that by replacing the averaged expressions M and W with F the EOM reduce to the leapfrog equations as expected.
The propagation of waves in the semi-implicit scheme turns out to be more complicated compared to the leapfrog or implicit scheme: given a wave vector k and a field amplitude A i (such that the Gauss constraint (85) is satisfied) the dispersion relation becomes polarization dependent, i.e. the scheme exhibits birefringence.
The amplitude A i of an arbitrary plane wave has to be split into a longitudinal and two momentumdependent transverse components where χ F/B i are dimensionless lattice momenta given by Eq. (63). In Appendix A we find that the components A L and A T,2 have the dispersion relation and the component A T,1 has a second different dispersion relation Analyzing the stability of the scheme using the two dispersion relations yields that it is stable if which, when requiring stability for all modes, reduces to If we consider the special case of a plane wave with a purely transverse amplitude and which propagates in the x 1 direction (i.e. setting the transverse momenta χ 2 = χ 3 = 0) we find that both dispersion relations agree: This dispersion becomes linear if we set a 0 = a 1 . This explicitly shows that the semi-implicit scheme for Abelian gauge fields allows for dispersion-free, stable propagation along the longitudinal direction if we use the "magic time-step".
The dispersion relations also agree if we set the longitudinal momentum χ 1 to zero for arbitrary transverse momenta. In general however, wave propagation in this scheme is bifractive. It would be interesting to see if there are alternative discretizations of the action, which allow for dispersion-free propagation without being bifractive.
The main result of this section is the action (83) which gives rise to the semi-implicit scheme. Here we used a combination of differently averaged field strengths, M x,i j and W x,1i , in the action. Our next goal is to generalize these expressions to non-Abelian gauge fields.

Non-Abelian gauge fields on the lattice
The continuum action for non-Abelian Yang-Mills fields is given by where the field strength is The constant g is the Yang-Mills coupling constant and A μ (x) = a A a μ t a is a non-Abelian gauge field, where t a are the generators of the gauge group. In the following we use the normalization tr t a t b = 1 2 δ ab . Through variation of the action we obtain the Gauss constraint and the equations of motion: where the gauge-covariant derivative acting on an algebra element χ is given by In real-time lattice gauge theory, instead of gauge fields A μ (x), we use gauge links (or link variables) U x,μ which are unitary matrices and interpreted as the Wilson lines connecting nearest neighbors on the lattice. U x,μ is the shortest possible Wilson line on the lattice starting at the site x and ending at x + μ. This is also reflected in the gauge transformations where V x is a gauge transformation defined at the lattice site x. Gauge links with negative directions are identified as In the continuum limit the gauge links can be approximated by where A x,μ = a A a x,μ t a . We use "lattice units" for the gauge fields, i.e. we absorb a factor of ga μ in the definition of the gauge field. The gauge links then read For the rest of this paper we drop the hat symbol and just remember to restore factors of ga μ whenever necessary. The smallest possible Wilson loops that can be constructed on the lattice are the so-called "plaquettes" which can also be written as The path traced by the plaquette is shown in Fig. 2 on the left. In the continuum limit the plaquettes can be identified with the field strength tensor Here, F x,μν contains a factor of ga μ a ν . Plaquettes represent a closed Wilson loop and therefore transform locally at the starting (and end) point x: By taking the trace of a plaquette one obtains a gauge invariant expression. In particular it holds that This leads to the standard Wilson gauge action [20] S where i denotes a sum over all spatial components. Using Eq. (108) it is clear that the action (109) is a discretization of the continuum Yang-Mills action where we made the split into temporal and spatial components explicit. Since it is built from gauge-invariant expressions, the Wilson gauge action (109) is invariant under lattice gauge transformations (100). At this point we remark that the continuum Yang-Mills action (110) and its discretization (109) look very different: while the Yang-Mills action is given in terms of squares of the field strength tensor, the Wilson gauge action is linear in plaquette variables. In terms of plaquettes it is not immediately clear how we might generalize our approach from Sect. 3. Fortunately, the action (109) can be written differently so that its functional form is more similar to (110). We define (see for instance p. 94 of [36]) which transforms non-locally For comparison to the plaquette U x,i j , the path traced by C x,i j is shown in Fig. 2 on the right. In the continuum limit C x,μν can be (up to constant factors) identified with the field strength F μν (x): expanding for small lattice spacing we find Most noteworthy is the exact relation with which we can identically rewrite the action as This functional form of the rewritten action (115) is now virtually the same as the continuum case (110). We will exploit this when generalizing the implicit (69) and semi-implicit schemes (83) to non-Abelian gauge fields.
Performing the variation of (109) or (115) is a bit more involved compared to the wave equation or Abelian gauge fields on the lattice. The degrees of freedom, the gauge links U x,μ , are not just completely general complex matrices but det U x,μ = 1.
This means that one has to perform a constrained variation of the degrees of freedom, which can either be done rather tediously by including the constraints explicitly in the action using Lagrangian multipliers or with the help of a variation that conserves the unitary and determinant constraint. An appropriate variation is given by where δ A x,μ = a δ A a x,μ t a is an element of the Lie algebra su(N), i.e. hermitian and traceless. An extended discussion can be found in the Appendix B.

Leapfrog scheme
Following the same procedure as in the case of Abelian gauge fields, we vary w.r.t. spatial components U x,i to obtain the discrete EOM and w.r.t. temporal components U x,0 to find the Gauss constraint.
Starting with the constraint we get (see Appendix C.1 for details) Here we introduced the shorthand P a (. . . ) to denote Varying the spatial link variables we obtain the EOM (see Appendix C.2 for more details) which upon using the definition of C x,i j can be written in the more familiar form As before, time evolution under the EOM conserves the Gauss constraint exactly (see Appendix D). In order to actually solve the equations we specify the temporal gauge which enables us to compute the spatial link variables of the next time-slice using past links and the temporal plaquette: The temporal plaquette U x,0i has to be determined from Eq. (122). For SU(2) this can be done explicitly [17,18]: We use the parametrization where σ a , a ∈ {1, 2, 3} are the Pauli matrices. The realvalued parameters u 0 , u a fulfill the constraint The EOM (122) can then be written as To obtain u 0 we take the positive branch of Eq. (126) assuming that the changes from one time slice to the next are "small", i.e. the temporal plaquette U x,0i will be "closer" to 1 than to −1.

Implicit scheme
Guided by what we learned from the Abelian case in Sect. 3, we would like to replace one of the C x,i j expressions in the Wilson action (115) with a time-averaged equivalent. At the same time we need to retain the gauge invariance of the action. Simply using the temporally averaged expression is not enough because C x+0,i j and C x−0,i j transform differently. A solution is to include temporal gauge links in order to "pull back" C x+0,i j and C x−0,i j to the lattice site x. This leads us to the definition of the "properly" averaged field strength which transforms like C x,i j , i.e.
This gauge-covariant averaging procedure can be generalized: consider an object X x,y that transforms like As an example, X x,y could be a Wilson line connecting points x and y along some arbitrary path. A time-averaged version of X x,y is given by where X x±0,y±0 is simply X x,y shifted up (or down) by one time step. It still transforms like Eq. (132), i.e.
Using this we can write Condensing our notation even further we write and It still holds that so the use of M x,i j instead of C x,i j does not change the accuracy of the scheme. Note that temporal gauge renders all temporal link variables trivial and Eq. (133) reduces to the simple time average.
The inclusion of temporal links in M x,i j breaks a symmetry on the lattice: the diagonal from x to x + i + j is now a preferred direction, which can be seen in Fig. 3. The loss of symmetry can be mitigated by also including terms like M x,i− j in the action. Therefore we propose the action where we explicitly include terms with negative spatial indices using the sum | j| over positive and negative components j to keep the action as symmetric as possible. The action is also invariant under time reversal, real-valued (see Appendix E for a proof) and gauge invariant. While we have not made any changes to the terms involving temporal plaquettes the spatial plaquette terms now include temporal links and therefore we will also obtain a modified Gauss constraint like in the semi-implicit scheme for Abelian fields. Performing the variation the same way we did for the leapfrog scheme, we obtain the Gauss constraint (see where |i| denotes the sum over positive and negative indices i. The left hand side (LHS) of (141) is the same as in the leapfrog scheme Eq. (119), but now there is also a new term on the right hand side (RHS) from varying the spatial part of the action. Performing the continuum limit for the Gauss constraint (after multiplying both sides with a 0 ), temporal link connections only at the start and end points. Therefore the diagonal from x to x + i + j is a preferred direction. This asymmetry can be repaired by including M x,i− j terms in the action (140). We can also see that in the continuous time limit (the vertically stacked time slices would merge into one) the paths traced by M x,i j and W x,1i would become identical to C x,i j the RHS term vanishes as O a 0 2 . This shows that the RHS is not a physical contribution, but rather an artifact of the implicit scheme. Note that the correct continuum limit of the constraint (and the EOM) is already guaranteed by the action.
In a similar fashion as before we perform the variation w.r.t. spatial link variables to get the discrete EOM. We find (see Appendix E.2) which is formally similar to the leapfrog scheme (121) with M x,i j in place of C x,i j , a sum over positive and negative components j (instead of just positive indices) and an additional factor of 1/2 to avoid overcounting. As a simple check one can replace M x,i j with C x,i j in Eq. (142) (only introducing an irrelevant error term quadratic in a 0 ) and recover Eq. (121). Compared to the leapfrog scheme, solving Eq. (142) is more complicated: it is not possible to explicitly solve for the temporal plaquette U x,i0 anymore because M x,i j on the RHS involves contributions from both "past" and "future" link variables. Completely analogous to the case of Abelian gauge fields on a lattice, we obtain an implicit scheme by introducing time-averaged field strength terms in the action.
Similar to what we did in Sect. 2.3, we propose to solve Eq. (142) iteratively using damped fixed point iteration. Start-ing from an initial guess for the future link variable U (0) x+0,i , for instance by performing a single evolution step using the leapfrog scheme, we iterate 1 a 0 a i 2 U a + P a U x,i−0 using U x,i j from the last iteration step to determine U (n+1) x,i0 . In the above equations we replaced P a (U x,i0 ) by an unknown variable U a . We first solve Eq. (143) for U a and then update P a U (n+1) x,i0 using Eq. (144). The parameter α is used as a damping coefficient to mitigate numerical instabilities induced by fixed point iteration. For SU(2) we construct the temporal plaquette U (n+1) using Eqs. (127) and (128) and using temporal gauge we update the link variables via Then we can repeat the iteration and keep iterating until convergence. This iteration scheme can be used to solve the EOM (142) until the Gauss constraint (141) is satisfied up to the desired numerical accuracy. Conversely, this means that unlike the leapfrog scheme, where the Gauss constraint (119) is always satisfied up to machine precision in a single evolution step, the implicit scheme, solved via an iterative scheme, only approximately conserves the Gauss constraint (141). However, in Sect. 5 we will show that using a high number of iterations the constraint can be indeed fulfilled to arbitrary accuracy. In practice however we find that a lower number of iteration is sufficient for stable and acceptably accurate simulations at the cost of small violations of the constraint.
It is also immediately obvious that solving the implicit scheme requires higher computational effort compared to the leapfrog scheme. Considering that one has to use the leapfrog scheme for a single evolution step once (as an initial guess) and then use fixed point iteration, where every step is at least as computationally demanding as single leapfrog step, it becomes clear that the use of an implicit scheme is only viable if increased stability allows one to use coarser lattices while maintaining accurate results.

Semi-implicit scheme
Using our knowledge from Sects. 3 and 4.2 we can now generalize the semi-implicit scheme to real-time lattice gauge theory. An appropriate generalization of the semi-averaged field strength (82) is given by where We also define W x,i1 ≡ −W x,1i . Using the time-averaging notation [see Eq. (133)] this can be written more compactly as where Note that which shows that the semi-averaged field strength W x,1i only differs from C x,1i by an irrelevant error term. Taking the Abelian limit (i.e. neglecting commutator terms) of Eq. (150) and expanding for small lattice spacing yields We find that the linear term of the gauge-covariant average (150) agrees with the expression we constructed in the Abelian semi-implicit scheme (78). While we had to include an "arbitrary" correction term in the Abelian scheme to fix gauge invariance, the link formalism of lattice gauge theory forces us to only consider closed paths constructed from gauge links in the action and thus naturally leads us to the "proper" averaging procedure. The Wilson line path traced by W x,1i , as compared to M x,i j , is shown in Fig. 3 on the right. As before, we keep the action as symmetric as possible by also including terms with negative transverse directions, i.e. W x,1−i . Inspired by the Abelian semi-implicit case (83), we define the new action as where the sum over i and j only run over transverse coordinates and x 1 is the longitudinal coordinate. The purely transverse part of the action uses the same terms as the implicit scheme, see Eq. (140). The longitudinaltransverse part is now given in terms of C x,1 j and W x,1 j analogous to Eq. (69). We have to explicitly include the hermitian conjugate in order to keep the action real-valued.
Varying with respect to temporal components yields the Gauss constraint (see Appendix F.1) where we use the shorthand Varying w.r.t. spatial links the discrete semi-implicit EOM read (see Appendix F.2) 1 a 0 a 1 2 P a U x,10 + U x,1−0 and for transverse components where |1| simply means summing over the terms with positive and negative longitudinal directions. We now have two sets of equations: one for longitudinal components, and two for transverse components. The equations of motion can be written more compactly by introducing the symbol where C can be exchanged for corresponding expressions with M or W and U can be exchanged for its temporally averaged version U . The longitudinal component of the EOM then reads 1 a 0 a 1 2 P a U x,10 + U x,1−0 and the transverse components are given by These equations can be solved numerically using damped fixed point iteration completely analogously to Sect. 4.2. First, one obtains an initial guess U (0) x+0,i for "future" link variables from a single leapfrog evolution step using Eqs. (122) and (124). Then iterate from n = 1 until convergence: 1. Compute the next iteration using damped fixed point iteration: in Eqs. (156) and (157) replace P a U x,10 → U a 1 and P a U x,i0 → U a i , solve for the unknown U 's and update the temporal plaquettes using P a U (n) x,10 = α P a U (n−1) x,10 and analogously for U (n) x,i0 and U a i . α is the damping coefficient. 2. For SU (2) we can reconstruct the full temporal plaquette from its components P a (U ) with the identity x, 10 and U = U (n) x,i0 . 3. Using U (n) x, 10 and U (n) x,i0 , compute the spatial links U 4. Repeat with n → n + 1.
As with the implicit scheme, our approach to solving the equations in the semi-implicit scheme is an iterative one. The Gauss constraint (154) is only approximately satisfied, depending on the degree of convergence. 2

Coupling to external color currents
Up until now we have only considered pure Yang-Mills fields.
In the continuum we can include external color currents by adding a J · A term to the action.
The equations of motion then read and due to gauge-covariant conservation of charge we have which is the non-Abelian continuity equation. On the lattice we can simply add a discrete J · A term to the action as well: where A a x,μ includes a factor of ga μ ("lattice units"). We also made the split into 3+1 dimensions explicit using J a 0 (x) ρ a x and J a i (x) j a x,i . The variation of S J simply reads which gives the appropriate contributions to the Gauss constraint and the EOM. In the leapfrog scheme the constraint now reads and the EOM read The constraint taken together with the EOM imply the local conservation of charge (see Appendix D) which is the discrete version of the continuity equation. For the implicit and semi-implicit schemes the procedure is the same: including the S J term simply leads to the appearance of ρ on the RHS of the Gauss constraint (see Eqs. (141) and (154)) and j x,i on the RHS of the EOM (see Eqs. (142) and (156), (157)). Due to conservation of the Gauss constraint without external charges, the continuity equation for the implicit and semi-implicit scheme is simply Eq. (171) as well. This implies that our treatment of the external currents in terms of parallel transport as detailed in our previous publication [14] does not require any modifications when using the newly derived schemes.

Numerical tests
In this last section we test the semi-implicit scheme on the propagation of a single nucleus in the CGC framework. For an observer at rest in the laboratory frame, the nucleus moves at the speed of light and consequently exhibits large time dilation. As the nucleus propagates, the interactions inside appear to be frozen and the field configuration is essentially static. On the lattice we would like to reproduce this behavior as well, but depending on the lattice resolution we run into the numerical Cherenkov instability, which leads to an artificial increase of the total field energy of the system. As previously stated the root cause of the instability is numerical dispersion: in the CGC framework a nucleus consists of both propagating field modes and a longitudinal current generating the field around it. It is essentially a non-Abelian generalization of the field of a highly relativistic electric charge. In our simulation the color current is modeled as an ensemble of colored point-like particles moving at the speed of light along the beam axis. The current is unaffected by dispersion, i.e. it retains its shape perfectly as it propagates. The field modes suffer from numerical dispersion, which over time leads to a deformation of the original longitudinal profile. The mismatch between the color current and the field leads to creation of spurious field modes, which interact with the color current non-linearly through parallel transport (color rotation) of the current. This increases the mismatch further and more spurious fields are created. As the simulation progresses this eventually leads to a large artificial increase of total field energy. The effects of the instability can be quite dramatic as seen in Fig. 4. The main difference to the numerical Cherenkov instability in Abelian PIC simulations is that in electromagnetic simulations the spurious field modes interact with the particles through the Lorentz force [27]. In our simulations we do not consider any acceleration of the particles (i.e. their trajectories are fixed), but interaction is still possible due to non-Abelian charge conservation (171), which requires rotating the color charge of the color current. Therefore our type of numerical Cherenkov instability is due to non-Abelian effects.
We now demonstrate that the instability is cured (or at least highly suppressed for all practical purposes) by using the semi-implicit scheme. We test the schemes ability to improve energy conservation the following way: we place a single  Fig. 4 The energy density of a right-moving nucleus averaged over the transverse plane as a function of the longitudinal coordinate x 1 after t = 2 fm/c. The results were obtained from the same simulations as in Fig. 5. We compare the performance of the leapfrog (LF) scheme (with N L = 256 and N L = 1024) to the semi-implicit (SI) scheme with N L = 256. In the most extreme example (LF 256) the nucleus becomes completely unstable due to the numerical Cherenkov instability. By eliminating numerical dispersion (SI 256) the nucleus retains its original shape almost exactly gold nucleus described by the McLerran-Venugopalan model in a simulation box of volume V = (3 fm) × (6 fm) 2 and set the longitudinal extent of the nucleus to roughly correspond to a boosted nucleus with Lorentz factor γ = 100. These parameters correspond to a similar setting as in our previous work [23]. We refer to [22] for a detailed description of the initial conditions. After setting up the initial condition we let the nucleus freely propagate along the longitudinal axis. As the simulation runs we record the total field energy where E a i (x) and B a i (x) are the color-electric and -magnetic fields at each time step. On the lattice the electric and magnetic fields are approximated using plaquettes: We compute the relative change (E(t) − E(0)) /E(0), which we plot as a function of time t. In the continuum we would have E(t) = E(0), but due to numerical artifacts and the Cherenkov instability this is not the case in our simulations.
The numerical results are shown in Fig. 5. We see that the leapfrog scheme leads to an exponential increase of the total energy over time, which can be suppressed using finer lattices. On the other hand, the semi-implicit scheme leads to better energy conservation even on a rather coarse lattice. Therefore, the resolution that is usually required to obtain accurate, stable results is lowered by using the semi-implicit   Fig. 5 The relative increase of the total field energy E(t) as a function of time t for the propagation of a single nucleus in a box of volume V = (3 fm) × (6 fm) 2 with a longitudinal length of 3 fm and a transverse area of (6 fm) 2 . The lattice size is N L × N 2 T with the transverse lattice fixed at N T = 256. Starting with the same initial condition, we evolve forward in time using the leapfrog (LF) and the semi-implicit (SI) method. In the case of the leapfrog scheme we vary the resolution along the beam axis using the number of longitudinal cells N L of the lattice. For N L = 256 the numerical Cherenkov instability leads to catastrophic failure, increasing the energy to many times its original value. The effect is suppressed when increasing the longitudinal resolution, but the instability is still present. In the case of the semi-implicit scheme it is possible to set N L = 256 and still obtain (approximate) energy conservation. After t = 2 fm/c the energy increase for the semiimplicit scheme is roughly 0.02%, compared to 1% for the leapfrog with N L = 2048. For the simulation using the semi-implicit scheme we used N i = 10 iterations and a damping coefficient of α = 0.45. The time step is set to the longitudinal lattice spacing. In the case of the leapfrog simulation we used a 0 = a 1 /4 scheme. However, using the new scheme might not always be economical: finer lattices suppress the instability as well and since the leapfrog scheme is computationally cheaper than the semi-implicit scheme, the leapfrog can be favorable in practice. On our test system (a single 256 GB node on the VSC 3 cluster) the simulation using the semi-implicit scheme (SI 256) takes ∼ 4 hours to finish, while the same simulation using the leapfrog with N L = 1024 (LF 1024) takes roughly ∼ 2.5 hours and with N L = 2048 (LF 2048) ∼ 10 hours. Even though energy conservation is not as good as SI 256, the longitudinal resolution is much better in comparison enabling us to extract observables with higher accuracy. It should be noted however that our implementation of the leapfrog scheme is already highly optimized, while the implementation of the semi-implicit scheme is very basic and should be considered as a proof of concept. Further optimizations and simplifications of the semi-implicit scheme might make it the better choice in most cases.
As a second test we look at the violation of the Gauss constraint. The leapfrog scheme (122) conserves its associated Gauss constraint (119) identically, even for finite time-steps a 0 . In numerical simulations this conservation is not exact due to floating point number errors, but the violation is zero up to machine precision. On the other hand, the semi-implicit We also compare to the constraint violation after only a single evolution step (red crosses), which does not differ much from the violation after a larger number of time steps. It is evident that the violation systematically converges towards zero (up to machine precision) as we increase the number of iterations N i scheme has to be solved iteratively and therefore the results depend on the number of iterations N i used in the fixed point iteration method. We define the relative violation of the Gauss constraint as the ratio of the absolute (squared) Gauss constraint violation to the total (squared) charge on the lattice. In the case of the leapfrog scheme this reads where the sum x runs over the spatial lattice of a single time slice at t. The numerator depends on the Gauss constraint of the scheme and has to be adjusted according to the implicit and semi-implicit method (either Eq. (141) or (154) including the charge density on the RHS as discussed in Sect. 4.4). In Fig. 6 we show how the Gauss constraint violation converges systematically towards zero as we increase the number of iterations. Therefore, even though we can not use an arbitrarily high number of iterations due to limited computational resources, the semi-implicit scheme conserves the Gauss constraint in principle. The same holds for the purely implicit scheme. In practice it is not necessary to satisfy the constraint up to high precision, as observables such as the energy density seem to converge much faster up to satisfying accuracy.

Conclusions and outlook
In this paper we derived new numerical schemes for realtime lattice gauge theory. We started our discussion based on two simpler models, namely the two-dimensional wave equation and Abelian gauge fields on the lattice. It turns out that using a discrete variational principle to derive numerical schemes for equations of motion is a very powerful tool: the use of time-averaged expressions in the discrete action yields implicit and semi-implicit schemes depending on how exactly the time-averaging is performed and what terms are replaced by their averages. We extended this concept to real-time lattice gauge theory, allowing us to make modifications to the standard Wilson gauge action that yield new numerical schemes, which have the same accuracy as the leapfrog scheme and, most importantly, are gauge-covariant and conserve the Gauss constraint. Finally, we demonstrated a peculiar property of the semi-implicit scheme: it allows for dispersion-free propagation along one direction on the lattice, thus curing a numerical instability that has plagued our simulations of three-dimensional heavy-ion collisions. Although all numerical tests in this work have been performed for the propagation of a single nucleus, we expect that the new semi-implicit scheme will improve simulations of nucleus-nucleus collisions in multiple ways. Primarily, using the new scheme we can be sure that the color fields of incoming nuclei have not been altered up until the collision event and all changes to the fields afterwards are solely due to interaction between the colliding nuclei during the collision event itself. This also helps to improve numerical accuracy in the forward and backward rapidity region: at later simulation times, when the now outgoing nuclei are well separated, all field modes (i.e. "gluons") with momenta at high rapidity must have been created in the collision and contributions from artificial modes emitted by the nucleus due to numerical Cherenkov radiation are strongly suppressed. Furthermore, from the dispersion relations (90) and (91) we can infer that these gluons with almost purely longitudinal momentum k L (and small transverse momentum k T ) exhibit a phase velocity approximately the speed of light (up to an error term quadratic in a 0 k T ). This means any interactions between high rapidity gluons produced in the collision and the finite-thickness color fields of the nuclei directly after the collision event can be considered physical and are not tainted by numerical dispersion as in our previous simulations. Consequently, it should be possible to extract space-time rapidity profiles of the local rest frame energy density of the Glasma as in [23] valid for larger ranges of rapidity as previously considered.
In conclusion, we hope that this new treatment of solving the Yang-Mills equations on the lattice allows us to perform better, more accurate simulations using more complex models of nuclei. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP 3 .

Appendix A: stability analysis of the semi-implicit scheme for Abelian gauge fields
In this section of the appendix we prove stability for the semiimplicit scheme for Abelian fields derived in Sect. 3.3. Using temporal gauge, A x,0 = 0, the equations of motion (86) and (87) read where W x,1i reduces to The Gauss constraint (85) reads Splitting the equations of motion into Laplacian terms and mixed derivative terms we find Using a plane wave ansatz we will use the Gauss constraint to first reduce the number of degrees of freedom and then compute the dispersion relation ω(k). Inserting the ansatz into the Gauss constraint we find where we use the dimensionless lattice momenta (63). The constraint equation can alternatively be written as where β is a momentum-dependent factor given by The temporal average A x,i reduces to a multiplication with a frequency dependent factor: where we used the shorthand c = cos ωa 0 . Inserting the plane wave ansatz into the EOM yields Note that both χ 0 and c depend on ω. After making use of the Gauss constraint the longitudinal EOM reads .14) and the two transverse equations read This system of equations can be written in matrix notation as an eigenvalue problem where the vector A is simply The eigenvectors of M are where we find two transverse, momentum dependent eigen- and ω T,2 a 0 = arccos 1 − χ 2 1 2 + χ 2 2 + χ 2 It turns out that the expressions for ω associated with the different eigenvectors are not the same, although A L and A T,2 share the same dispersion relation. We interpret this as numerical (or artificial) birefringence. Given a momentum k, the amplitude of an arbitrary wave has to be split into two components: a part which is projected into the plane spanned by A L and A T,2 which oscillates with ω L = ω T,2 and a part parallel to A T,1 which oscillates with frequency ω T,1 . We require the propagation of a wave to be stable, i.e. we require the frequencies ω to be real-valued. This is guaranteed if the arguments of the arccos expressions in Eqs. This concludes the proof that the semi-implicit scheme, even though exhibiting peculiar wave propagation phenomena, is stable.

Appendix B: variation of gauge links
We introduce the infinitesimal variation of a gauge link variable where the variation of the gauge field δ A x,μ is traceless and hermitian. In the continuum limit δ A x,μ becomes the infinitesimal variation of the gauge field A μ (x). The infinitesimal variation δU x,μ preserves the unitary of gauge links. Let U x,μ = U x,μ + δU x,μ then we find The determinant is also unaffected for infinitesimal variations. det The variation δU x,μ therefore preserves the constraints and allows us to vary the action without Lagrange multipliers, which dramatically simplifies the derivation of equations of motion.

Appendix C: variation of the leapfrog action
In this section of the appendix we give a derivation of the discrete equations of motion obtained from the standard Wilson action (109) using the constraint preserving variation of link variables of the previous section. To make the calculation more organized, we first split the action into two parts: a part containing temporal plaquettes S E [U ] ("E" for electric) and a part containing spatial plaquettes S B [U ] ("B" for magnetic). We write To go from the first to the second line, we applied a shift x → x − i in the right term of the first line and made use of the cyclicity of the trace. In the third line we simply used the definition of the variation of gauge links. The variation of the action therefore reads where we used Replacing the C x,i j terms with link variables we find The variation then reads (C.43) We set δS = 0 and after canceling some constants we find the discrete EOM 1 a 0 a i 2 P a U x,i0 + U x,i−0 which can also be written as (C.45)

Appendix D: conservation of the Gauss constraint in the leapfrog scheme
We now explicitly show that the leapfrog EOM preserve the associated Gauss constraint. We use the identity for the fundamental representation of SU(N) which can be shown using the Fierz identity for the generators of su(N) a t a i j t a kl = where i, j, k, l are fundamental representation matrix indices. Using a shorthand we can write a t a P a (X ) = [X ] ah , (D.48) where "ah" denotes the anti-hermitian traceless part of X . The constraint and the equations of motion then read (including external charges) We take (D.50) and sum over i. Due to U x,i j ah being antisymmetric in the index pair i, j we find Doing the same at x − i and parallel transporting from (D.52) Subtracting the above two equations gives Using antisymmetry we have and where we used the Gauss constraint in the last line to replace the temporal plaquette terms with charge densities. This yields the gauge-covariant continuity equation If there are no external charges, we simply have the conservation of the Gauss constraint: assume that the constraint in the previous time slice holds, i.e. then the EOM guarantee that it will also hold in the next one, i.e.

Appendix E: variation of the implicit action
We now consider the action The electric part is the same as in the leapfrog scheme. However, the magnetic part S B [U ] now contains temporal gauge links and gives a new contribution to the Gauss constraint. Before we vary this action, we must verify that S B [U ] is indeed real-valued. While it is easy to see that the original leapfrog action is real-valued because of the obvious hermicity of C x,i j C † x,i j , the term C x,i j M † x,i j is not hermitian in general. Still, we can show that the action is real: we have where we can rewrite Here we used the cyclicity of the trace. Then using a shift x → x − 0 we have This means that under the sum over x and the trace, the expression C x,i j M † x,i j is indeed real-valued and by extension S B [U ] is real-valued as well. We have also shown that the time-average in M x,i j can be "shifted" to the other term C x,i j under the sum and trace, i.e.
This is a useful property that we will need in the following derivation. On the other hand, the variation of the magnetic part involves terms like δ t C x,i j M † x,i j , which we now discuss explicitly. First, we make use of Then, after some algebra we find and δ t C x,i j C † x,i j − h.c.
The variation of the magnetic part therefore reads x,a,i,| j| x,a,|i|,| j| 1 8 1 a i a j 2 δ A a x,0 P a C (+0) x,i j C † x,i j . (E.74) In the last line we consolidated terms with index i and −i into a single term using the sum |i| . Taking the result for δ t S E [U ] from the leapfrog scheme, we find the Gauss constraint