RRothkopf I NSTITUTE FOR M ATHEMATICS AND P HYSICS T HEORETICAL P HYSICS D IVISION

RESEARCH ARTICLE - PREPRINT

Conserving Lattice Gauge Theory for FiniteSystems

Alexander Rothkopf

Correspondence:[email protected] of Science and Technology,University of Stavanger, 4021Stavanger, NorwayFull list of author information isavailable at the end of the article

Abstract

In this study I develop a novel action for lattice gauge theory for ﬁnite systems,which accommodates non-periodic boundary conditions, implements the properintegral form of Gauss’ law and exhibits an inherently symmetric energymomentum tensor, all while realizing automatic O ( a ) improvement. Taking themodern summation-by-parts formulation for ﬁnite diﬀerences as starting pointand combining it with insight from the ﬁnite volume strategies of computationalelectrodynamics I show how the concept of a conserving discretization can berealized for non-Abelian lattice gauge theory. Major steps in the derivation areillustrated using Abelian gauge theory as example. Keywords: lattice gauge theory; boundary conditions; summation by parts

Revealing the dynamical properties of matter under extreme conditions, be it inthe context of ultra-high temperature and density, or ultra-strong interparticle cou-pling, lies at the heart of modern theoretical physics. Non-perturbative access tothe quantum many-body physics of such systems is called for and lattice gaugetheory has established itself as a vital tool [1]. It provides the conceptual andcomputational basis to predict phenomenologically relevant observables over a vastrange of energies, one pertinent example being lattice Quantum-Chromodynamics(QCD). Lattice QCD predicts from ﬁrst principles not only the properties of sta-ble hadrons (see e.g. [2, 3]) and resonances [4] but captures the thermodynamicproperties of a strongly interacting many-body system, such as the quark-gluonplasma (QGP) [5, 6]. On the other hand at lower energies, recent experimental ad-vances in the manipulation of conﬁned light ﬁelds have made accessible a regime ofnon-perturbatively large couplings in cavity quantum electrodynamics (QED) (forrecent reviews see [7, 8]), which deﬁes the standard methods of perturbation theory.A lattice ﬁeld theoretical approach to this deep-strong coupling (DSC) regime, inparticular, in the presence of large number of interacting dipoles promises urgentlyneeded insight (see also [9]) into an as of yet unexplored realm.In each case, the study of ﬁnite quantum systems and the physical consequencesof non-trivial boundary conditions play a central role. At high temperatures, theunexplained signatures of collectivity in small collision systems (see [10] for a review) a r X i v : . [ h e p - l a t ] F e b othkopf Page 2 of 9 urge theory to provide non-perturbative insight into strongly-interacting many-bodysystems of ﬁnite extent (see e.g. [11, 12]). At low energies the conﬁnes of the reﬂectivecavities used to manipulate the photon-atom interactions leaves an essential imprinton spectra both on the single particle level [13] as well as on collective excitations.In order to make decisive progress in these ﬁelds, non-perturbative theory supportis essential, which however can only be provided by a proper lattice gauge theoryfor ﬁnite systems.According to Wilsons seminal work [14], the simplest gauge invariant action fornon-Abelian ﬁelds on a lattice with periodic boundary conditions (PBC), spatialand temporal lattice spacings a s , a t and grid points N s , N t can be expressed [15] as S = g (cid:80) x (cid:2) (cid:80) i a s a t ReTr[ P i − − (cid:80) i,j a t a s ReTr[ P ij − (cid:3) , in terms of link variables U µ,x = exp[ ia µ A µ,x ], combined in elementary plaquettes P × µν,x = U µ,x U ν,x + a µ ˆ µ U † µ,x + a ν ˆ ν U † ν,x = e ia µ a ν ˜ F µν,x + O ( a ) , ˜ F µν = ∆ F µ A ν,x − ∆ F ν A µ,x + i [ A µ,x , A ν,x ] . (1)The forward ﬁnite diﬀerences (FD) are deﬁned as ∆ F µ φ ( x ) = ( φ ( x + a µ ˆ µ ) − φ ( x )) /a µ and the generators T a of the gauge group enter as A µ,x = A aµ,x T a .The past decades have seen signiﬁcant progress in the formulation of lattice gaugetheory beyond eq. (1). Symanzik’s improvement program [16] laid out how to sys-tematically reduce discretization errors by adding higher dimensional operators,leading to highly improved actions, such as in [17–22]. On the other hand optimizedalgorithms exploit the structure of such actions to realize eﬃcient Monte-Carlo sim-ulations of the corresponding Euclidean path integral [23]. Two challenges exist in deploying the discretization of eq. (1) in a ﬁnite system,where PBC are inapplicable and translational invariance is broken: (A) a naive ﬁrstorder FD scheme does not accurately reproduce the true solution of a diﬀerentialequation if it remains restricted to the interior of the solution domain. Let us take alook at the Abelian Gauss law ∇ · E = 0 in the interior of a ﬁnite circular capacitor,whose end-plates are located on the boundary in the z-direction, held at a ﬁnitepotential. The accurate solution based on the electric potential in the centered yzplane is shown as red arrows in the left half of the left panel of ﬁg. 1. Setting thissolution as Dirichlet boundary for a naive backward FD approximation ∆ B · E =0 we obtain the blue arrows instead. The approximation is unable to access theforward facing boundary and thus cannot sustain ﬁeld strength close to it. The resultimproves when using a naive central FD ∆ C µ φ ( x ) = ( φ ( x + a µ ˆ µ ) − φ ( x − a µ ˆ µ )) / a µ in ∆ C · E = 0. As seen from the green arrows in the right panel of ﬁg. 1, ﬁeldstrength is generated close to the forward boundary but the ﬁeld lines show anartiﬁcial staggered pattern.The correct treatment of non-trivial boundary conditions in FD approximationsis reﬂected in their implementation of the discrete counterpart to integration byparts, the so called summation by parts (SBP) property (for recent reviews see[24, 25]). Using the trapezoid rule for integration (cid:82) L dxf ( x ) g ( x ) ≈ T N [ f x g x ] , one othkopf Page 3 of 9 Figure 1 (left) Electric ﬁeld (red arrows) from the capacitor potential (background) vs. backwardFD solution (blue arrows). (right) The ﬁeld strength from a naive central FD (green arrows) andthe more accurate solution based on a SBP discretization (red arrows). must require the ﬁrst order FD approximation to fulﬁll T N [(∆ SBP f x ) g x ] ! = −T N [ f x (∆ SBP g x )] + f N g N − f g . (2)The lowest order implementation of such an SBP operator is obtained by combiningthe central FD stencil in the interior with the forward FD and backward FD stencilat the boundaries. The resulting operator is anti-symmetric (corresponding to ahermitean momentum operator on the whole domain) and when used in ∆ SBP · E = 0yields the signiﬁcantly more accurate solution for the capacitor denoted by the redarrows in the right panel of ﬁg. 1.(B) The second challenge pertains to the fact the naive forward FD approximationof eq. (1) introduces signiﬁcant artifacts (in particular anisotropy) in the energymomentum tensor (EMT). This issue manifests itself already in simple cases, such asin the presence of static charges. The concept of an electric ﬁeld becomes ambiguousin non-Abelian gauge theory, where one must instead resort to computing the ﬁeldstress tensor Θ ij , i.e the spatial components of the gauge independent EMT T αβ = π (cid:0) g αµ F µλ F λβ + g αβ F µλ F µλ (cid:1) [26]. The stress tensor yields the force density f = ∇ · Θ + ∂ S /∂t acting on a charge after diﬀerentiation and S denotes the Poyntingvector, which is irrelevant in the static case [27]. The force ﬁeld lines are computedvia diagonalizing the stress tensor in Minkowski time. It features two negative andone positive eigenvalue, all of the same magnitude and its eigenvectors (EV) indicatethe direction of the forces acting on its faces. Placing a charge at ( x ) and anti-charge at ( x ), centered in x and y and separated by a distance in z, the true stresstensor in the centered yz plane, due to symmetry, exhibits two eigenvectors in thatplane and one perpendicular to it. The EV associated with the positive eigenvalueindicates the direction of the force lines, as shown in the analytic solution (redarrows) in the left panel of ﬁg. 2If instead we compute the stress tensor from the solution of the discrete Gausslaw ∆ B · E = a s (cid:2) δ xx − δ xx (cid:3) we ﬁnd that in the yz plane only a few of its EV othkopf Page 4 of 9 z Figure 2 (left) Force ﬁeld lines of the stress tensor between an Abelian unit charge andanti-charge. (center) The EVs obtained from a backward FD solution of Gauss’ law. Gray arrowsdenote the only vectors truly in the yz plane, blue arrows the yz component of those predominantlyin the zy plane. Analytic solution is shown in light red (right) solution from the central FD withpoint sources (green arrows) and the conserving solution with distributed sources (red arrows). (gray arrows in the center panel of ﬁg. 2) are truly parallel to it. If one plots thevectors with more than 50% of their magnitude in that plane, one obtains the bluearrows. Compared to the true solution given as light red arrows, a clear asymmetryis visible in the forward FD solution, i.e. the force ﬁeld is not accurately reproduced.Solving Gauss’ law with the central FD discretization we ﬁnd the symmetric greenarrows in the right panel, which however exhibit an artiﬁcial pattern of force linesoriented 90 degreed to each other.The reason for the failure of the naive central FD approximation lies in its inabil-ity to reproduce the integral form of Gauss’ law, a fact well known in computationalelectrodynamics [28, 29]. Indeed the third equality in Q = (cid:82) dV q = (cid:82) dV ( ∇ E ) = (cid:82) ∂V d A · E does not hold in general for the central FD. As is standard in moderncomputational electrodynamics, a genuinely conserving ﬁnite volume discretization[30], i.e. one that obeys both the diﬀerential and integral Gauss law can be con-structed by starting from discretizing the Gauss law with the midpoint rule (cid:90) x i +1 / x i − / dx (cid:90) y i +1 / y i − / dy (cid:90) z i +1 / z i − / dz (cid:0) dE x dx + dE y dy + dE z dz (cid:1) = (cid:90) d xδ (3) ( x − x ) . (3)The resulting equations on the LHS can be brought into the form of a central ﬁnitediﬀerence by adding multiple shifted versions of eq. (3). This leads us to (cid:88) i ∆ C i E i ( x ) = 18 a (cid:104) (cid:88) i (cid:0) δ x + a ˆ i , x + δ x − a ˆ i , x (cid:1) − δ x , x (cid:105) . (4) othkopf Page 5 of 9 We see that in order for the central FD approximation to implement the integralGauss law we must employ a spatially distributed source (c.f. smearing techniques inconventional lattice QCD). Equation (4) amounts to a ﬁnite volume discretization ofGauss’ law. Solving the above yields the symmetric and signiﬁcantly more accuratesolution shown as red arrows on the right of ﬁg. 2.We conclude from the Abelian Gauss law that a proper lattice gauge theoryfor ﬁnite systems must both implement the summation by parts property and beformulated with distributed sources. In the remainder of this letter I construct alattice action that implements these prerequisites.

In the preceding section we used Gauss’ law as guide towards an appropriate de-scription for ﬁnite systems and discovered that it requires implementing the sum-mation by parts property (A) and the presence of distributed sources (B). To con-struct such an SBP discretization for non-Abelian gauge theory we need to ﬁnd thecentral FD counterpart to eq. (1), operating at global order O ( a ) in the latticespacing. We thus need links of the corresponding order ¯ U µ,x = exp (cid:2) igA µ,x + a ˆ µ (cid:3) =exp (cid:2) ig (cid:0) A µ,x + A µ,x + a ˆ µ (cid:1)(cid:3) + O ( a ). The central FD requires access to gauge ﬁelds inthe forward and backward direction in each dimension. The product of eight links(see ﬁg. 3) centered around the current spacetime point x includes exactly thesecontributions P × µν,x = ¯ U µ,x − a ˆ µ − a ˆ ν ¯ U µ,x − a ˆ ν ¯ U ν,x + a ˆ µ − a ˆ ν ¯ U ν,x + a ˆ µ × ¯ U † µ,x + a ˆ ν ¯ U † µ,x − a ˆ µ + a ˆ ν ¯ U † ν,x − a ˆ µ ¯ U † ν,x − a ˆ µ − a ˆ ν =exp (cid:2) iga µ a ν ¯ F µν,x (cid:3) + O ( a )¯ F µν,x = ∆ C µ A ν,x − ∆ C ν A µ,x + i [ A µ,x , A ν,x ] (5)A gauge invariant action with the correct classical continuum limit may now beconstructed from P × µν,x in the interior of the ﬁnite volume V S × = (cid:88) x/ ∈ ∂V a t a s (cid:104) a t a s (cid:88) i ReTr (cid:2) − P × i,x (cid:3) − a s (cid:88) ij ReTr (cid:2) − P × ij,x (cid:3)(cid:105) (6)Lattice discretization errors start at one higher power of the lattice spacing as inthe Wilson action, making S × automatically O ( a ) improved. Note the diﬀerenceto the Symanzik improved actions of Refs. [17–19, 22], which contain an admixtureof asymmetric Wilson plaquettes in the interior and thus do not correspond to aSBP compatible central ﬁnite diﬀerence.Next we combine the central FD of the interior with the forward and backwardFD on the boundaries, which is achieved by using the forward and backward Wilsonplaquettes there in addition to S × , leading ﬁnally to˜ F µν,x = ∆ SBP µ A ν,x − ∆ SBP ν A µ,x + i [ A µ,x , A ν,x ] . (7) othkopf Page 6 of 9 This concludes the solution to challenge (A) for non-Abelian gauge theory. Let menote that simply replacing the Wilson plaquette by the cloverleaf sum does not dothe trick: computing ∂S CL /∂A a leads to the same backward FD scheme for theGauss law that one obtains from the naive plaquette action itself. P × as productof four shifted P × ’s adds the exponents, not the exponentials. Wilson plaquette clover leaf this study x x x μ μ μ νa μ a ν Figure 3

Plaquette actions: (left) Wilson (center) clover-leaf (right) stand-alone centered × plaquette. Let’s now turn to challenge (B). In contrast to computational electrodynamics,which is formulated on the level of the equations of motion, we construct an ac-tion based formulation, applicable to the quantum realm. Thus one cannot simplyintroduce spatially distributed dynamical sources by hand but needs to prescribehow such a term arises from coupling to matter ﬁelds. Let me demonstrate thestrategy for static matter ﬁelds, in which only the covariant derivative in time di-rection plays a role [31]. For static fermions, we may decouple their four-spinorinto two Pauli spinors Φ = ( ψ, χ ). Let us focus on the action for ψ , as that of χ follows from complex conjugation: S ψ = a t a s (cid:80) x iψ † x D ψ x . In order to obtainthe shift pattern of eq. (4) required for a conserving discretization, we will haveto modify the covariant derivative D from its standard form. Taking inspirationfrom a recently developed FD operator for the simulation of open quantum systemdynamics [32] we may amend the standard temporal derivative on an isotropiclattice by an expression that extends it into one of the spatial directions j as∆ RN − SBP t ( j ) φ ( x ) = ( φ ( x + a ˆ0 − a ˆ j ) − φ ( x − a ˆ0 − a ˆ j ) − φ ( x − a ˆ0+ a ˆ j )+ φ ( x + a ˆ0+ a ˆ j )) / a .One symmetric and gauge covariant formulation of the fermion action thus reads S ψ = a t a s i (cid:88) x,i ψ † x ¯ D i ) ψ x = a t a s i (cid:88) x,i ψ † x a (cid:104) (8) (cid:0) U ,x U † j,x + a ˆ0 − a ˆ j + U † j,x − a ˆ j U ,x − a ˆ j (cid:1) ψ x + a ˆ0 − a ˆ j − (cid:0) U † j,x − a ˆ j U † ,x − a ˆ0 − a ˆ j + U † ,x − a ˆ0 U † jx − a ˆ0 − a ˆ j (cid:1) ψ x − a ˆ0 − a ˆ j − (cid:0) U † ,x − a ˆ0 U j,x − a ˆ0 + U j,x U † ,x − a ˆ0+ a ˆ j (cid:1) ψ x − a ˆ0+ a ˆ j +( U j,x U ,x + a ˆ j + U ,x U j,x + a ˆ0 (cid:1) ψ x + a ˆ0+ a ˆ j (cid:105) In the construction of the classical equation of motion in the presence of fermionsas outlined in Ref. [33], the charge density for ψ is obtained from the expectationvalue of the fermion two-point function q x = g Tr (cid:2) (cid:104) [ ψ † x , ψ x ] (cid:105) . Expressing the time othkopf Page 7 of 9 derivative as shown above will introduce symmetric forward and backward shiftsin this term, which can be combined into the conserving discretization of eq. (4).For truly static charges on the other hand, the spatial distribution of source termsin Equation (4) can be implemented directly through a modiﬁcation of the sourceterms within the Wilson loop. I have argued and presented numerical evidence in the Abelian theory that in orderto capture the physics of non-trivial boundary conditions, the naive forward ﬁnitediﬀerence discretization underlying the Wilson action must be replaced with onethat respects the summation by parts property. In addition, even in the presenceof trivial periodic boundary conditions, the asymmetric Wilson action, as well asthe clover leaf, both lead to a backward ﬁnite diﬀerence Gauss law, which as Ihave shown imprints asymmetries into the force ﬁeld lines encoded in the energymomentum tensor. The standard Symanzik improvement, while reducing the magni-tude of the discretization error, does not alleviate the asymmetry of the underlyingdiscretization. The naive central ﬁnite diﬀerence discretization however does notsuﬃce for a faithful reproduction of ﬁeld line geometry and instead I have arguedthat a ﬁnite volume formulation using distributed sources is called for. To addressboth issues, I thus proposed a novel central ﬁnite diﬀerence discretization of thegauge action based on a stand-alone 2-by-2 plaquette centered around the point ofevaluation. Combined with forward and backward plaquettes on the boundary, asummation by parts discretized ﬁeld strength tensor ensues. I furthermore proposeto implement the distributed source structure by a modiﬁcation of the derivativeterms in the fermion action, which shifts the derivatives into adjacent planes. S × is easily deployed in existing simulation codes for Euclidean time. For its usein real-time applications, which are most conveniently formulated in a Hamiltonianlanguage, we may straight forwardly derive a continuous time variant. While theleap-frog type equations of motion from the naive Wilson action automaticallypreserve the diﬀerential Gauss-law up to rounding errors, the spatially extendednature of S × leads to a non-trivial consistency condition in the Dirac-Bergmansense, which has to be accommodated by a modiﬁed RATTLE algorithm [34], whichis work in progress.The novel action (interior 2 × × Acknowledgements

The author thanks Sigbjørn Hervik, Jan Nordstr¨om, Jan Pawlowski and AndersTranberg for valuable discussion and acknowledges funding from the Research Coun-cil of Norway under the FRIPRO Young Research Talent grant 286883. othkopf Page 8 of 9

Competing interests

The author declares that he has no competing interests.

Author’s contributions • A. Rothkopf: conceptualization, implementation, writing

References

1. Detmold, W., Kronfeld, A., Meiß ner, U.-G.: Topical issue on opportunities for lattice gauge theory in the era ofexascale computing. The European Physical Journal A (11), 192 (2019). doi:10.1140/epja/i2019-12942-8.Accessed 2021-02-142. [Extended Twisted Mass Collaboration], Alexandrou, C., et al. : Unpolarized and Helicity Generalized PartonDistributions of the Proton within Lattice QCD. Physical Review Letters (26), 262001 (2020).doi:10.1103/PhysRevLett.125.262001. Accessed 2021-02-143. Borsanyi, S., et al. : Ab initio calculation of the neutron-proton mass diﬀerence. Science (6229), 1452–1455(2015). doi:10.1126/science.12570504. Brice˜no, R.A., Dudek, J.J., Young, R.D.: Scattering processes and resonances from lattice QCD. Reviews ofModern Physics (2), 025001 (2018). doi:10.1103/RevModPhys.90.025001. Accessed 2021-02-145. Bazavov, A., Karsch, F., Mukherjee, S., Petreczky, P., USQCD Collaboration: Hot-dense lattice qcd. TheEuropean Physical Journal A (11), 194 (2019). doi:10.1140/epja/i2019-12922-0. Accessed 2021-02-146. Borsanyi, S., et al. : Calculation of the axion mass based on high-temperature lattice quantum chromodynamics.Nature (7627), 69–71 (2016). doi:10.1038/nature20115. Accessed 2021-02-147. Kockum, F., et al. : Ultrastrong coupling between light and matter. Nature Reviews Physics (1), 19–40 (2019).doi:10.1038/s42254-018-0006-2. Accessed 2021-02-148. Forn-D´ıaz, P., Lamata, L., Rico, E., Kono, J., Solano, E.: Ultrastrong coupling regimes of light-matterinteraction. Reviews of Modern Physics (2), 025005 (2019). doi:10.1103/RevModPhys.91.025005. Accessed2021-02-149. Jaako, T., Garcia-Ripoll, J.J., Rabl, P.: Quantum Simulation of Non-Perturbative Cavity QED with TrappedIons. Advanced Quantum Technologies (4), 1900125 (2020). doi:10.1002/qute.201900125. Accessed2021-02-1410. Nagle, J.L., Zajc, W.A.: Small System Collectivity in Relativistic Hadronic and Nuclear Collisions. AnnualReview of Nuclear and Particle Science (1), 211–235 (2018). doi:10.1146/annurev-nucl-101916-123209.Accessed 2021-02-1411. Mogliacci, S., Kolbe, I., Horowitz, W.: Geometrically conﬁned thermal ﬁeld theory: Finite size corrections andphase transitions. Physical Review D (11), 116017 (2020). doi:10.1103/PhysRevD.102.116017. Accessed2021-02-1412. Panero, M.: Geometric eﬀects in lattice QCD thermodynamics. PoS LATTICE2008 , 175 (2008).doi:10.22323/1.066.017513. Casanova, J., et al. : Deep Strong Coupling Regime of the Jaynes-Cummings Model. Physical Review Letters (26), 263603 (2010). doi:10.1103/PhysRevLett.105.263603. Accessed 2021-02-1414. Wilson, K.G.: Conﬁnement of quarks. Physical Review D (8), 2445–2459 (1974).doi:10.1103/PhysRevD.10.2445. Accessed 2021-02-1415. Klassen, T.R.: The anisotropic Wilson gauge action. Nuclear Physics B (1), 557–575 (1998).doi:10.1016/S0550-3213(98)00510-0. Accessed 2021-02-1416. Symanzik, K.: Continuum limit and improved action in lattice theories: (I). Principles and φ theory. NuclearPhysics B (1), 187–204 (1983). doi:10.1016/0550-3213(83)90468-6. Accessed 2021-02-1417. Garcia Perez, M., Gonzalez-Arroyo, A., Snippe, J.R., van Baal, P.: Instantons from over - improved cooling.Nucl. Phys. B , 535–552 (1994). doi:10.1016/0550-3213(94)90631-918. Beinlich, B., Karsch, F., Laermann, E.: Improved actions for QCD thermodynamics on the lattice. Nucl. Phys.B , 415–436 (1996). doi:10.1016/0550-3213(95)00681-819. de Forcrand, P., Garcia Perez, M., Stamatescu, I.-O.: Topology of the SU(2) vacuum: A Lattice study usingimproved cooling. Nucl. Phys. B , 409–449 (1997). doi:10.1016/S0550-3213(97)00275-720. Snippe, J.R.: Computation of the one loop Symanzik coeﬃcients for the square action. Nucl. Phys. B ,347–396 (1997). doi:10.1016/S0550-3213(97)00270-821. Bilson-Thompson, S.O., Leinweber, D.B., Williams, A.G.: Highly improved lattice ﬁeld strength tensor. AnnalsPhys. , 1–21 (2003). doi:10.1016/S0003-4916(03)00009-522. Langfeld, K.: A Novel improved action for SU(3) lattice gauge theory (2004). hep-lat/040301823. Jansen, K.: Lattice QCD: A Critical status report. PoS LATTICE2008 , 010 (2008). doi:10.22323/1.066.0010.0810.563424. Sv¨ard, M., Nordstr¨om, J.: Review of summation-by-parts schemes for initial–boundary-value problems. Journalof Computational Physics , 17–38 (2014)25. Fern´andez, D.C.D.R., Hicken, J.E., Zingg, D.W.: Review of summation-by-parts operators with simultaneousapproximation terms for the numerical solution of partial diﬀerential equations. Computers & Fluids ,171–196 (2014)26. Yamamoto, A.: Lattice QCD in curved spacetimes. Phys. Rev. D (5), 054510 (2014).doi:10.1103/PhysRevD.90.054510. 1405.666527. Jackson, J.D.: Classical Electrodynamics. Wiley, ??? (1998)28. Taﬂove, A., Hagness, S.C.: Computational Electrodynamics: The Finite-diﬀerence Time-domain Method.Artech House antennas and propagation library. Artech House, ??? (2005).https://books.google.no/books?id=n2ViQgAACAAJ29. Bruce Langdon, A.: On enforcing Gauss’ law in electromagnetic particle-in-cell codes. Computer PhysicsCommunications (3), 447–450 (1992). doi:10.1016/0010-4655(92)90105-8. Accessed 2021-02-1430. LeVeque, R.J.: Finite volume methods for hyperbolic problems (2002). doi:10.1017/CBO9780511791253 othkopf Page 9 of 9

31. Brambilla, N., Pineda, A., Soto, J., Vairo, A.: Eﬀective-ﬁeld theories for heavy quarkonium. Reviews of ModernPhysics (4), 1423–1496 (2005). doi:10.1103/RevModPhys.77.1423. Accessed 2021-02-1432. ˚Alund, O., et al. : Trace preserving quantum dynamics using a novel reparametrization-neutralsummation-by-parts diﬀerence operator. Journal of Computational Physics , 109917 (2021).doi:10.1016/j.jcp.2020.109917. Accessed 2021-02-1433. Kasper, V., Hebenstreit, F., Berges, J.: Fermion production from real-time lattice gauge theory in theclassical-statistical regime. Physical Review D (2), 025016 (2014). doi:10.1103/PhysRevD.90.025016.Accessed 2020-12-0634. Andersen, H.C.: Rattle: A “velocity” version of the shake algorithm for molecular dynamics calculations. Journalof Computational Physics (1), 24–34 (1983). doi:10.1016/0021-9991(83)90014-1. Accessed 2021-02-1435. Christiansen, S.H., Hu, K.: Finite Element Systems for vector bundles : elasticity and curvature.arXiv:1906.09128 [cs, math] (2020). Accessed 2021-02-1236. Schubel, M.D.: Discretization of diﬀerential geometry for computational gauge theory. PhD thesis, University ofIllinois at Urbana-Champaign (2018)37. Christiansen, S.H., Halvorsen, T.G.: Second order gauge invariant discretizations to the Schr \\