Gauge invariant canonical symplectic algorithms for real-time lattice strong-field quantum electrodynamics

A class of high-order canonical symplectic structure-preserving geometric algorithms are developed for high-quality simulations of the quantized Dirac-Maxwell theory based strong-field quantum electrodynamics (SFQED) and relativistic quantum plasmas (RQP) phenomena. With minimal coupling, the Lagrangian density of an interacting bispinor-gauge fields theory is constructed in a conjugate real fields form. The canonical symplectic form and canonical equations of this field theory are obtained by the general Hamilton’s principle on cotangent bundle. Based on discrete exterior calculus, the gauge field components are discreted to form a cochain complex, and the bispinor components are naturally discreted on a staggered dual lattice as combinations of differential forms. With pull-back and push-forward gauge covariant derivatives, the discrete action is gauge invariant. A well-defined discrete canonical Poisson bracket generates a semi-discrete lattice canonical field theory (LCFT), which admits the canonical symplectic form, unitary property, gauge symmetry and discrete Poincaré subgroup, which are good approximations of the original continuous geometric structures. The Hamiltonian splitting method, Cayley transformation and symmetric composition technique are introduced to construct a class of high-order numerical schemes for the semi-discrete LCFT. These schemes involve two degenerate fermion flavors and are locally unconditional stable, which also preserve the geometric structures. Admitting Nielsen-Ninomiya theorem, the continuous chiral symmetry is partially broken on the lattice. As an extension, a pair of discrete chiral operators are introduced to reconstruct the lattice chirality. Equipped with statistically quantization-equivalent ensemble models of the Dirac vacuum and non-trivial plasma backgrounds, the schemes are expected to have excellent performance in secular simulations of relativistic quantum effects, where the numerical errors of conserved quantities are well bounded by very small values without coherent accumulation. The algorithms are verified in detail by numerical energy spectra. Real-time LCFT simulations are successfully implemented for the nonlinear Schwinger mechanism induced e-e+ pairs creation and vacuum Kerr effect, where the nonlinear and non-perturbative features captured by the solutions provide a complete strong-field physical picture in a very wide range, which open a new door toward high-quality simulations in SFQED and RQP fields.


Introduction
Quantum electrodynamics (QED) at extreme conditions is becoming more and more important, as the relativistic quantum effects are becoming dominant mechanism in many branches of modern physics. With the development of high power laser technology, e.g. chirped pulse amplification (CPA), the peak intensity above 10 22 W·cm −2 is available by 1∼10 PW lasers, which is far stronger than the direct ionization threshold of 10 16 ∼ 10 18 W·cm −2 [1,2]. When the matter is exposed in such intense laser beams, strong ionization can be generated and large relativistic quantum plasmas (RQP) will be produced [3][4][5][6][7][8]. Next generation 10∼PW laser projects, such as the extreme light infrastructure (ELI) and the high power laser energy research facility (HiPER), aimed to approach the Schwinger threshold of 10 29 W·cm −2 or realize the fast ignition [1,2,5]. The Schwinger mechanism induced creation and following annihilation of fermion pairs play a fundamentally important role in modern high energy density physics (HEDP), astrophysics, and strong-field JHEP02(2021)127 quantum electrodynamics (SFQED) [3][4][5][7][8][9][10]. Although the direct experimental verification of electron-positron (e-e + ) pair creation under the Schwinger limit in laboratory is still expected to realize in near future, the e-e + RQP is already an important target for astronomical observers, such as the magnetosphere of an X-ray pulsar [11][12][13]. The typical magnetic field of X-ray pulsars is 10 12 G, and the effective temperature of X-ray pulsars is 10 KeV. In such a environment, the magnetic energy approaches to the rest energies of electron and positron, and it is higher than the thermal energy. As a result, the relativistic quantum effects lead to anharmonic cyclotron absorption features observed in spectra of X-ray pulsars [11][12][13][14]. Effective and accurate non-perturbative methods are needed in understanding these SFQED and RQP phenomena. Among a group of semi-analytical and numerical methods, the lattice quantum field theory (LQFT) is an advanced theoretical tool to study relativistic quantum effects both in vacuum and plasmas.
As a quantum gauge field theory on the discrete lattice in Euclidean space-time, the LQFT first developed by Wilson has been widely used in quantum chromodynamics (QCD) to describe the strong interactions, such as the quark confinement and the quark-gluon plasmas (QGP) [15][16][17][18][19]. Based on numerical path integrals and large-scale Monte Carlo (MC) simulations, the lattice quantum chromodynamics (LQCD) brings many significant results, such as the QCD phase transition and the hadron spectroscopy [18,19]. By using the Schwinger-Keldysh time contours, the LQCD can even be expanded to simulate non-equilibrium statistical systems [20,21]. The LQCD can not only capture the basic quantum loop effects, but also provide us with a well-defined non-perturbative theory of QCD [19]. As a post-MC technique, the tensor network (TN) techniques provide an alternative approach to simulate the lattice gauge theories (LGT), which can be efficiently extended to real-time evolution of out-of-equilibrium systems [22][23][24][25][26]. When it comes to phenomena with high occupation numbers and weak coupling, e.g. SFQED with non-trivial backgrounds and RQP, the classical relativistic field equations can be treated as good approximations to describe the dynamics of particles, where the quantum fluctuations can be introduced by constructing a statistically quantization-equivalent ensemble [27][28][29][30][31][32][33][34][35][36][37][38][39][40][41]. Based on this real-time lattice quantum electrodynamics (LQED) method in classical statistic regime, some interesting phenomena have been numerically studied, such as the pair creation of fermions beyond the Schwinger limit [33], the real-time dynamics of string breaking [32], the chiral magnetic effects [37], and the e-e + pair production in laser-plasma interactions (LPI) [40]. The LGT simulations are even reconstructed to implement on optical lattice based quantum simulators in recent time, which show great vitality [42,43].
When implementing a real-time LQED simulation in classical statistic regime, a stable and high fidelity numerical algorithm is needed to obtain reliable and accurate results. When it comes to the U(1) gauge field, there are many popular schemes for Maxwell's equations, such as the finite-difference time-domain (FDTD) method, the finite-element (FE) method, and the method of moments (MoM) etc., which are widely used in computational electrodynamics (CED) [44][45][46][47][48]. When it comes to the bispinor field, the numerical calculations of Dirac equation may encounter more difficult, e.g. fermion doubling problem, which will bring pseudo-fermion modes on the lattice [49]. There are several stable Dirac solvers, such as the time-splitting spectral (TSSM) method, quantum lattice Boltz-JHEP02(2021)127 mann (QLBM) technique, summation-by-parts-simultaneous approximation term (SBP-SAT) method, and time-dependent Galerkin (TDG) method etc. [50][51][52][53][54][55][56][57][58][59][60][61]. Although these algorithms have different advantages in part, an unified scheme with almost perfect performance is still a beautiful goal. Because of the nonlinearity and the multi-scale nature of the Dirac-Maxwell equations, high-quality simulations of real-time LQED face challenges. For instance, the numerical errors of conserved quantities can coherently accumulate, though these errors may be very small in each numerical step. The breakdown of conservation laws over a long simulation time amounts to pseudophysics. The structure-preserving geometric algorithms first developed by Feng et al. for classical Hamiltonian systems have excellent performance in long-term simulations [62][63][64][65][66][67][68][69][70][71][72], which are widely used in many complex systems, especially in geophysics and plasma physics [73][74][75][76][77][78][79][80][81][82][83][84][85][86]. In this work, we construct a class of structure-preserving geometric algorithms for Dirac-Maxwell theory. The algorithms preserve the symplectic and unitary structures, which also admit the U(1) gauge symmetry. The continuous Poincaré symmetry is reduced to a discrete subgroup, and there are only two degenerate fermion flavors exist on the lattice. The algorithms provide a powerful numerical tool for Dirac-Maxwell theory based real-time LQED simulations.
In section 2, a canonical field theory of Dirac-Maxwell systems is constructed in a conjugate real fields form to describe the fermion-photon interactions. The canonical symplectic form on cotangent bundle is obtained explicitly. In section 3, a semi-discrete lattice canonical field theory (LCFT) is constructed via the discrete exterior calculus (DEC) [87][88][89], the pull-back and push-forward gauge covariant derivatives, and a discrete canonical Poisson bracket. The LCFT admits symplectic and unitary structures, and also admits U(1) gauge symmetry on the lattice. In section 4, a class of high-order structure-preserving geometric algorithms are constructed for the LCFT. These schemes are locally unconditional stable, which also preserve the geometric structures and symmetries. The fermion doubling and chirality problems are discussed in detail. An ensemble model based field quantization procedure is reconstructed to simulate the Dirac vacuum and RQP. In section 5, the algorithms are detailedly verified and successfully used to simulate Schwinger mechanism induced e-e + pairs creation and vacuum Kerr effect. Numerical results show good properties in secular simulations. In section 6, we give a brief discussion about the advantages and attentions in implementing a real-time LCFT simulation by our algorithms, and show the outlook of applications in SFQED and RQP researches.

