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

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 one-dimensional 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.

Figure 1

Discretization of the $txy$ submanifold of space time using a rectangular lattice. The discrete function ${\varphi }_v$ lives on the vertexes (blue squares). For example, ${\varphi }_{i,j,k}^n=\varphi (t_n,x_i,y_j,z_k)$ lives on the vertex $(n,i,j,k)$. The discrete 1-form $A_e$ lives along edges (red circles). For example, the $t$ component $A_{i+1,j,k}^{n+1/2}=A^0(t_n+\mathrm{\Delta }t/2,x_{i+1},y_j,z_k)$ lives along the timelike edge connecting vertexes $(n,i+1,j,k)$ and $(n+1,i+1,j,k)$, and the $x$ component $A_{i+1/2,j,k}^n=-A^1(t_n,x_i+\mathrm{\Delta }x/2,y_j,z_k)$ lives along the spacelike edge connecting vertexes $(n,i,j,k)$ and $(n,i+1,j,k)$. The discrete 2-form $F_f$ lives on faces (green crosses). For example, electric field $E_{i+1/2,j,k}^{n+1/2}=E_x(t_n+\mathrm{\Delta }t/2,x_i+\mathrm{\Delta }x/2,y_j,z_k)$ lives on the timelike face spanned by vertexes $(n,i,j,k),(n+1,i,j,k),(n+1,i+1,j,k)$, and $(n,i+1,j,k)$; magnetic field $B_{i+1/2,j+1/2,k}^{n+1}=B_z(t_{n+1},x_i+\mathrm{\Delta }x/2,y_j+\mathrm{\Delta }y/2,z_k)$ lives on the spacelike face spanned by vertexes $(n+1,i,j,k),(n+1,i,j+1,k),(n+1,i+1,j+1,k)$, and $(n+1,i+1,j,k)$.

Figure 2

Coupling pattern of ${\varphi }_v$ (blue squares), $A_e$ (red circles), and $F_f$ (green crosses) in the $tx$ submanifold. (a) The discretized KG equation [Eq. (16)] couples ${\varphi }_v$ with its nearest neighbors though $A_e$. (b) The discretized Gauss's law [Eq. (17)] couples $E_{i-1/2}$ and $E_{i+1/2}$ through ${\varphi }_v$, centered around the common timelike edge. (c) The Maxwell-Ampère's law [Eq. (19)] couples $E^{n+1/2}$ to $E^{n-1/2}$ through ${\varphi }_v$ and $B^n$ (not depicted here), centered around the common spacelike edge. (d) The Lorenz gauge condition couples values of $A_e$ that share the same vertex.

Figure 3

Time evolution scheme for discrete KGM equations using the Lorenz gauge. As initial conditions, the values of ${\varphi }_v(n=0)$ and ${\varphi }_v(n=1)$ are given (blue squares), and so are $A_e(n=0)$ and $A_e(n=1/2)$ (red circles). Then the Gauss's law [Eq. (17)] is used to calculate $A_e(n=1)$. On entering the time loop, the first step is to calculate $A^{n+1/2}$ using the Lorenz gauge condition. The second step is to use the KG equation [Eq. (16)] to calculate ${\varphi }^{n+1}$ and, concurrently, use the Maxwell-Ampère's law [Eq. (19)] to calculate $A^{n+1}$. The time loop is advanced by $n\to n+1$ and then repeat.

Figure 4

Power spectra (color) of the transverse electric field $E_y$ (a) and the longitudinal electric field $E_x$ (b) are well traced by tree-level dispersion relations (black lines) up to the grid resolution. The power spectra are averaged over an ensemble of 100 simulations with statistically equivalent initial conditions. In these simulations, immobile ion background is homogeneous. The charge $q=0.3$, such that the fine structure constant $q^2/4\pi \approx 1/137$ is physical. The unperturbed background plasma density is extremely high, such that the plasma frequency ${\omega }_p=0.85m$ can be shown on the same scale as $m$. The resolution $m\mathrm{\Delta }x=0.04$ and $m\mathrm{\Delta }t=0.02$. The number of spatial grid point is $L=512$, and the total number of time steps, including the initial conditions, is $T=1024$. The dashed gray lines are the light cone.

Figure 5

Charge density (a), (b) and energy density (c), (d) of the $\varphi$ field. When the gamma-ray laser (${\omega }_0=0.7m$) is relativistic $(a\approx 1)$, but not strong enough to produce Schwinger pairs $(E_x\approx 0.3E_c)$, “electrons” are expelled by the laser ponderomotive force, accelerated by the wakefield, and splashed from the plasma boundaries (a), (c). On the other hand, when the laser field exceeds the Schwinger threshold $(a\approx 16,E_x\approx 5E_c)$, copious pairs are produced when the laser interacts with plasma waves (b), (d). The spin-0 “electrons” are initially confined by a smooth immobile neutralizing background, with a density plateau $n_0=m^3$ and a Gaussian off-ramp $\sigma =20/m$. The trajectories of the pulse center (black lines) and the pulse half widths (dashed lines) are well traced by geometric optics. Both the charge density (normalized by $e{m}^3$) and the energy density (normalized by $m^4$) are averaged over an ensemble of size 200. The resolutions are such that $m\mathrm{\Delta }x=0.04$ and $m\mathrm{\Delta }t=0.005$.

Figure 6

Total energy density of EM fields (a), (b), and the power spectral density of its transverse components (c), (d). The inserts show the initial (blue) and final (red) spectra of EM waves. When $a\approx 1$ ($E_x\approx 0.3E_c$) is below the Schwinger field, the laser excites plasma waves and is Raman scattered (a), (c). The time evolution of the main pulse is well traced by geometric optics (dashed lines). On the other hand, when the laser field $a\approx 16$ ($E_x\approx 5E_c$) is above the Schwinger field, a noticeable amount of energy is lost due to pair production (b), and the $k$ spectrum is substantially broadened (d). The field energy density is normalized by the Schwinger field $E_c^2$ and is averaged over an ensemble of size 200. The resolutions are such that $m\mathrm{\Delta }x=0.04$ and $m\mathrm{\Delta }t=0.005$. The dotted gray lines mark where the geometric-optics trajectory of the pulse center crosses the plasma plateau boundaries.

Figure 7

Evolution of total charge (a) and total energy (b) in the numerical example discussed in Sec. 3b, where periodic boundary conditions are used. The total charge remains constant up to machine precision, both when $E<E_c$ (cyan) and $E>E_c$ (blue). When $E<E_c$ is below the Schwinger field, a small amount of energy is transfered from the electromagnetic field (magenta) to the charged field (cyan) due to wakefield acceleration and plasma wave excitation, while the total energy (gray) remains constant. In contrast, when $E>E_c$, a large amount of laser energy (red) is consumed by pair production. The energy of the charged field (blue) significantly increases, while the total energy (black) remains constant. The total charge $Q^{n+1/2}={\sum }_s{J}_s^{n+1/2}$ is normalized by the total ion charge, and the total energy ${\mathcal{U}}^{n+1/2}={\sum }_s{\mathcal{H}}_s^{n+1/2}$ is normalized by $m^3/\mathrm{\Delta }x$. The vertical dotted gray lines mark the time when the laser pulse center enters and leaves the plasma plateau boundaries.