# Simulations of relativistic-quantum plasmas using real-time lattice scalar QED

## Abstract

Real-time lattice quantum electrodynamics (QED) provides a unique tool for simulating plasmas in the strong-field regime, where collective plasma scales are not well-separated from relativistic-quantum scales. As a toy model, we study scalar QED, which describes self-consistent interactions between charged bosons and electromagnetic fields. To solve this model on a computer, we first discretize the scalar-QED action on a lattice, in a way that respects geometric structures of exterior calculus and U(1)-gauge symmetry. The lattice scalar QED can then be solved, in the classical-statistics regime, by advancing an ensemble of statistically equivalent initial conditions in time, using classical field equations obtained by extremizing the discrete action. To demonstrate the capability of our numerical scheme, we apply it to two example problems. The first example is the propagation of linear waves, where we recover analytic wave dispersion relations using numerical spectrum. The second example is an intense laser interacting with a 1D plasma slab, where we demonstrate natural transition from wakefield acceleration to pair production when the wave amplitude exceeds the Schwinger threshold. Our real-time lattice scheme is fully explicit and respects local conservation laws, making it reliable for long-time dynamics. The algorithm is readily parallelized using domain decomposition, and the ensemble may be computed using quantum parallelism in the future.

## I Introduction

Lattice QED, a scheme usually used to study vacuum quantum electrodynamics, can also be used to simulate plasmas. By adding dynamical background fields, we extend lattice QED into a valuable tool for plasma physics, especially when plasmas are dense or when fields are strong. Under these extreme conditions where collective QED effects are important, the commonly adopted plasma kinetic model, which arises as the geometrical optics approximation of the relativistic-quantum world Ruiz and Dodin (2015), is no longer sufficient. An example where QED effects are important is the production of electron-positron pairs when intense lasers interact with plasma targets Liang *et al.* (1998); Gahn *et al.* (2000); Liang *et al.* (2015); Sarri *et al.* (2015). To describe such phenomena in semiclassical framework, source terms must be inserted into quantum kinetic or fluid equations Berezhiani *et al.* (1992); Kluger *et al.* (1998); Schmidt *et al.* (1998); Roberts *et al.* (2002); Hebenstreit *et al.* (2010), which can then be solved by numeric integration Hebenstreit *et al.* (2008, 2009) or QED-PIC simulations Duclous *et al.* (2011); Nerush *et al.* (2011); Ridgers *et al.* (2012). However, prefabricated source terms take little account of the interplay between coexisting processes Schützhold *et al.* (2008), nor can they describe quantum interference, through which the created pairs may be entangled. Therefore, while semiclassical approximations may be applicable for long-wavelength lasers, large errors are expected when fields, such as those of x-ray lasers, evolve on scales comparable to intrinsic QED scales. Moreover, in semiclassical treatments, there is no obvious way to subtract both energy and momentum from fields once pairs are produced. Although errors may be tolerable when fields dominate particles, disrespecting energy-momentum conservation will likely have nonphysical consequences, after a large number of pairs are generated.

To model plasmas where QED processes have no clear scale separation from classical processes, a faithful description can only be provided on the relativistic-quantum level. While lattice simulations may be unfamiliar tools for plasma physics, they have been used extensively in quantum chromodynamics (QCD) to describe the strong interaction Wilson (1974) and quark-gluon plasmas Bass *et al.* (1999); Satz (2000). In conventional lattice-QCD simulations, quantum correlation functions are computed using numerical path integrals, from which observables are extracted as coefficients of scaling laws Creutz (1980). This scheme can be analytically continued to imaginary time to describe statistical systems in thermal equilibrium Yagi *et al.* (2005). For out-of-equilibrium systems, real-time simulations can be carried out using the Schwinger-Keldysh time contours Schwinger (1961); Keldysh (1965). The above formulations, based on numerical path integrals, are capable of capturing genuine quantum loop effects, but are numerically expensive. Fortunately, the computational cost can be dramatically reduced when the occupation numbers of quantum states are high and when the coupling is weak. This is precisely the case for plasma physics, where a large number of particles are present, and the coupling coefficient is small. In this classical-statistic regime, tree-level effects dominate loop effects Aarts and Berges (2002); Mueller and Son (2004); Jeon (2005); Berges and Gasenzer (2007); Berges *et al.* (2014), and the quantum system can be adequately described by time-advancing the classical field equations with an ensemble of statistically equivalent initial conditions Aarts and Smit (1999); Gelis and Tanji (2013); Polkovnikov (2003); Borsanyi and Hindmarsh (2009). Based on this approach, lattice spinor-QED simulations have been carried out to demonstrate production of fermion pairs from the vacuum by a prescribed external electric field in one spatial dimension Hebenstreit *et al.* (2013a, b); Kasper *et al.* (2014). However, the role of background plasmas during pair production has not been investigated. By incorporating a nonperturbative amount of background particle fields, we turn real-time lattice simulations into numerical tools useful for plasma physics.