JHEP02(2021)127
Where the Dirac bispinor ψ is a 4 components complex field. With the Minkowski metric g µν = g µν = diag(+, −, −, −), the contravariant 4-vectors of coordinate and U(1) gauge field can be given by x µ = (ct, x), A µ = (φ, A), and the Maxwell gauge field strength tensor The Dirac conjugate bispinorψ = ψ + γ 0 , where the superscript + means Hermitian and the 4 × 4 matrices γ µ belong to a Clifford algebra {γ µ , γ ν } = 2g µν . With minimal coupling, the gauge covariant derivative D µ = ∂ µ + i e c A µ , and the Feynman dagger / D = γ µ D µ means Dirac contraction, where the charge e, reduced Planck constant and light speed c have their usual meanings. In the Dirac representation, the γ µ matrices are given by [93,94], Where the Pauli matricesσ i have their usual forms, The Dirac bispinor and it's Hermitian can be rewritten in a conjugate real fields form as, (2.7) By introducing the Dirac matrices [93,94], we can rewrite the Lagrangian density (2.1) as, Where E = −Ȧ/c − φ and B = × A are electric and magnetic strengths of the U(1) gauge field, and the superscript · means derivative with respect to time.Ĥ =Ĥ R + iĤ I = −i cα · − eα · A + eφ + βmc 2 is the Hamiltonian operator of the Dirac equation, which JHEP02(2021)127 can be given by, (2.11) By substituting eqs. (2.10)-(2.11) into eq. (2.9) and integrating it in full Minkowski spacetime manifold, we obtain the action functional S = T Ldt = T V Ld 4 x. Where the Lagrangian functional can be given by, On the tangent bundle T G of the configuration manifold G = (ψ R , ψ I , A, φ), the Hamilton's principle δS = 0 gives rise to the classical dynamical equations of the fields, (2.14) Where the bilinear form J = ecψ + αψ is the Dirac current density [94], and the gauge field components are restricted by the Gauss's law c 2 φ + ·Ȧ = −4πecψ + ψ.

Hamiltonian field theory
The cotangent bundle of the configuration manifold G can be defined as T * G = (ψ R , ψ I , A, φ, δL/δψ R , δL/δψ I , δL/δȦ, δL/δφ), where the variational derivatives of Lagrangian functional with respect to field components are given by,

JHEP02(2021)127
The Hamiltonian functional H can be obtained via the Legendre transformation of L, which is a map T G → T * G, On the cotangent bundle T * G, we can obtain a 2-form field, Where F and G are arbitary functionals on T * G. The dynamical equations of a field theory with canonical symplectic structure can be generated by the Hamiltonian functional of this field theory,Ḟ = {F, H} . (2.20) By taking the total variation of the Hamiltonian functional (2.17) with fixed boundary, we obtain,

JHEP02(2021)127
Where the Dirac current density can be expanded as,

JHEP02(2021)127
The canonical equations (2.23)-(2.32) equal to the dynamical equations (2.13)-(2.14), which means that the Hamiltonian field theory on T * G is an equivalent theory to the Lagrangian field theory on T G, both of which describe the intrinsic geometric structures of the interacting particles.

