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 finite 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 finite differences as starting pointand combining it with insight from the finite 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 first 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 confined light fields have made accessible a regime ofnon-perturbatively large couplings in cavity quantum electrodynamics (QED) (forrecent reviews see [7, 8]), which defies the standard methods of perturbation theory.A lattice field 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 finite 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 finite extent (see e.g. [11, 12]). At low energies the confines of the reflectivecavities 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 fields, non-perturbative theory supportis essential, which however can only be provided by a proper lattice gauge theoryfor finite systems.According to Wilsons seminal work [14], the simplest gauge invariant action fornon-Abelian fields 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 finite differences (FD) are defined 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 significant 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 efficient Monte-Carlo sim-ulations of the corresponding Euclidean path integral [23]. Two challenges exist in deploying the discretization of eq. (1) in a finite system,where PBC are inapplicable and translational invariance is broken: (A) a naive firstorder FD scheme does not accurately reproduce the true solution of a differentialequation 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 finite circular capacitor,whose end-plates are located on the boundary in the z-direction, held at a finitepotential. 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 fig. 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 field 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 fig. 1, fieldstrength is generated close to the forward boundary but the field lines show anartificial staggered pattern.The correct treatment of non-trivial boundary conditions in FD approximationsis reflected 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 field (red arrows) from the capacitor potential (background) vs. backwardFD solution (blue arrows). (right) The field strength from a naive central FD (green arrows) andthe more accurate solution based on a SBP discretization (red arrows). must require the first order FD approximation to fulfill 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 significantly more accurate solution for the capacitor denoted by the redarrows in the right panel of fig. 1.(B) The second challenge pertains to the fact the naive forward FD approximationof eq. (1) introduces significant 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 field becomes ambiguousin non-Abelian gauge theory, where one must instead resort to computing the fieldstress 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 differentiation and S denotes the Poyntingvector, which is irrelevant in the static case [27]. The force field 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 fig. 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 find that in the yz plane only a few of its EV othkopf Page 4 of 9 z Figure 2 (left) Force field 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 fig. 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 field is not accurately reproduced.Solving Gauss’ law with the central FD discretization we find the symmetric greenarrows in the right panel, which however exhibit an artificial 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 finite volume discretization[30], i.e. one that obeys both the differential 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 finitedifference 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 finite volume discretization ofGauss’ law. Solving the above yields the symmetric and significantly more accuratesolution shown as red arrows on the right of fig. 2.We conclude from the Abelian Gauss law that a proper lattice gauge theoryfor finite 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 finite 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 find 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 fields inthe forward and backward direction in each dimension. The product of eight links(see fig. 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 finite 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 differenceto 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 finite difference.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 finally 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 fields. Let me demonstrate thestrategy for static matter fields, 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 modification 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 finitedifference 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 finite difference Gauss law, which as Ihave shown imprints asymmetries into the force field 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 finite difference discretization however does notsuffice for a faithful reproduction of field line geometry and instead I have arguedthat a finite volume formulation using distributed sources is called for. To addressboth issues, I thus proposed a novel central finite difference 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 field strength tensor ensues. I furthermore proposeto implement the distributed source structure by a modification 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 differential 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 modified 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 difference. 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 confined thermal field theory: Finite size corrections andphase transitions. Physical Review D (11), 116017 (2020). doi:10.1103/PhysRevD.102.116017. Accessed2021-02-1412. Panero, M.: Geometric effects 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.: Confinement 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 coefficients 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 field 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 differential 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. Taflove, A., Hagness, S.C.: Computational Electrodynamics: The Finite-difference 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.: Effective-field 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 difference 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 differential 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 \\