In this paper, we demonstrate plasma effects during pair productions using the scalar-QED model, where the electromagnetic (EM) fields evolve in a self-consistent manner, instead of being imposed from boundary and initial conditions. Using the scalar-QED model, we avoid the fermion doubling problem when discretizing the Dirac field Chandrasekharan and Wiese (2004) and focus on unambiguous plasma contributions. The scalar-QED model governs interactions between EM fields and spin-0 charged bosons, such as charged pions or Cooper pairs. Although laboratory plasmas are typically made of spin-1/2 charge fermions, classical plasma physics takes no account of particle spin-statistics at all. Therefore, to demonstrate that lattice simulations are useful for plasma physics, it is sufficient to study scalar QED, which has been used to describe laser-plasma interaction in relativistic-quantum regime Eliasson and Shukla (2011); Shi *et al.* (2016). In the classical-statistics regime, scalar QED is governed by Klein-Gordon-Maxwell (KGM) equations

(1) | |||||

(2) |

where we have used the natural units . In the above equations, is the complex scalar field, describing spin-0 bosons of species , whose charge is and mass is . The real-valued 1-form is the gauge field, describing spin-1 bosons, and is the field strength tensor. Charged bosons couple to the gauge field through the covariant derivative , and the gauge field couples with charged fields through the gauge-invariant current density

(3) |

where denotes complex conjugation of . By the famous Klein paradox Klein (1929), the charged scalar field cannot be interpreted as the probability amplitude of a single particle. A more appropriate interpretation is that the classical field is intrinsically a many-particle field, which can be represented as , where is the symmetrized many-body wave function, and the integration is carried out on the many-body configurations space Shi *et al.* (2016). Regardless of the interpretation, we can solve the KGM equations as coupled partial differential equations, whose solutions model the tree-level behavior of charged bosons interacting with EM fields.

This paper is organized as follows. In Sec. II, we develop a variational algorithm for solving the KGM equations. In Sec. III, we apply this algorithm to two example problems in plasma physics. The first example is the propagation of linear waves, where we compare numerical spectra with analytical dispersion relations. The second example is wakefield acceleration and pair production, when intense lasers interact with a 1D plasma slab. Conclusion and discussion are given in Sec. IV. In Appendix A, we discuss local conservation laws underlying our algorithm. In Appendix B, we summarize an explicit numerical scheme using the Lorenz gauge condition.

## Ii Variational algorithm

In the continuum, the KGM equations can be derived from the action , where the Lorentz-invariant and U(1)-gauge-invariant scalar-QED Lagrangian

(4) |

Here we have omitted the species subscript , and the summation of charged species is implied. By Noether’s theorem, the U(1)-gauge symmetry of the action

(5) |

implies charge conservation , where the current density is given by Eq. (3). Similarly, by the Lorentz symmetry, energy and momentum are also conserved , where the gauge invariant stress-energy tensor

(6) | |||||

Here, is the Minkowski metric with characteristics . This scalar-QED theory, omitting the self-coupling, is the underlying model of our algorithm on the discrete spacetime lattice. In fact, a variational algorithm for solving the KGM equations has already been developed in the numerical analysis community Christiansen and Halvorsen (2011), which shows superior charge conservation property when gauge symmetry is respected. In this paper, we rederive the variational algorithm in arbitrary gauge, using local energy conservation to justify the choice of Yee-type action Yee (1966) over Wilson-type action Wilson (1974), and emphasize on the application of such algorithm to plasma physics.

### A Discretization of fields and action