Gauge and Poincaré invariances
The QED is U(1) gauge and Poincaré invariant. Based on the Noether's theorem, the canonical field theory constructed in section 2 for Dirac-Maxwell systems admits charge, energy-momentum and angular momentum conservation laws [90,91].
The U(1) gauge symmetry means that the action and dynamical equations are invariant under the U(1) gauge transformation, where the gauge parameter θ is an arbitrary scalar field. It is convenient to verify the U(1) gauge symmetry of the canonical field theory by substituting eq. (2.33) into the Lagrangian density (2.9) or canonical equations (2.23)-(2.32). With an infinitesimal gauge transformation δ(A µ , ψ, ψ + ) = (∂ µ θ, i e c θψ, −i e c θψ + ) and using the dynamical equations (2.13)-(2.14), we can obtain the charge conservation law via δS = 0, which can be explicitly written as, The Poincaré symmetry consists of two parts, which are translation symmetry and Lorentz covariance. The translation symmetry means that the action and dynamical equations are invariant under the space-time translation, Where µ is an arbitrary translation parameter, and the local field components admit . By substituting an infinitesimal translations δ(x µ , A µ , ψ, ψ + ) = ( µ , 0, 0, 0) into δS = 0 and using the dynamical equations (2.13)-(2.14), we can obtain the energy-momentum conservation law as, Where the energy-momentum tensor is given by, The Lorentz covariance means that the action and dynamical equations are invariant under the Lorentz transformation,

DEC based discretization
The first step to construct a lattice field theory is discretization. As a differential geometry based numerical framework, DEC defines a class of complete operational rules and differential forms on a discrete differential manifold, which form a cochain complex [87][88][89].
To construct a semi-discrete LCFT for Dirac-Maxwell systems, the space-like submanifold of the Minkowski space-time manifold is discretized by using a rectangular lattice (other lattices are also viable). Then the scalar field A 0 = φ, which is a 0-form on the space-like submanifold, naturally lives on the vertex of the lattice, Where the subscript J indicates lattice label which traverses all lattice points, and (x i , y j , z k ) is the coordinate of the lattice vertex. The U(1) gauge and electric field 1forms A = A i dx i and Y = Y i dx i naturally live along the edges of the lattice, In the above discretization, a half integer index indicates along which edge does the field resides, where ∆x, ∆y and ∆z are lattice periods. In DEC framework, the magnetic field 2-form dA lives on the face center of the lattice. By using the Hodge dual operator * , we obtain the discrete charge 3-form −d * Y and current 2-form d * dA on the volume and face centers of the dual lattice respectively. Where the coordinate of a form on dual lattice

JHEP02(2021)127
is translated by (∆x, ∆y, ∆z)/2 after the Hodge operation, which means the primary-dual lattice generated by the Hodge star is a staggered lattice. The discrete gradient d , curl d ×, and divergence d · operators in DEC framework can be defined as [85,87], When it comes to the bispinor field, the fermion doubling is a serious problem in LGT, especially in LQCD simulations. Nielsen-Ninomiya no-go theorem states that the discretization of the Dirac equation on a regular space lattice forbids a single chirally invariant fermion flavor without breaking one or more of the following assumptions: translation invariance, locality, and Hermiticity [49]. There are several strategies to solve this problem, e.g. Wilson's momentum-dependent mass term and Kogut-Susskind staggered fermion [19,95]. In particular, the Kähler fermion constructed by P. Becher and H. Joos is the first geometric theory based lattice fermion, which is proved equivalent to a staggered fermion [96]. Here, in the unified DEC framework, we treat the bispinor components as different differential forms on the space-like submanifold, which are ψ 1 , ψ 2 dx∧dz, ψ 3 dz, ψ 4 dx,ψ 1 dx∧dy,ψ 2 dy∧dz,ψ 3 dx∧dy∧dz andψ 4 dy, Where the double sampled bispinor components live on the vertex, edge, faces and volume centers of the lattice respectively. By using the Hodge dual operator, we find the dual relations ψ 1 * ψ 3 , ψ 2 * ψ 4 , ψ 3 * ψ 1 and ψ 4 * ψ 2 , which naturally generate a staggered checkerboard-like lattice. The DEC dual relations could lead to a lattice bispinor field which involves two degenerate flavors. The DEC based discretization of the LCFT on the space-like submanifold is shown in figure 1.
Based on the discretization of fields, we can construct a discrete Poisson bracket, which admits bilinearity, anticommutativity, product rule, and Jacobi identity. The fields can be JHEP02(2021)127 reconstructed as, Then, the variational derivative with respect to A is [81,83], (3.15) and the variational derivatives with respect to Y , ψ iR /ψ iR and ψ iI /ψ iI have similar expressions. Here, V = x ∧ y ∧ z is the volume form on the lattice. Based on eq. (3.15), the canonical Poisson bracket (2.19) is discretized as, With the discrete canonical Poisson bracket (3.16), the functionals on the discrete cotangent bundle form a complete Poisson algebra. Then a semi-discrete LCFT can be generated by this discrete canonical Poisson bracket with a proper Hamiltonian functional on T * G d .

Pull-back and push-forward gauge covariant derivatives
The guage 1-form defines the guage connection on the U(1) bundle, which enables parallel transport bispinor on the Minkowski manifold. In order to construct a gauge invariant semi-discrete LCFT, we introduce a pair of discrete gauge covariant derivatives for different bispinor components, which can be recognized as Wilson lines in the DEC framework [40].
, and D zψ4 , the pull-back gauge covariant derivative D < is used along the relative gauge connections, e.g., The other D < components can be given in a similar form.
, and D yψ4 , the push-forward gauge covariant derivative D > is used along the relative gauge connections, e.g., The other D > components can be given in a similar form. By using d operator, the semi-discrete gauge transformation can be defined as,

