Strong-field quantum electrodynamics (SF-QED) plays a crucial role in ultraintense laser–matter interactions and demands sophisticated techniques to understand the related physics with new degrees of freedom, including spin angular momentum. To investigate the impact of SF-QED processes, we have introduced spin/polarization-resolved nonlinear Compton scattering, nonlinear Breit–Wheeler, and vacuum birefringence processes into our particle-in-cell (PIC) code. In this article, we provide details of the implementation of these SF-QED modules and share known results that demonstrate exact agreement with existing single-particle codes. By coupling normal PIC simulations with spin/polarization-resolved SF-QED processes, we create a new theoretical platform to study strong-field physics in currently running or planned petawatt or multi-petawatt laser facilities.
Apart from these kinetic effects, spin/polarization effects also arise with the possibility of generating polarized high-energy particle beams or when particles traverse large-scale intense transient fields in laser–plasma interactions. Classically, the spin of a charged particle will precess around the instantaneous magnetic field, i.e., ds/dt ∝ B × s, where s denotes the classical spin vector.7 In storage rings, owing to radiation reaction, the spin of an electron/positron will flip to the direction parallel/antiparallel to the external magnetic field in what is known as the Sokolov–Ternov effect8 (an unpolarized electron beam will be polarized to a degree of ∼92.5%), and a similar process also occurs in NCS.9–11 Some recent studies have shown that with specific configurations, for example, when elliptically or linearly polarized lasers scatter with high-energy electron bunches (or plasmas), the polarization degree of the electrons can reach 90% and be used to diagnose transient fields in plasmas.12,13 Meanwhile, the photons created by NCS can be polarized, and when these polarized photons decay into electron/positron pairs, the contribution to the probability from polarization can reach ∼30%,14 and will be inherited by the subsequent QED cascade. For example, in laser–plasma/beam interactions, the polarization degree for linearly polarized (LP) photons is about 60% or higher, and for circular polarized (CP) γ-photons, it can reach 59% when longitudinally polarized primaries are employed.11,15–17
Analytical solutions in the case of ultraintense laser-matter interactions are scarce owing to the high nonlinearity and complexity of the problem. Moreover, the microlevel processes such as ionization, recombination, and Coulomb collisions, coupled with the complicated configurations of lasers and plasmas, make explicit derivations almost impossible. Fortunately, computer simulation methods provide alternative and more robust tools to study those unsolvable processes, even in more realistic situations.18 In general, simulation methods for laser–plasma (ionized matter) interactions can be categorized as kinetic or fluid simulations: specifically, kinetic methods include the Fokker–Planck (F–P) equation (or the Vlasov equation for the collisionless case) and the particle-in-cell (PIC) method, while fluid methods mainly use the magnetohydrodynamic (MHD) equations.19 Among these methods, both F–P and MHD discretize the momentum space of particles and are prone to the nonphysical multistream instability, which may obscure the real physics, such as the emergence of turbulence, physical instabilities, etc. In comparison, the PIC method can provide much more detailed information on the discrete nature and intrinsic statistical fluctuations of the system, regardless of the stiffness of the problem. Therefore, the PIC method has been widely used in the simulation of ultraintense laser–plasma interactions.18–20
Thanks to emerging PIC simulation methods, the development of parallelism, and large-scale cluster deployment, simulations of laser–plasma wakefield acceleration, laser ion acceleration, THz radiation, as well as SF-QED, have become accessible for general laser–plasma scientists.18,21–24 However, the spin and polarization properties of the plasma particles and QED products have not been widely considered in mainstream studies, owing to a lack of appropriate algorithms. In some recent studies, spin- and-polarization resolved SF-QED processes have been investigated in laser–beam colliding configurations, and it has been shown that these processes are prominent in generating polarized beams.10,11,14,16,17,25 Locally constant approximations of the relevant probabilities can be readily introduced into any PIC code.
In this paper, we briefly review the common PIC simulation algorithms and present some recent implementations in spin/polarization averaged/summed QED. The formulas and algorithms for spin/polarization-dependent SF-QED processes are given in detail and have been incorporated into our PIC code SLIPs (“spin-resolved laser interaction with plasma simulation code”). The formulas and algorithms presented in this paper, especially the polarized version, can be easily adopted by any other PIC code and used to simulate the ultraintense laser–matter interactions that are already relevant or will become so in near-future multi-petawatt (PW) to exawatt (EW) laser facilities,26 such as Apollo,27,28 ELI,29 SULF,30 and SEL. Throughout the paper, Gaussian units will be adopted, and all quantities are normalized as follows: time t with 1/ω (i.e., t′ ≡ t/(1/ω) = ωt), position x with 1/k = λ/2π, momentum p with mec, velocity v with c, energy ɛ with mec2, EM fields E and B with mecω/e, force F with mecω, charge q with e, charge density ρ with k3e, and current density J with k3ec, where λ and ω = 2πc/λ are the reference wavelength and frequency, respectively.
II. PIC ALGORITHM
Simulation of laser–plasma interactions involves two essential components: the evolution of the EM field and the motion of particles. The corresponding governing equations are the Maxwell equations (with either A–ϕ or E–B formulations) and the Newton–Lorentz equations. Therefore, the fundamentals of PIC codes consist of four kernel parts: force depositing to particles, particle pushing, particles depositing to charge and current densities, and solving the Maxwell equations; see Fig. 1. Here, we review each part briefly (these algorithms are used in SLIPs) and refer to the standard literature or textbooks for more details.18,19
A. Particle pushing
B. Field solving
The standard finite-difference method in the time domain for the Maxwell equations is to discretize field variables on a spatial grid and advance forward in time. Here, following the well-known Yee-grid approach,34 we put E and B as in Fig. 2(a), which automatically satisfies the two curl equations. For lower-dimensional simulations, extra dimensions are squeezed, as shown in the 2D example in Fig. 2(b). In these dimension-reduced simulations, field components in the disappeared dimensions can be seen as uniform, i.e., the gradient is 0.
Using Esirkepov’s method of current deposition,35 the current is calculated from the charge density via charge conservation, i.e., ∂tρ + ∇ · J = 0. Once the initial condition obeys Gauss’s law, ∇ · E = ρ, this law is automatically embedded. This can be verified by taking the gradient of Eq. (14): 0 = ∇ · (∇ × B) = ∂t(∇ · E) + ∇ · J = ∂t(∇ · E − ρ), i.e., the temporal variation in the violation of Gauss’s law is 0. Therefore, in the field solver, only the two curl equations are solved. We take the Ey and Bz components as examples here:
Here, the lower indices with i, j, k denote the spatial discretization and upper indices with n indicate the time discretization. The time indices are assigned using the leapfrog algorithm; see Sec. II F.
C. Current deposition
D. Force deposition
We deposit the updated field variables from the Maxwell solver to the particles for calculating acceleration or further SF-QED processes. The field deposition to the particles follows a similar procedure as the charge density deposition. For each particle at position xr, we find its nearest grid point and its nearest half grid point (i, j, k)h = floor(xr/Δx), where Δx = (Δx, Δy, Δz) is the spatial grid size. We then weight the field to the particle by summing over all nontrivial terms of W(i, j, k) · F(i, j, k), where W(i, j, k) is the particle weighting function (see Sec. II E for more details) on the grid (half grid) (i, j, k) and F(i, j, k) is the field component of E or B on the spatial grid with proper staggering according to Fig. 2.
E. Particle shape function
F. Time ordering
In SLIPs, the simplest forward method is used to discretize all differential equations that are reduced to first order with respect to time.18 To minimize the errors introduced by the discretization, some variables are updated at integer time steps and others at half-integer time steps. For example, the EM field variables E and B are updated alternately at integer and half-integer time steps, and the position x and momentum p of particles are updated alternately as well; see Fig. 3. The leapfrog updating is also applied to the current deposition and field interpolation.
III. QED ALGORITHM
This section presents some SF-QED processes (with unpolarized and polarized versions) that are relevant for laser–plasma interactions. The classical and quantum radiation corrections to the Newton–Lorentz equations, namely, the Landau–Lifshitz equation and the modified Landau–Lifshitz equation, and their discretized algorithms are reviewed first. The classical- and quantum-corrected equations of motion (EOM) for the spin, namely the Thomas–Bargmann–Michel–Telegdi equation and its radiative version, and their discretized algorithms, are reviewed next. NCS with unpolarized and polarized version and their Monte Carlo (MC) algorithms are reviewed. NBW pair production with unpolarized and polarized versions and their MC implementations are presented as well. Finally, the implementations of high-energy bremsstrahlung and vacuum birefringence under the conditions of weak pair production (χγ ≲ 0.1) are briefly discussed.
A. Radiative particle pusher
Charged particles moving in strong fields can emit either classical fields or quantum photons. This leads to energy/momentum loss and braking of the particles, i.e., radiation reaction. A well-known radiative EOM for charged particles is the Lorentz–Abraham–Dirac (LAD) equation.36 However, this equation suffers from the runaway problem, since the radiation reaction terms involve the derivative of the acceleration. To overcome this issue, several alternative formalisms have been proposed, among which the Landau–Lifshitz (LL) version is widely adopted.37 The LL equation can be obtained from the LAD equation by applying iterative and order-reduction procedures,38,39 which are valid when the radiation force is much smaller than the Lorentz force. More importantly, in the limit of ℏ → 0, the QED results in a planewave background field are consistent with both the LAD and LL equations.40,41 Depending on the value of the quantum nonlinear parameter χe (defined in Sec. III A 1), the particle dynamics can be governed by either the LL equation or its quantum-corrected version.1,23,37,42
1. Landau–Lifshitz (LL) equation
2. Modified Landau–Lifshitz (MLL) equation
3. Algorithms for the radiative pusher
B. Spin dynamics
The consideration of electron/positron spin becomes crucial in addition to the kinetics when plasma electrons are polarized or when there is an ultrastrong EM field interacting with electrons/positrons and γ-photons. The significance of this aspect has been highlighted in the recent literature, particularly in the context of relativistic charged particles in EM waves and laser–matter interactions.53,54 This issue can be addressed either by employing the computational Dirac solver55 or by utilizing the Foldy–Wouthuysen transformation and the quantum operator formalism, such as through the reduction of the Heisenberg equation to a classical precession equation.56,57 However, these approaches are not directly applicable to many-particle systems. Here and throughout this paper, the spin is defined as a unit vector S. In the absence of radiation, the electron/positron spin precesses around the magnetic field in the rest frame and can be described by the classical Thomas–Bargmann–Michel–Telegdi (T-BMT) equation. This equation is equivalent to the quantum-mechanical Heisenberg equation of motion for the spin operator or the polarization vector of the system.7,56,57 When radiation becomes significant, the electron/positron spin also undergoes flipping to quantized axes, typically aligned with the magnetic field in the rest frame. By neglecting stochasticity, this effect can be appropriately accounted for by incorporating the radiative correction to the T-BMT equation, which is analogous to the quantum correction to the LL equation.
1. Thomas–Bargmann–Michel–Telegdi (T-BMT) equation
2. Radiative T-BMT equation
3. Algorithms for simulating spin precession
Figure 6 presents a comparison between the T-BMT and radiative T-BMT equations for different cases: Lorentz equation + T-BMT equation (A), Lorentz equation + radiative T-BMT equation (B), LL equation + radiative T-BMT equation (C), and MLL equation + radiative T-BMT equation (D). The evolution of each spin component depends on different terms. In our setup, the magnetic field is along the z direction, and so the spin precession occurs in the x–y plane, affecting Sx and Sy. The radiation reaction mainly affects Sz. In the case without radiation reaction (case A), Sx and Sy oscillate owing to precession and are conserved in Fig. 6(d). In the case with only spin radiation reaction (case B), Sx is strongly damped by the term (dS/dt)R. Sy and Sz oscillate owing to the combined effects of precession and radiation reaction, as shown in Figs. 6(a) and 6(b). When both spin and momentum radiation reactions are included (case C), the particle momentum and energy decrease, i.e., γe decreases, which lowers the spin radiation reaction term (dS/dt)R(χe) and the damping of Sx and Sz [see Fig. 6(c) for a comparison of cases B, D, and C in terms of Sz amplitude]. Simultaneously, the precession term (dS/dt)T ∝ B/γe grows with decreasing γe, which amplifies the oscillation of Sy, as shown by the contrast between B (Lorentz), D (MLL), and C (LL) in Fig. 6(b).
C. Nonlinear Compton scattering (NCS)
When the radiation is strong (χe ≳ 0.1), its stochastic nature can no longer be neglected in the laser–beam/plasma interactions. Also, the photon dynamics should be taken into account. In this regime, the full stochastic quantum process is required to describe the strong radiation, i.e., nonlinear Compton scattering (NCS).2,60,61 Therefore, the radiation reaction and photon emission process will be calculated via MC simulation based on the NCS probabilities. The electron/positron spin and the polarization of the NCS photons will be also included in the MC simulations.
1. Spin-resolved/summed NCS
Finally, the polarization of the emitted photon is determined under the assumption that the average polarization is in a mixed state. The basis for the emitted photon is chosen as two orthogonal pure states with Stokes parameters (, , )/, where . The probabilities of photon emission in these states, , are given by Eq. (49). A stochastic procedure is defined using the fourth random number r4: if , the polarization state will be chosen; otherwise, the polarization state will be assigned as . Here, and .
Between photon emissions, the electron dynamics in the external laser field are described by the Lorentz equation dp/dt = −e(E + β × B) and are simulated using the Boris rotation method, as shown in Eqs. (5)–(10). Owing to the smallness of the emission angle for an ultrarelativistic electron, the photon is assumed to be emitted along the parental electron velocity, i.e., pf ≈ (1 − ωγ/|pi|)pi. Besides, in this simulation, interference effects between emissions in adjacent coherent lengths (lf ≃ λL/a0) are negligible when the employed laser intensity is ultrastrong, i.e., a0 ≫ 1. Therefore, the photon emissions occurring in each coherent length are independent of each other.
Examples of the electron dynamics and spin can be seen in Fig. 8: clearly, the average value matches the MLL equations for dynamics and the MLL + radiative T-BMT equations for spins. The beam evolution is also shown in Fig. 9. The energy spectra of electrons and photons, as well as the photon polarization, can be seen in Fig. 10.
2. Definition and transformation of Stokes parameters
D. Nonlinear Breit–Wheeler (NBW) pair production
When the energy of a photon exceeds the rest mass of an electron–positron pair, i.e., ωγ ≥ 2mec2, and the photon is subjected to an ultraintense field a0 ≫ 1, the related nonlinear quantum parameter χγ can reach unity. Here, and is approximately equal to in the colliding geometry. In this scenario, the photon can decay into an electron–positron pair through the nonlinear Breit–Wheeler pair production (NBW) process (ωγ + nωL → e+ + e−).2 In Refs. 25, 64, and 70, 71 a spin- and polarization-resolved NBW MC method was proposed, and here we follow the methods described in detail in Ref. 72.
1. NBW probability
2. MC algorithm
The algorithm for simulating pair creation with polarization is illustrated in Fig. 11. At every simulation step Δt, a pair is generated with positron energy ɛ+ = r1ɛγ when the probability density P ≡ d2Wpair/dɛ+dt · Δt of pair production is greater than or equal to a random number r2 within the range [0, 1). Here, d2Wpair/dɛ+dt is computed using Eq. (75). The momentum of the created positron (electron) is parallel to that of the parent photon, and the energy of the electron ɛ− is determined as ɛγ − ɛ+. The final spin states of the electron and positron are determined by the four probability densities P1,2,3,4, each representing spin parallel or antiparallel to the SQA, where P1,2,3,4 is computed from Eq. (64). Finally, a random number r3 is used to sample the final spin states for the electron and positron. Note that here all random numbers are sampled uniformly from [0, 1), as in the NCS algorithm. An example of the production of secondary electrons and positrons resulting from a collision between a laser and an electron beam is illustrated in Fig. 12.
E. High-energy bremsstrahlung
The bremsstrahlung implementation is based on direct MC sampling. Given an incident electron with energy Ee and velocity v, the probability of triggering a bremsstrahlung event is calculated as Pbr = 1 − eΔs/λ, where Δs = vΔt, v = |v| is the incident particle velocity, Δt is the time interval, λ = 1/nσ(Ee), n is the target particle density, and σ(Ee) is the total cross-section. A random number r1 is then generated and compared with Pbr. If r1 < Pbr, then a bremsstrahlung event is triggered. The energy of the resulting photon is determined by generating another random number r2, which is then multiplied by σbr(Ee) to obtain the energy ratio y through σ(y, Ee) = σ(Ee)r2. Finally, a photon with energy ℏω = Eey and momentum direction k/|k| = v/|v| is generated. To improve computational efficiency, low-energy photons are discarded by setting a minimum energy threshold. This probabilistic approach is similar to the method used to calculate the random free path.76 The implementation of Bethe–Heitler pair production follows a similar process.
The implementation of bremsstrahlung emission was tested using Geant4 software,74 which is widely used for modeling high-energy particle scattering with detectors. In this study, we utilized electron bunches of 100 MeV and 1 GeV with 105 primaries, colliding with a 5 mm Au target with Z = 79 and ρ = 19.3 g/cm3 and a 5 mm Al target with Z = 13 and ρ = 2.7 g/cm3. We disabled the field updater and weighting procedure in the PIC code, and enabled only the particle pusher and bremsstrahlung MC module. The electron and photon spectra were found to be in good agreement with the Geant4 results, except for a slightly higher photon emission in the high-energy tail (which is due to the difference in the cross-section data). Figure 13 displays the spectra of electrons and photons from a 100 MeV electron bunch normally incident onto the Al and Au slabs, and similar distributions for a 1 GeV electron bunch are shown in Fig. 14.
F. Vacuum birefringence
For an example of the VB effect, see Fig. 16.
IV. FRAMEWORK OF SLIPS
These physical processes have been incorporated into a spin-resolved laser–plasma interaction simulation code, known as SLIPs. The data structure and framework layout are illustrated in Figs. 17 and 18.
As depicted in Fig. 17, SLIPs utilizes a toml file to store simulation information, which is then parsed into a SimInfo structure that includes domainInfo, speciesInfo, boundaryInfo, laserInfo, pusherInfo, and other metadata. Subsequently, this metadata are employed to generate a SimBox that comprises all ParticleList and Fields, and to define the FieldSolver and EOMSolver and initialize QED processes.
The internal data structure of SLIPs is constructed using the open-source numerical library, Armadillo C++.83,84 String expressions are parsed using the ExprTk library.85 The data are then dumped using serial-hdf5 and merged with external Python scripts to remove ghost cells.
The spin-resolved processes, i.e., those tagged as Spin-QED in the diagram in Fig. 18, are implemented in conjunction with the Lorentz equation. In the coding, the Spin-QED part is arranged as a sequential series of processes. For example, Lorentz and T-BMT are followed by radiative correction, VB, NBW, and NCS with bremsstrahlung: and Bremss.
V. POLARIZED PARTICLE SIMULATIONS
A. Polarized electron/positron simulation
To simulate the generation of spin-polarized electrons, we utilized an elliptically polarized laser with an intensity a0 = 30, a wavelength λ0 = 1 μm, and an ellipticity ay,0/ax,0 = 3%. This laser was directed toward an ultrarelativistic electron bunch with an energy of 10 GeV, which was produced through laser-wakefield acceleration. The resulting polarized electrons are depicted in Fig. 19, and show good agreement with the previously published results in Ref. 25.
B. Polarized γ-photons via NCS
The polarization state of emitted photons can be determined in spin/polarization-resolved NCS. Here, following Ref. 25, we utilized a linearly polarized (LP) laser to collide with an unpolarized electron bunch to generate LP γ-photons. Additionally, we used an LP laser to collide with a longitudinally polarized electron bunch to generate circularly polarized (CP) γ-photons, which were also observed in a previous study.12 The final polarization states of LP and CP γ-photons are presented in Figs. 20 and 21, respectively.
C. Laser–plasma interactions
Finally, we present a simulation result demonstrating the interaction between an ultraintense laser with a normalized intensity a0 = 1000 and a fully ionized 2 μm thick aluminum target. Note that this configuration, previously examined in Ref. 86 with a thickness of 1 μm, employs a thicker target in the present study to enhance the SF-QED processes. When the laser is directed toward a solid target, the electrons experience acceleration and heating due to the laser and plasma fields. As high-energy electrons travel through the background field, they can emit γ-photons via NCS. The EM field distribution and number densities of target electrons, NBW positrons, and NCS γ-photons are shown in Fig. 22, all of which show good consistency with Ref. 86. The laser is linearly polarized along the y direction, indicating that the polarization frame is mainly in the y–z plane with two polarization bases and , where denotes the momentum direction of the photon. The polarization angle-dependence observed in this study is consistent with previous results in the literature. However, the average linear polarization degree is ∼60% , as illustrated in Figs. 23(b) and 23(d). Notably, low-energy photons contribute primarily to the polarization, as demonstrated in Figs. 23(a) and 23(c). Additionally, during the subsequent NBW process, the self-generated strong magnetic field couples with the laser field dominating the positrons’ SQA. As a result, the positrons’ polarization is aligned with the z direction, contingent on their momentum direction, as shown in Fig. 24. These findings constitute a novel contribution to the investigation of polarization-resolved laser–plasma interactions.
Computer simulation techniques for laser–plasma interactions are constantly evolving, not only in terms of the accuracy of high-order or explicit/implicit algorithms, but also in the complexity of new physics with more degrees of freedom. The rapid development of ultraintense laser techniques not only provides opportunities for experimental verification of SF-QED processes in the high-energy-density regime (which serves as a micro-astrophysics laboratory), but also presents challenges to theoretical analysis. The introduction of Spin-QED into widely accepted PIC algorithms may address this urgent demand and pave the way for studies in laser-QED physics, laser–nuclear physics (astrophysics), and even physics beyond the Standard Model.
The work is supported by the National Natural Science Foundation of China (Grant Nos. 12275209, 12022506, and U2267204), the Open Foundation of the Key Laboratory of High Power Laser and Physics, Chinese Academy of Sciences (Grant No. SGKF202101), the Foundation of Science and Technology on Plasma Physics Laboratory (Grant No. JCKYS2021212008), and the Shaanxi Fundamental Science Research Project for Mathematics and Physics (Grant No. 22JSY014).
Conflict of Interest
The authors have no conflicts to disclose.
Feng Wan: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Chong Lv: Formal analysis (equal); Resources (equal); Validation (equal); Writing – review & editing (equal). Kun Xue: Software (equal); Writing – review & editing (equal). Zhen-Ke Dou: Software (equal); Writing – review & editing (equal). Qian Zhao: Writing – review & editing (equal). Mamutjan Ababekri: Conceptualization (equal); Writing – review & editing (equal). Wen-Qing Wei: Writing – review & editing (equal). Zhong-Peng Li: Writing – review & editing (equal). Yong-Tao Zhao: Conceptualization (equal); Writing – review & editing (equal). Jian-Xing Li: Conceptualization (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal).
The data supporting this study's findings are available from the corresponding author upon reasonable request.