To solve the continuous system numerically, let us discretize the spacetime manifold using a rectangular lattice. Then the scalar field , namely, a function on the spacetime manifold, naturally lives on the vertexes of the discrete manifold

(7) |

where is the coordinate of the vertex. In comparison, the gauge 1-form naturally lives along the edges of the discrete spacetime manifold. For example, the - and -components

(8) | |||||

(9) |

where and . The minus sign comes from the Minkowski metric , which lower the index . In the above discretization, a half-integer index indicates which edge does the field resides along. For example, resides along the edge connecting vertices and , and is therefore the component of . Notice that since is a 1-form living along edges, only one of its four indices can take half-integer values, while the other three indices must take integer values. Moreover, to each edge of the lattice, the discrete 1-form only assigns the component of that is parallel to this edge (Fig. 1), while other components of are not assigned.

Now that we have discretized the fields, the gauge-covariant derivatives can be computed using the Wilson’s lines Wilson (1974). Since the gauge-covariant derivatives are 1-forms, they also lives along edges when discretized. For example, the - and -components of the pull-back gauge covariant derivatives are

(10) | |||||

(11) |

These pull-back covariant derivatives transform under the discrete U(1)-gauge symmetry [Eqs. (26)-(28)] as [Eq. (29)]. Analogously, one can define push-forward covariant derivatives, which we shall not use in this paper. The gauge field serves as the 1-form defining the connection on the U(1)-bundle, which enables parallel transport on the spacetime manifold.

To compute the field strength tensor , notice that is the curvature 2-form and hence lives on faces of the lattice upon discretization. For example, the time-like component is the electric field in the -direction, which can be computed by

(12) |

This component lives on the time-like face spanned by four vertices , , , and . Analogously, we can compute the space-like components of . For example, is the magnetic field in the -direction

(13) | |||||

This z-component of the magnetic field lives on the space-like face spanned by four vertices , , and . Notice that the sign of the discretized is determined by the orientation of the face. Since the above discretization respects geometric structures of exterior calculus, the Bianchi identities, namely, and the Faraday’s law, are automatically and exactly satisfied (Appendix A).

Using the discrete gauge-covariant derivatives and the discrete field strength, the action can be discretized by

(14) |

where and are the discrete fields. Here the subscript denotes vertexes, and denotes edges. In the discrete action, is the volume 4-form, and the summation runs over all cells of the lattice. In each unit cell, the discrete Lagrangian function

(15) |

where summations over unique , , and are implied. Here, the subscript denotes faces. Notice that in contrast to what is typically done in lattice gauge theories, here we have directly used the discrete field strength , instead of the Wilsonian plaquettes . We made this choice so that a discrete version of the local energy conservation law for EM fields is exactly satisfied when the coupling goes to zero (Appendix A) .

### B Equations of motion for discrete fields

Having discretized the action, the classical equation of motion (EOM) for the discrete field can be obtained by extremizing . Taking variation with and set , a discrete version of equation of Eq. (1) is

(16) | |||||

where the time index is explicit, the vertex-centered spatial index is abbreviated as , and summations in directions are implied. By taking variation with , we can obtain the EOM for , which is the complex conjugation of the above equation. The finite difference equation (16) is centered around vertexes, and couples with its eight nearest neighbors though , as illustrated by Fig. 2(a) in the -submanifold.

To find the equation for the electric field, take variation of with respect to the time-like component . Setting , we obtain a discrete version of the Gauss’ law, centered along time-like edges

(17) |

The charge density 1-form is the hodge dual of the charge density 3-form , discretized by

(18) |

When there are multiple charged species, the right-hand side (RHS) should sum over charge densities of all species. In Fig. 2(b), we illustrate the coupling pattern of the above finite difference equation.

To find the equations involving components of the magnetic field, take variation of with respect to the space-like components . For example, by setting , we obtain an equation advancing the x-component of the electric field in time by

(19) |

Here, is the abbreviated index for the body center, is the Levi-Civita symbol, and summations over repeated indexes are implied. The current density 1-form is the hodge dual of the current density 3-form . The hodge dual gives rise to a negative sign, so that the x-component of the current density is discretized by

(20) |

The finite difference equation (19) is the discrete version of the Maxwell-Ampère’s law centered around space-like edges, whose coupling pattern is illustrated in Fig. 2(c). When computing on the RHS, summation over charged species is implied.