JHEP02(2021)127
It shows that after a gauge transformation, the pull-back and push-forward gauge covariant derivatives get an unified phase, which ensures the semi-discrete Lagrangian density of the bispinor is gauge invariant.
The semi-discrete Lagrangian density of the U(1) gauge field is also gauge invariant in the DEC framework, which can be directly verified. As a result, the semi-discrete action functional admits gauge symmetry.

Semi-discrete canonical field theory
With the DEC and discrete gauge covariant derivatives, the Hamiltonian functional (2.17) is discreted as, (3.24) Where the superscript means 1-bispinor momentum, 2-bispinor mass-energy, and 3-U(1) gauge field respectively. The discrete Hamiltonian functionals are given by, By substituting the discrete Hamiltonian functional (3.24) into the discrete Poisson bracket (3.16), we obtain the canonical equations of the semi-discrete LCFT for Dirac-Maxwell systems. Here, we introduce the Hamiltonian splitting method and generate three linear canonical subsystems [83].

43)
(3.47) Where the Wilson line components cos / sin x J are defined as, In these equations, the translations of lattice index J ± 1 along the relevant gauge connections. The subsystem generated by H (2) d is given by, The second subsystem can be solved exactly when φ J (t) is given explicitly.

JHEP02(2021)127
The subsystem generated by H (3) d is given by, Because the semi-discrete action functional admits U(1) gauge symmetry, the semidiscrete LCFT for Dirac-Maxwell systems is gauge invariant. Because there are no explicit lattice coordinates in the semi-discrete action functional, it admits translation symmetry, then the semi-discrete LCFT for Dirac-Maxwell systems is translation invariant. We should emphasize that when the DEC lattice is fixed, the Lorentz boosts are forbidden, and there are only parity, time-reversal, and discrete rotations with angles (lπ/2, mπ/2, nπ/2) exsit in the discrete subgroup of SO(3, 1).

Gauge invariant canonical symplectic algorithms
Based on the Hamiltonian splitting method, the three linear canonical subsystems can be solved independently. The solution maps of the subsystems will be combined in various ways to give desired structure-preserving geometric algorithms for the semi-discrete LCFT.
For the subsystem generated by H (1) d , the canonical equations (3.28)-(3.47) can be rewritten as, Where Ξ(A) is an skew-symmetric matrix, which is also an infinitesimal generator of the symplectic group. To preserve the unitary property of the bispinor field, we adopt the symplectic mid-point method for this subsystem, and the one step map M D (∆t) : can be given by, Eq. (4.3) is a linear algebraic equation whose solution can be written as, Where Cay(S) denotes the Cayley transformation. It is well-known that Cay(S) is a symplectic rotation transformation when S in the Lie algebra of the symplectic rotation group [69]. As a result, the one step solution map M D (∆t) : is symplectic and unitary for bispinor field. Once ψ n+1 and ψ n+1 are known, Y n+1 can be calculated explicitly via eq. (4.5). Thus, is a second order symplectic scheme, which also preserves the unitariness of the bispinor field. For the subsystem generated by H (2) d , the canonical equations (3.51)-(3.60) can be solved exactly when φ J (t) is given explicitly. Because the LCFT is gauge invariant, we can get an explicit φ J (t) by adopting some gauge conditions, such as the temporal gauge φ = 0. Here, the one step solution map M M (∆t) : (ψ iR ,ψ iR , A, ψ iI ,ψ iI , Y ) n → (ψ iR ,ψ iR , A, ψ iI ,ψ iI , Y ) n+1 can be given by, (4.16) Where ω ± J (τ ) = [eφ J (τ ) ± mc 2 ]/ are eigen-frequencies of the fermion and anti-fermion. Because eqs.
is a symplectic scheme, which also preserves the unitariness of the bispinor field.
For the subsystem generated by H d , the canonical equations (3.61)-(3.64) can be rewritten as, Where Q is a constant matrix which belongs to the Lie algebra of the symplectic group. We also use the second order symplectic mid-point rule for this subsystem, and the one step (4.20) (4.21) A second order symplectic symmetric method can be constructed by the following symmetric composition, From a 2l-th order symplectic symmetric method M 2l (∆t), a 2(l + 1)-th order symplectic symmetric method can be constructed as [69], Obviously, the high order algorithms for the LCFT are symplectic and unitary structurepreserving, which is also gauge invariant.

Doubler, chirality and numerical stability
Fermion doubling problem is induced by a bad numerical dispersion of the free Dirac equation, which introduces pseudo cones in the Brillouin zone (BZ) of the lattice [53,60]. As a result, pseudo-fermion modes are excited on the lattice and the fermion velocity can faster than c. We can prove that the nonphysical doublers are suppressed in this work and there are only two degenerate fermion flavors exist on the DEC lattice. There are several equivalent approaches to achieve the doubler modes, such as the Poincaré-Hopf theorem based topological methd used for searching the chiral fermion doublers [97]. Here we directly derive the numerical dispersion or inverse Dirac propagator of the lattice fermion constructed in this work. A plane wave mode of the bipsinor field on the DEC lattice can be given by, (4.25)

JHEP02(2021)127
Where the superscript ∨ means amplitude. By substituting eq. (4.25) into the H (1) d subsystem, and setting the U(1) gauge field A µ = 0, we obtain the numerical dispersion of a mass free fermion as, (4.26) Once (∆t, ∆x) → (0, 0), the continuum limit of eq. (4.26) can be obtained as, (4.28) Where the momentum k is restricted in the lattice BZ k ∈ [−π/∆x, π/∆x]. Eq. (4.28) is also the numerical dispersion of the staggered fermion [98]. It shows that there is only one cone centered at k = 0 in the lattice BZ, and the Kramers and time-reversal symmetries are preserved. We should emphasis that due to the Dirac field is double sampled into eight complex field components, there are two degenerate fermion flavors, and the doubler degree is two. For free Dirac field these two flavors are uncoupled. When the gauge field is non-trivial, flavor mixing emerges, as the two flavors interact with different-valued gauge field on the DEC lattice. In this work, the strong gauge field will be treated as a classical field without fluctuation, then the flavor mixing is tolerable. The contribution of flavor mixing to gauge field is reflected in the ψψ andψψ terms in eqs. (3.45)-(3.47). The Nielsen-Ninomiya no-go theorem forbids perfect fermions on a regular lattice. The reduce of fermion doublers in our work is achieved by giving up on partial chiral symmetry. In the Dirac representation, the chiral operator can be given by γ 5 = iγ 0 γ 1 γ 2 γ 3 , and the chiral symmetry of free fermions can be defined by an anti-commutator { / D, γ 5 } = 0. Based on the Hamiltonian (3.25) of semi-discrete lattice field theory, the spatial Dirac contraction of discrete gauge covariant derivative is obtained as, (4.29) -21 -