In order to advance the above finite difference equations in time, we need to fix a gauge to eliminate the extra degree of freedom. Since the discrete action is U(1)-gauge invariant (Appendix A), we can choose any gauge. For example, one convenient choice is the Lorenz gauge . When discretized, the Lorenz gauge condition allows time advance in a very simple way [Eq. (43), Fig. 2(d)]. Another convenient choice is the temporal gauge . When discretized, remains zero on all time-like edges.

Having obtained discrete equations and fixed the gauge, an explicit time advance scheme can be constructed (Appendix B). We first initialize the simulation by giving values of at both and for every spatial lattice points . This is necessary because the KG equation [Eq. (1)] is a second order partial differential equation and therefore needs two initial conditions. Similarly, we need to give initial values of at and . Second, we use the discrete Gauss’ Law to calculate at by solving a system of linear equations [Eq. (40), Fig. 2b]. Notice that the continuous version of this equation can be written as , where the RHS is known. Because the unknowns on the left-hand side (LHS) involve only first order spatial derivative, the discrete Gauss’ law couples less number of points than the discrete Poisson’s equation, and is therefore easier to solve. Third, we use the Lorenz gauge condition to advance [Eq. (41), Fig. 2d]. Fourth, we use the discrete KG equation to calculate in terms of and [Eq. (42), Fig. 2a]. This involves exponentiation of and , whose values are already known at this point. Simultaneously, we can compute in terms of at previous time steps [Eq. (43), Fig. 2c], using the discrete Maxwell-Ampère’s law. Notice that the discrete Gauss’ Law is preserved during time advance, which is a consequence of the discrete local charge conservation law [Eq. (A)]. Finally, having computed both and at , we can move forward in the time loop by updating , with proper boundary conditions supplied. In similar fashion, explicit time advance schemes can be constructed when other gauge conditions are used.

## Iii Numerical examples

In this section, we demonstrate our numerical scheme using two examples. The first example is the propagation of linear waves, and the second example is laser-plasma interaction in one spatial dimensional (1D).

### A Linear wave spectra

To test our code implementation, we compare numerical spectra with analytical linear wave dispersion relations Hines and Frabkel (1978); Kowalenko *et al.* (1985); Eliasson and Shukla (2011); Shi *et al.* (2016). For small-amplitude waves, the dispersion relation constrains the wave frequency as function of the wave vector . In unmagnetized cold scalar-QED plasma, the tree-level dispersion relation of the transverse EM wave is

(21) |

where is the total plasma frequency, and is the plasma frequency of individual charged species . The other eigenmode is the longitudinal electrostatic wave, whose dispersion relation is

(22) |

where the cold plasma susceptibility

(23) |

The dispersion relation of electrostatic wave contains three branches. The gapless branch is the acoustic wave, the low frequency branch is the Langmuir mode, and the high frequency branch is the pair mode. While acoustic mode and Langmuir mode exist in classical plasmas, the pair mode only exists in relativistic-quantum plasmas Fuda and Furlani (1982). The pair mode can be excited when gamma photons inelastically scatter in high density plasmas, creating longitudinal oscillations in which virtual pairs are created and annihilated to carry the wave quanta.

We compute the numerical spectra in a single species plasma, in which immobile ions serve as homogeneous neutralizing background. To initialize the simulation so that a broad spectrum of linear waves are excited, the initial values of are given using small amplitude white noise with mean and standard deviation . Assuming the charged field is initially free, then its initial conditions can be given using the free field expansion , where is Minkowski inner product, and is the relativistic energy corresponding to momentum . From the above expansion, the momentum space distribution functions for particles and antiparticles are and the , respectively. Consider the simple example where the plasma is initially homogeneous and constituted of cold particles, namely, and , where is the background plasma density. Then, the free charged field , where is some random phase. When discretized, this free field corresponds to the initial conditions and . An ensemble of statistically equivalent initial conditions can then be constructed by randomly sample the phase and the gauge field.

After advancing the initial conditions in time using periodic boundary conditions, numerical spectra can be read out from simulations by taking discrete Fourier transforms of electric field components. Since the unmagnetized plasma is isotropic, it is sufficient to read out the dispersion relation in the tx-submanifold. In this submanifold, the spectra of either or correspond to the dispersion relation of transverse EM modes, and the spectrum of corresponds to the dispersion relation of longitudinal electrostatic modes. The ensemble-averaged power spectrum of [Fig. 4(a)] is indistinguishable from that of , and is well-traced by the analytical dispersion relation (black line) of the transverse EM wave [Eq. (21)], until where the spatial resolution is no longer sufficient. Similarly, the ensemble-averaged power spectrum of [Fig. 4(b)] is localized near three bands, corresponding to the cold acoustic mode, the Langmuir mode and the pair mode [Eq. (22)]. That the analytical dispersion relations are recovered by numerical spectra indicates that our solutions faithfully capture the propagation of linear waves up to the grid resolution.

### B Laser-plasma interaction

Having verified our code implementation, let us study laser-plasma interaction as another example, which can no longer be described self-consistently under the classical framework once the laser wavelength becomes too short and the field strength becomes too large. Before discussing our simulations in this relativistic-quantum regime, it is helpful to recall what happens in the classical regime Kruer (1988). Classically, when the plasma slab is under-dense, namely when the laser frequency , much of the laser will travel through the plasma slab, with some reflection and inverse Bremsstrahlung absorption. In an initially quiescent slab, the laser will propagate uneventfully, if its frequency stays away from the two-plasmon-decay resonance, and its intensity is not strong enough to grow instabilities within the pulse duration. Beyond nonlinear wave instabilities, when the laser field becomes relativistically strong, namely when , the ponderomotive force of a short laser pulse can expels a significant fraction of plasma electrons and form wakefield Pukhov and Meyer-ter Vehn (2002). The wakefield can then accelerate particles, generating energetic beams of particles and radiations trailing the laser pulse. When the beams are energetic enough, they may produce gamma photons through synchrotron radiation or Bremsstrahlung. The virtual gamma photons may then decay into electron-positron pairs through the trident process Bjorken and Chen (1967). Alternatively, the on-shell gamma photons may interact with ion potentials and produce pairs through the Breit-Wheeler process Breit and Wheeler (1934). Finally, when the laser field becomes even stronger, namely when , pairs may also be produced directly through the Schwinger process Schwinger (1951).

Many aspects of laser-plasma interaction can be studied using our new numerical tool. Here, to validate that our code can capture genuine relativistic-quantum effects, we select parameters in our 1D simulations to demonstrate transition from wakefield acceleration to Schwinger pair production as we increase the laser intensity. Notice that in 1D, the phase space is highly constrained. Using periodic boundary conditions in directions transverse to laser propagation, Schwinger pair production by laser fields is suppressed. This is because when transverse fields try to pull pairs apart, their wave functions are enforced to be the same by the periodic boundary condition, which prevents pairs from emerging out of vacuum fluctuations. Therefore, in 1D simulations, Schwinger pair production requires longitudinal field . To generate beyond the Schwinger field through wakefield, the plasma density must be extremely high. Heuristically, to produce on-shell pairs, the critical electric field needs to separate the pair by Compton wavelength within the Compton time , namely, . In the wave-breaking regime, , so the inequality requires that the plasma density be high enough such that the plasma frequency . In reality, at those densities, we should treat the electron Fermi degeneracy to capture the full physical effects. However, simulating instead a high-density bosonic plasma is just a toy model that tests our code, with the density picked so high that we can already see laser Schwinger pair production in 1D simulations.

With this basic understanding of how laser pair production happen in 1D, we choose parameters to suppress the trident and Breit-Wheeler processes, by treating ions as immobile homogeneous neutralizing background, so that there is no spiky ion potentials from which energetic “electrons” and gamma photons can scatter. The smooth ion background provides an electrostatic potential that initially confines the “electrons”. We initialize the charged boson wave function according to , where is the background ion density with a plateau of width and Gaussian off-ramps with . For density of the bosonic plasma to be high enough to enable pair production, we take so that the plasma frequency is enormous. The above wave function is a linear superposition of many eigenstates of the system. In our simulations, we let the wave function evolve to statistically stationary states through phase mixing, before we start to draw samples at random time intervals. The sampled wave functions are then used as initial conditions for , which are combined with initial values of a Gaussian pulse to construct an ensemble. The linearly-polarized Gaussian pulse is initialized in the vacuum region with zero carrier phase , where and . For the laser to be able to transmit the high-density plasma slab, we use a gamma-ray laser with frequency , for which semiclassical treatments are far from valid. The laser envelope is slowly varying (), and has full width at half maximum about twice the plasma skin depth. When intense laser pulse propagates, it can excite plasma waves, from which the laser can be Raman scattered.