JHEP02(2021)127
Then we can derive the chiral anti-commutator with discrete derivative as, The residual term in eq. (4.30) shows that the break of chiral symmetry is only induced by the non-balance between z-direction discrete gauge covariant derivatives. This partial chirality is gauge dependent and the complete chiral symmetry can naturally recover in the continuous limit. We indicate that the broken symmetry comes from the non-commutation between discrete derivative and continuous chiral operator, and we can define a complete discrete chirality on the DEC lattice by introducing a pair of discrete chiral operators, Where the shifting operators CS + translates lattice indexes of the 2nd, 3rd, 6th and 7th bispinor components as J →J + 1 along z-direction, and CS − translates lattice indexes of the 1st, 4th, 5th and 8th bispinor components as J →J − 1 along z-direction. Then Γ ± are inverse operators. By using the discrete chiral operators, we can define the discrete chiral fermions as, Where the chiral projection operators are given by P L = 1 2 (1 − Γ + ) and P R = 1 2 (1 + Γ − ). By substituting eqs. (4.32)-(4.33) into eq. (3.25), we obtain, Then the discrete dynamical equations of free chiral fermions are obtained by the discrete variational principle, We should emphasize that the complete discrete chirality is only a lattice analogue, and the free chiral fermions defined by continuous chiral operator are coupled on the DEC lattice, where the nonphysical coupling will vanish in the continuous limit. The numerical dispersion (4.26) shows that the energy ω is always real for arbitrary lattice periods (∆t, ∆x), which means that the solution map M D (∆t) of H well-known unconditional stable CED scheme [47]. As a result, the high order solution map M (∆t) for the LCFT is locally unconditional stable.
To implement the algorithms, Jacobian inversions are needed, as the Cayley transformation brings several linear algebraic equations. The Krylov subspace theory provides us with many efficient linear solvers, such as the generalized minimum residual (GMRES) method, the incomplete Cholesky conjugate gradient (ICCG) method, and the biconjugate gradient stabilized (BICGSTAB) method, which can be used to solve the large sparse matrix equation. Based on these efficient linear solvers, the algorithms can be conveniently implemented via standard parallel strategies.

Field quantization
In classical statistic regime, the quantization of Dirac field can be simulated via a statistically quantization-equivalent ensemble, which can reconstruct pairs of anticommuting creation and annihilation operators in a statistic sense. The low-cost fermions strategy is such an ensemble model widely used in real-time LGT simulations [29,30]. Along this approach, we introduce a unified ensemble model of vacuum and plasmas into the LCFT to realize real-time LCFT simulations for abundant SFQED and RQP phenomena. Based on the standard canonical quantization procedure, the Dirac bispinor field can be quantized as [90,91,99], Where {a s p , a s + p } = (2π ) 3 δ(p−p )δ s,s and {b s p , b s + p } = (2π ) 3 δ(p−p )δ s,s are creation and annihilation operators for fermions and antifermions respectively, with spin index s = ± 1 2 . u s (p) and v s (p) are relevant eigen spinors of free particles, which can be normalized as, Where U s are normalized Pauli spinors U s+ U s = δ s,s . E is the positive definite energy norm of an on shell fermion. Then the orthogonality relations u + s (p)u s (p) = v + s (p)v s (p) = δ s,s and u + s (p)v s (−p) = 0 can be obtained directly. To achieve the Fermi-Dirac statistics via a classical Dirac field ensemble, we replace the creation and annihilation operators with a class of stochastic variables and reconstruct a pair of stochastic Dirac spinors as, To describe the Dirac vacuum, we can assume these stochastic variables admit the same amplitude distribution ( (2π ) 3 − σ 2 , σ 2 ) and uniform phase distribution U(−π, π] in the momentum space. To describe a single specie and spin polarized plasma background, we can set distributions of the stochastic variables ξ(p) and η(p) admit the Pauli blocking density. Then the ensemble averaged bilinear covariant gives rise to the background plasma density. Eqs. (4.40)-(4.41) ensure the statistically equivalence between ensemble model and field quantization for the Dirac vacuum and non-trivial plasma backgrounds [30,40]. On a DEC lattice, the stochastic Dirac spinors are discreted as, Then the non-trivial correlators are where the discrete momentum space is given by p ∈ [−π /∆x, π /∆x] with lattice spacing 2π /∆xN x .
Based on this ensemble model, the Lagrangian density eq. (2.1) on T G = (ψ M R , ψ M I , ψ F R , ψ F I , A, φ,ψ M R ,ψ M I ,ψ F R ,ψ F I ,Ȧ,φ) can be rewritten as, Where the gender conjugate g.c. means commutation between the stochastic bispinor pairs, and : · : means normal product. Then the Hamiltonian functional on can be obtained via the Legendre transformation, (4.45) Ne ·, where the ensemble capacity Ne is the number of systems in a given ensemble. The minimum Ne should be larger than the lattice degree of freedom, which guarantees all the lattice modes can be sampled.

JHEP02(2021)127
To keep the physical constraints, the sampling of gauge field configuration should adapt the self-consistent field condition at initial time. We should emphasize that the gauge field can be treated as a classical field without quantum fluctuation only valid in very high occupation states. The SFQED and RQP phenomena always satisfy this condition. When it comes to weak field problems, e.g. some quantum optics and quantum electronics phenomena, the strong-field condition breaks and the quantization of gauge field should also be taken into consideration.