With the above setup, the laser pulse simply travels through the plasma with some refraction and reflections when the laser field is weak (). More interesting phenomena happen when the laser field becomes strong. For example, when is relativistically strong but the resulting is below the Schwinger field, our simulation recovers what happens in classical plasmas McKinstrie and Startsev (1996); Naumova *et al.* (2004); Geyko *et al.* (2009). First, let us look at what happens to charged particles. After the laser enters the plasma, beams of “electrons” are formed in the forward direction by both ponderomotive snow-plow and laser wakefield acceleration. At the same time, some “electrons” are splashed in the backward direction from strongly-driven plasma boundaries (Fig. 5a, c).
Next, for the laser pulse, its center (solid black lines) and half widths (dotted black lines) are well-traced by geometric optics in the -space (Fig. 6a), as well as in the -space (Fig. 6c, dashed white line), because the background plasma is smooth on the laser wavelength scale. Beyond geometric optics, as the laser travels through the plasma slab, ponderomotive expulsion of “electrons” cause the laser pulse to adiabatically loose a small amount of energy in the form of frequency redshift (Figs. 6a,c and 7b). In addition, the laser excites plasma waves, from which the laser is Raman-scattered in both forward and backward directions. In the insert of Fig. 6c, the final spectrum (red) shows distinctive Raman scattering peaks at up to , and second harmonics peaks and in the forward direction. In the backward direction, peaks at and can also be identified unambiguously.

When we increase the laser field beyond the Schwinger threshold (). For example, when (), a large amount of pairs are produced (Figs. 5b, d). A very small fraction of pairs are produced and trapped in the laser wakefield, forming low-luminosity “electron” (negative charge density, blue) and “positron” (positive charge density, red) beams that leave the plasma slab from its right boundary. In contrast, a much larger fraction of pairs are produced when the backscattered EM wave, whose intensity is near the Schwinger threshold (Fig. 6b), interacts with forward-propagating plasma waves. “Positrons” produced in this way form high-luminosity collimated beams, leaving the plasma slab from its left boundary. Apart from these beams, many “positrons” never manage to leave the plasma slab. These trapped “positrons” have large probabilities to annihilate with “electrons” in the highly constrained 1D phase space.
Due to pair creation and particle acceleration, the laser initially looses a significant amount of energy, until pair creation and annihilation roughly balance (Figs. 6b,c and 7b). At that point, the -spectrum of the laser is substantially broadened (Fig. 6d). Such a spectral broadening is expected from general wave action considerations Wilks *et al.* (1988); Dodin and Fisch (2010), which predict frequency up-shift due to pair creation, and frequency down-shift due to pair annihilation and plasma expulsion. In the insert of Fig. 6d, the final EM wave spectrum (red) shows distinctive annihilation bumps near integer multiples of “electron” rest mass. These annihilation peaks are very broad since “electrons” and “positrons” annihilate with large kinetic energy. Finally, notice that no pair is produced when the laser travels through the vacuum region, which is expected in 1D. It is remarkable that very rich physics can already be captured by simply solving the classical field equations with proper initial conditions.

To extract observables from simulations, the charge density (Figs. 5a, b) is computed using Eq. (18), which includes no contribution from background ions. Therefore, negative charge (blue) indicates “electron” density in excess of “positron” density, whereas positive charge (red) indicates the contrary. The energy density of the charged field (Figs. 5c, d) and the EM fields (Figs. 6a, b) are computed using Eqs. (35) and (36), respectively. To compute -spectra of EM waves (Figs. 6c, d), notice that a monochromatic EM wave satisfies . Upon discretization, this relation remains exactly satisfied if we take and , where is the positive solution of the local numerical dispersion relation . In the discrete version of , it is necessary that we take , and center on time-like faces . A similar relation holds for the and components, which are subdominant in our simulations. Using these momentum-space Faraday’s law, the -spectrum of right-propagating EM waves () and left-propagating EM waves () can be separated from the spatial Fourier transforms of electric and magnetic fields.

Results presented in Figs. 5-7 are averaged over an ensemble of 200 simulations with statistically equivalent initial conditions. The ensemble average starts to show convergence for tens of realizations. In these simulations, temporal gauge is used, and periodic boundary conditions are employed for both and . We choose resolutions and , high enough that the fastest dynamics is resolved and the simulation results converge. The 1D box is large enough such that the laser does not transit the spatial domain before we terminate the simulations.

## Iv Discussion and Summary

In this paper, we develop an algorithm for solving the Klein-Gordon-Maxwell’s equations [Eqs. (1) and (2)], which can be used to model bosonic plasmas in the relativistic-quantum regime. This algorithm is derived by first discretizing the action [Eq. (14)] in a way that respects the U(1)-gauge symmetry. We then take variations with respect to the discrete fields to find their classical equations of motion. The resultant variational algorithm guarantees that the Bianchi identities, namely, and the Faraday’s law, are automatically and exactly satisfied. The remaining equations of motions are the discrete Gauss’s law [Eq. (17)], which can be used to initialize the simulation; the discrete Klein-Gordon’s equation [Eq. (16)], which can be used to advance the charged field; and the discrete Maxwell-Ampère’s law [Eq. (19)], which can be used to advance the gauge field. After fixing a gauge, explicit scheme for advancing the discrete fields in time can be constructed (Appendix B, Fig. 3). Our variational scheme respects local conservation laws (Appendix A), and can be easily parallelized using domain decomposition. Moreover, the numerical scheme we have developed can be inherently mimicked by quantum systems with local couplings Wiese (2013); Martinez *et al.* (2016), which can be efficiently realized using quantum parallelism Feynman (1986); Lloyd (1996).

The numerical scheme we have developed has a number of advantages over conventional methods for simulating plasmas. As comparison, the two conventional methods that can fully simulate kinetic effects are the particle-in-cell (PIC) solvers and the Vlasov solvers.
The PIC solvers represent particles in the continuum, while representing EM fields on grids. Particles feel EM fields through interpolations, and EM fields feel particles through depositions. Using proper smoothing functions, these two steps can preserve gauge symmetry and symplectic structures, thereby respect local conservation properties when used in geometrical algorithms Squire *et al.* (2012); Xiao *et al.* (2013, 2015); Qin *et al.* (2015). Nevertheless, interpolation and deposition introduce artificial collisions, which are absent in physical systems.
In the other scheme, the Vlasov solvers, EM fields are represented on the three-dimensional space, while particles are represented in the six-dimensional phase space. Particles are directly forced by fields on spatial grids, while the fields feel particles though velocity space integrals, which requires resolving three extra dimensions with substantial computational cost.
In contrast, our algorithm represents both particles and EM fields on the same grid. Therefore, there is no need for interpolations and depositions as in the case of PIC solvers, nor is there need for resolving extra velocity space dimensions as in the case of Vlasov solvers.
Our algorithm folds the phase space dynamics of charged particles into the complex plane, and enables modeling of relativistic and quantum dynamics in regimes that cannot be described using semiclassical treatments. In the example of linear waves (Sec. III.A), we show that relativistic-quantum wave dispersion relations can be recovered. Moreover, using the example of a gamma-ray laser interacting with a dense plasma slab (Sec. III.B), we show that our algorithm naturally allows pair production when the laser intensity exceeds the Schwinger threshold.

Of course, the advantages of our algorithm come at an expense. The expense comes from the necessity of resolving the relativistic-quantum scales, which are much smaller than scales that classical plasma physics typically deals with. The coarsest resolution needed in relativistic quantum plasma simulations is determined by the lowest energy scale of the problem, which is the rest mass of electrons MeV, corresponding to time scale of s, and spatial scale of m. This resolution requirement can be seen from the discrete KG equation [Eq. (42)], in which we must have in order for . Moreover, since we are solving a system of hyperbolic partial differential equations, the Courant–Friedrichs–Lewy (CFL) condition must be satisfied. Finally, it is worth noting that high resolution is required for large gauge fields. Since the gauge field appears through the Wilson’s lines in complex exponentials [Eqs. (10) and (11)], the discrete theory is invariant under the gauge transformation . Therefore, the discrete gauge field lives on the torus , which has a very different topology than . Therefore, the step size must be small enough such that , in order to avoid exciting topological modes that are absent in the continuous theory.