Energy spectra
To verify the canonical symplectic structure-preserving geometric algorithms constructed in this work, we implement the code to obtain a class of numerical energy spectra of the Dirac-Maxwell theory based LCFT. As benchmarks, the analytical dispersion relations of linearized scalar QED are introduced to compare with these numerical energy spectra [6,14]. The dispersion relation of free Dirac fermions are given by [90], which means that there are two fermion modes sharing a gap 2mc 2 / . The Dirac doublecone of positive and negative states is a basic property of relativistic particles. If there is no strong background magnetic field, e.g. the vacuum and unmagnetized plasmas, the Klein-Gordon-Maxwell (KGM) theory based scalar QED is a good toy model of the Dirac-Maxwell fields theory, both of which admit the same branches of linearized dispersion relations. The 1/2-spin effects only modify the mode structures. The tree-level electromagnetic mode dispersion relation of the scalar QED can be given by [14], Where ω p is the plasma frequency of background fermions. When it comes to the vacuum, ω p = 0 and eq. (5.2) reduces to the light cone. The tree-level dispersion relation of electrostatic mode can be given by [6], Eq. (5.3) shows that the electrostatic mode consists of four branches. In vacuum, two gapless branches relate to the fermions moving with self gauge fields, which are known as Langmuir modes. The other two gapped branches are pair modes, and the half gap 2mc 2 / means that if the photon energy ω > 2mc 2 , the fermion pairs will be generated and the quanta of these pair plasmas have finite group velocities. The pair mode is also known as Zitterbewegung effect in relativistic quantum mechanics, which is described as the interference between positive and negative states of a fermion on the Compton space-time scale.
To implement real-time LCFT simulations, the natural units are used, where the constants = c = 1, the elementary charge e = 0.0854, so that the fine structure constant JHEP02(2021)127 To calculate the energy spectra, a uniform 256 × 256 × 256(N x ×N y ×N z ) DEC lattice is introduced, where the spatial lattice periods ∆x = ∆y = ∆z = 1, and the temporal lattice period ∆t = 0.5. In all directions, the periodic boundary is used to naturally introduce an infrared truncation. To initialize the simulations so that a broad spectrum of linear waves are excited, the bispinor field is given using small amplitude unbiased white noise with standard deviation σ = 1 × 10 −6 , and the temporal gauge is adopted explicity. After a N t = 512 steps simulation, the numerical spectra of the LCFT can be read out from simulation results by taking multi-dimensional fast Fourier transforms (FFT) of bispinor and electric field components. By abandoning the gauge field in simulation, we can obtain the dispersion relation of a free fermion via the same procedure.  0-spin KGM theory. Figure 3 (d) shows a weak self gauge field dressed fermion mode, which can be well traced by the free fermion dispersion. That the analytical dispersion relations are recovered by numerical spectra indicates that our solutions faithfully capture the propagation of linear waves up to the lattice resolution.
To illustrate the advantages of our algorithms, we implement a long-term simulation and record the numerical errors of the conserved quantities. The simulation domain is a uniform 10 × 10 × 10 DEC lattice, and the periodic boundary is used in all directions. ∆x = ∆y = ∆z = 1, and ∆t = 0.05. At initial time, a unbiased white noise with standard deviation σ = 3×10 −2 is introduced into the bispinor field. After a million steps simulation, the relative numerical errors of the total Hamiltonian and total probability are plotted in figure 4. We find that after a extremely long-term simulation, the numerical errors of conserved quantities are bounded by small values without coherent accumulation. The excellent conservation property in these numerical solutions comes from the preservation of geometric structures and symmetries by using our algorithms. The conservation is a footstone to implement secular simulations for nonlinear multi-scale SFQED and RQP phenomena.