In summary, we develop a variational algorithm for solving the Klein-Gordon-Maxwell equations, which described tree-level scalar QED in the classical-statistical regime and may be solved using quantum computers. We demonstrate that remarkably rich physics are contained in solutions to classical field equations, which can be used to model high-density plasmas interacting with short-wavelength electromagnetic fields. Our work uses scalar QED as a toy model to make the case that real-time lattice QED is a powerful tool for plasma physics in the strong-field regime, where relativistic-quantum scales overlap with collective plasma scales. The applications of lattice spinor-QED to laboratory relevant plasma conditions remain to be demonstrated in the future.

###### Acknowledgements.

The authors are grateful to Qun Wang and Sebastian Meuren for helpful discussions. This rsearch is supported by NNSA Grant No. DE-NA0002948 and DOE Research Grant No. DEAC02-09CH11466. Jianyuan Xiao is supported by National Magnetic Confinement Fusion Energy Research Project (2015GB111003, 2014GB124005), National Natural Science Foundation of China (NSFC-11575185, 11575186, 11305171), JSPS-NRF-NSFC A3 Foresight Program (NSFC-11261140328), Chinese Scholar Council (201506340103), Key Research Program of Frontier Sciences CAS (QYZDB-SSW-SYS004), and the GeoAlgorithmic Plasma Simulator (GAPS) Project.## Appendix a Geometric Identities and Local Conservation Laws

When discretizing the gauge 1-form and calculating the curvature 2-form in Sec. II.A, geometric structures of discrete exterior calculus are respected. Consequently, the identity holds for the discrete exterior derivative. In components, this Bianchi identity can be written as . One nontrivial identity, corresponding to all indexes being spatial, is . When discretized, this identity becomes

(24) |

The other nontrivial identity, corresponding to two spatial indexes and one temporal index, is the Faraday’s law , whose discrete version is

(25) |

The above finite difference equations are automatically satisfied in our algorithm by geometric constructions, in contrast to standard elecromagnetic algorithms, such as the Yee’s algorithm Yee (1966), in which the Faraday’s law needs to be solved as a dynamical equation.

In addition to the above geometric identities, we also have a number of local conservation laws. The first is local charge conservation, which is a direct consequence of local U(1)-gauge symmetry. Under the continuous U(1)-gauge transformation

(26) | |||||

(27) | |||||

(28) |

where is any real-valued function living on vertexes. These transformations leave the discrete face-centered field strength tensor invariant, while transforming the pull-back covariant derivative by

(29) |

Since the action is invariant, we can use the classical field equations and write

(30) |

Substituting in the infinitesimal transformation , the above identity is equivalent to the discrete charge conservation law

Here, the charge density is given by Eq. (18), and the current density is given by Eq. (20), and the sign is due to the Minkowski metric. It is straightforward to check that the above discrete charge conservation law is compatible with the discrete Gauss’ law [Eq. (17)] and the discrete Maxwell-Ampère’s law [Eq. (19)]. In the numeric example discussed in Sec. III.B, the total charge is constant up to the machine precision (Fig. 7a), both when the laser field is below () and above () the Schwinger field.

Moreover, the discrete action is invariant under translations on the discrete spacetime manifold. Although the symmetry group in this case is discrete and hence the Noether’s theorem does not immediately apply, we do have local energy conservation laws for the charged field and EM fields separately when their coupling vanishes. Using the classical field equations and , as well as the Bianchi identity, we have the following identity

(31) | |||||

where the vertex-centered time covariant derivative

(32) |

After rearranging terms, the above identity gives rise to the local energy conservation law

(33) |

where the sign is again due to the Minkowski metric. The energy density can be separated into three terms

(34) |

where the energy density of the charged field is

(35) | |||||

and the energy density of the EM fields is

(36) |

The energy density correction can take many different forms, each has a corresponding error term at finite-resolution. As expected, the energy density is U(1)-gauge invariant, so is the momentum density, which can be split into two terms

(37) |

The momentum density of the charged field is

(38) |

and the momentum density of the EM fields is