Schwinger mechanism induced e-e + pairs creation
The Schwinger mechanism induced creation and annihilation of electron and positron pairs are genuine phenomena in SFQED, which can not be described via classical theories [100]. In section 5.1, the numerical spectra show that the pair mode can be found once the energy of the γ photon exceeds double electron rest energy. Schwinger effect states that when the photon wavelength is not very short, the e-e + pair can also be generated once the gauge field strength is extremely strong [10,90]. The typical electrostatic field strength of the Schwinger limit is E S = 1.32 × 10 16 V/cm, and the equivalent magnetic field strength JHEP02(2021)127 and laser intensity are of orders 10 9 T and 10 29 W·cm −2 respectively [10]. Beyond the Schwinger threshold, the virtual e-e + pairs can be pulled apart from quantum fluctuations on the Compton space-time scale and large on shell e-e + pairs can be created from the vacuum. Although there are some other QED mechanisms can create on shell e-e + pairs in the vacuum, e.g. the Breit-Wheeler process, the Schwinger effect becomes the dominate process once the U(1) gauge field becomes a low frequency and extremely strong field. The Feynman diagram of Schwinger effect is shown in figure 5.
To simulate the Schwinger mechanism induced e-e + pair creation, we set a longitudinal quasi-static electric field whose strength is normalized by the Schwinger limit E S = m 2 /e. Numerical experiments are implemented on a 1 × 1 × 256 uniform DEC lattice, and the periodic boundary is used in all directions. e/m = 0.2, ∆x = ∆y = ∆z = 0.05/m, and ∆t = 0.5∆x i . At initial time, an ensemble model based Dirac vacuum state is introduced by sampling a class of stochastic Dirac spinors ψ M (x, 0) and ψ F (x, 0) as n + p = n − p = 0, where the ensemble capacity N e = 512. The gauge field A µ is sampled in the temporal gauge, where A(x, 0) = 0 and Y (x, 0) = (0, 0, −E S /4π) are given as an initial condition. By setting Y (x, 0) = 0, we can also simulate the quantum fluctuations of Dirac vacuum. After a 20000 steps simulation, the numerical results are recorded, which include the electric field evolution, pair production rate, Hamiltonian transfer, and spectral density. Figure 6 illustrates the numerical evolution of normalized Hamiltonians, where the vacuum energy has been renormalized by the normal product. The blue solid line show us the decaying oscillation of the U(1) gauge field. During this process, on shell e-e + pairs are continuously created and driven, and the energy of photons is continuously transfered into the fermions energy. With the growth of fermion density, the pair plasma frequency increases and the chirp feature of the plasma oscillation can be distinctly recognized from the Hamiltonian of gauge field. This is a nonlinear phenomenon, as the production of e-e + pair can be effectively suppressed by the radiation reaction, which means that the energy of gauge field will be absorbed by self generated pair plasmas and the pair production will reach a saturation level. This nonlinearity can also be read out from the fermion Hamiltonian (red dashed line). The lower envelope of this curve told us that the production rate of the e-e + pair in the simulation domain has a very high level in the early time of the vacuum breakdown (t·m ∼ 50), and then the pair production rate is saturated with a relatively stable plasma frequency after a long time evolution (t·m ∼ 500). During this nonlinear Schwinger process, the total Hamiltonian in the simulation domain is perfectly JHEP02(2021)127 conserved (black solid line). In summary, the photon energy is continuously absorbed by new on shell e-e + pairs, and the energy between ptoton and pair plasmon is exchanged cycle by cycle with a chirped plasma frequency.
Moreover, we plot the guage connection A z , electric field E z , net charge Q and Dirac current density J z in figure 7 to demonstrate the dynamical properties of U(1) gauge field and fermions. The numerical evolution of gauge field shown in figure 7 (a)-(b) illustrates the dissipative anharmonic effects induced by the separation and recombination of nonlinear ee + pair oscillators. From the Hamiltonian curves plotted in figure 6, we already know that after hundreds of Compton periods 1/m, the pair production will be effectively suppressed by the nonpertubative field backreaction. In figure 7 (b), we find a slowly varying electric field amplitude E z ≈ 0.68E S at the end of this simulation. During this time, more than half of the gauge field energy has been consumed to create e-e + pairs, and the consumed energy is transferred and stored in the pair plasmas. Gauge symmetry induced charge conservation law can be found in figure  The chirped electric field exhibits dissipative anharmonic feature, which can be used as a probe to diagnose the state of e-e + pair oscillators. In this simulation (t ∼ 500/m), the electric field amplitude approaches 0.68E S , which means that more than half of the gauge field energy has been consumed to create e-e + pairs. All QED phenomena admit the charge conservation law, which means that the net charge Q in Schwinger process must keep 0. It can be found in subfigure (c), where a numerical noise induced extremely small net charge is well conserved in the simulation. (d) shows the pair plasma oscillation induced ensemble current density. With the nonlinear increase of fermion density, The current amplitude experienced a growth process with a decay rate. Obviously, the chirped oscillations of current and elctric field are well matched with a π/2 phase difference. (e/m = 0.2, ∆x i = 0.05/m, ∆t = 0.5∆x i , N t = 2 × 10 4 , N x = N y = 1, N z = 256, N e = 512). the oscillations of current and elctric field are well matched with a π/2 phase difference. Different from the classical oscillator, the nonlinear increase of pair plasma density gives rise to a nonlinear current amplitude growth and an oscillation frequency blue shift.
To compare the Schwinger effect with different background gauge field strengths, we implement a class of simulations under different initial electric fields (normalized by E S ), and plot the associated fermion Hamiltonian curves in figure 8. The lower envelope of the curve with E z = 0.6 demonstrate that the Schwinger effect can be effectively cut off when the background electric field is lower than half E S . In this situation, the extremely tenuous fermion density gives rise to a very low plasma frequency, and the energy transfer from photons to fermions is very slow. On the contrary, the lower envelope of the curve with E z = 1.4 shows us that the nonlinear suppression effect in Schwinger process can be significantly enhanced when the background electric field is far stronger than E S . Finally, we can make a brief summary that the Schwinger mechanism induced fermion pairs production is inherently a nonlinear and non-perturbative phenomenon, which exhibits abundant anharmonic, non-equilibrium and self-modulation features. To illustrate the complete physics of this process, real-time LGT simulation or other non-perturbative methods are needed. Due to the symmetric and geometric structure-preserving nature, Our algorithm provide an efficient, accurate, stable and conservative approach to implement real-time LCFT simulations to study this kind of complicated SFQED phenomena.

Vacuum Kerr effect
The Schwinger mechanism induced pair plasmas are strongly polarized, which means that the Dirac vacuum is strongly polarized under an extreme electric field. The Kerr effect states that when a dielectric medium is polarized by an external electric field, it will exhibit birefraction property, for the refractive index parallel to external electric field is modulated. Be treated as an QED analogue of the classical polarized dielectric medium, the polarized Dirac vacuum may also exhibit birefraction property, and a Kerr-like effect can be expected to be observed in the SFQED regime. A schematic of the vacuum Kerr effect is shown in figure 9.
To simulate the vacuum Kerr effect, we set a transverse quasi-static electric field whose strength approaches E S , and then introduce a weak linear polarized free electron laser (FEL) beam as incident wave. Numerical experiment is implemented on a 40 × 40 × 20 uniform DEC lattice, and the periodic boundary for both fermion and gauge field is used in x and y directions. To cut off the longitudinal radiations, we introduce the second order Mur's boundary in z direction to simulate a open space. e/m = 0.2, ∆x = ∆y = ∆z = 0.1/m, and ∆t = 0.5∆x i . At initial time, an ensemble model based pair plasma state is introduced by sampling the stochastic bispinors ψ M (x, 0) and ψ F (x, 0) as n + p = n − p = 0.1, where the ensemble capacity N e = 3200. The gauge field A µ is sampled in the temporal gauge, where A(x, 0) = 0 and Y (x, 0) = (0, 0, −E S /4π) are given as an initial condition. To excite an incident FEL plane wave, we set a total-scattered fields boundary in the z = 5∆z (source) plane, and the laser frequency is given by ω = 0.2πm. After a 400 steps simulation, the numerical results are recorded, where the trace of magnetic field vector in the z = 15∆z (target) plane demonstrates the polarization state of the beam.  found in figure 10 (b), where the trace of magnetic field vector in the source plane draws a perfect line segment. When propagating throw the strongly polarized Dirac vacuum, the x and y components of the FEL beam have different phase velocities. Then there is a phase difference between x and y modes in the target plane, and the beam polarization will transferred into a elliptical-like state, which can be found in figure 10 (d). From the numerical results of Schwinger mechanism induced e-e + pairs creation, we know that the pair plasma density and the background electric field strength are time dependent. As a result, the Vacuum Kerr effect admits a varying coefficient and the phase difference is time dependent. The insets of figure 10 (d) demonstrate this feature, where the second cycle (C2) of the beam admits a larger semi-major axis than the first cycle (C1), and the trace gyrocenter is drifting perpendicular to the major axis. If the pair plasma is dense enough, more significant polarization conversion features can be expected to be observed.
In summary, the numerical experiments illustrated in section 5 cover relativistic QED wave structures, fermion pairs creation and annihilation effects, self-consistent interactions between fermion plasmas and U(1) gauge field, nonlinear and non-perturbative nature of strong-field physics. The good properties of these numerical solutions ensure the structurepreserving real-time LCFT simulation method is expected to be a unified first-principle based theoretical tool in studying SFQED and RQP phenomena.

Conclusion and outlook
In this paper, we developed a class of high-order canonical symplectic structure-preserving geometric algorithms for simulating the quantized Dirac-Maxwell theory based SFQED and RQP. We constructed a canonical field theory of the Dirac-Maxwell systems, and obtained the canonical symplectic form and Poisson algebra admitted by this field theory. Based on the Noether's theorem, this field theory admits charge, energy-momentum JHEP02(2021)127 and angular momentum conservation laws via the gauge and Poincaré symmetries. In DEC framework, we constructed a LCFT which is a good semi-discrete analogue of the continuous canonical field theory. The U(1) gauge field is discreted to form a cochain complex which guarantees the Bianchi identities of the U(1) gauge theory. With the Hodge dual relations, the bispinor field components are discreted as eight different differential forms, which naturally generate a staggered checkerboard-like lattice. Two kinds of discrete gauge covariant derivatives, i.e. pull-back and push-forward, are used to construct a gauge invariant semi-discrete action. A well-defined discrete Poisson bracket is constructed, which admits bilinearity, anticommutativity, product rule, and Jacobi identity. With the previous numerical techniques, the semi-discrete LCFT is gauge invariant, which also preserves the canonical symplectic and unitary structures. By using the Hamiltonian splitting method, we obtained three linear subsystems which can be solved independently,

JHEP02(2021)127
and constructed a class of high-order structure-preserving geometric algorithms via the Cayley transformation and symmetric composition technique. The algorithms preserve the gauge symmetry and geometric structures of the semi-discrete LCFT. We proved that the numerical dispersion of the mass free fermions subsystem has only one Dirac double-cone centered at the origin of lattice BZ, and there are only two degenerate fermion flavors in our algorithms. The locally unconditional stable property can also be obtained from the numerical dispersion. The structure-preserving and unconditional stable properties make the scheme superior to conventional Wilson and staggered fermions. To simulate the quantization of Dirac field to achieve a correct Fermi-Dirac statistics, we introduce an unified statistically quantization-equivalent ensemble model to describe the Dirac vacuum and non-trivial plasma backgrounds. Although the algorithms are unconditional stable, it does not means that the lattice periods can be chosen arbitrary large values. On the one hand, some basic physics can not the captured if the lattice periods exceed the typical space-time scales, such as the Compton scale. On the other hand, large lattice periods will lead to sparse matrices with very large condition numbers, which are very expensive for matrix inversion. Additionally, the topology of the U(1) gauge field is changed into a torus by the Wilson lines, which means that the lattice periods should not be too large to avoid topological modes.
The numerical energy spectra of the LCFT were calculated and compared with the analytical dispersion relations of linearized scalar QED. Simulation results show that the relativistic quantum wave dynamics and the vacuum responses can be captured perfectly. The gapless lower branches (Langmuir) of electrostatic mode relate to the fermions moving with self gauge fields. The gaped higher branches of electrostatic mode are QED pair modes, where the virtual fermion pairs in quantum fluctuations are generated. As the quanta of pair plasmas, the pair plasmons have finite group velocities, which means the virtual pairs created and annihilated on the Compton space-time scale are very different from the classical plasmas. The nonlinear Schwinger effect was also simulated to illustrate the power of our algorithms. To simulate the quantum fluctuations, we introduced an ensemble of statistically quantization-equivalent initial conditions via random momentum and phase, which can be used as a statistical model of the quantized Dirac vacuum. With a uniform strong field, the pair production rate is obtained with a nonlinear suppression, which means that the energy of the gauge field will be absorbed by self generated pair plasmas and the pair creation will reach a saturation level. This nonlinear property of Schwinger mechanism can only be resolved by non-perturbative methods, such as the quantum particle-in-cell (PIC) and real-time LQED methods [5,30,101]. Our algorithms provide a more accurate and efficient solver for simulating these SFQED and RQP problems because of the advanced conservation performance in secular simulations and good unconditional stable property. We also simulated the vacuum Kerr effect, where the vacuum response can be resolved. After propagating through a strongly polarized vacuum area, a linear polarized FEL beam transferred into an elliptical polarized one, where the vacuum birefraction property was well traced. Because the vacuum polarization state is strongly affected by the evolution of background electric field, the birefraction index of polarized Dirac vacuum is time dependent. This dynamical property can only be resolved JHEP02(2021)127 by nonlinear non-perturbative methods. All simulations implemented in this work show a common property that the numerical errors of conserved quantities, e.g. total Hamiltonian and charge, are bounded by a very small value after a long-term simulation. This advantage enables us to simulate nonlinear multi-scale problems dominated by relativistic quantum effects, such as a high energy FEL beam interacting with RQP and the magnetosphere of an X-ray pulsar.
In summary, the gauge invariant canonical symplectic structure-preserving geometric algorithms constructed in this work provide us with a powerful first-principle based theoretical tool to implement quantized Dirac-Maxwell theory based real-time LCFT simulations. Because of the nonlinear and non-perturbative nature of this approach, it can bring abundant physics from the interacting fields. With well-designed field quantization models, this method opens a new door toward high-quality simulations in SFQED and RQP fields.