Laser–matter interactions can trigger strong-field quantum-electrodynamics (SF-QED) processes when the laser intensity I₀ reaches or exceeds 10²² W/cm².1,2 For example, when the laser intensity is of the order of 10²¹–10²² W/cm², i.e., the normalized peak laser field strength parameter a₀ ≡ eE₀/m_ecω₀ ∼ 10, electrons can be accelerated to GeV energies3,4 (with Lorentz factor γ_e ∼ 10³ or higher) in a centimeter-long gas plasma, where −e and m_e are the charge and mass of the electron, E₀ and ω₀ are the electric field strength and angular frequency of the laser, and c is the speed of light in vacuum (here, for convenience, it is assumed that ω₀ = 2πc/λ₀ and that the wavelength of the laser is λ₀ = 1 μm). When the laser is reflected by a plasma mirror and collides with the accelerated electron bunch, the transverse electromagnetic (EM) field in the electron’s instantaneous frame can reach the order of a′ ≃ 2γa₀ ∼ 10⁴–10⁵. Such a field strength is close to the QED critical field strength (Schwinger critical field strength) $E_{\text{Sch}}\equiv m_e^2{c}^3/e\hslash$, i.e., a_Sch = m_ec²/ℏω₀ ≃ 4.1 × 10⁵, within one or two orders of magnitude. In this regime, the probabilities of nonlinear QED processes are comparable to those of linear ones, and depend on three parameters as W(χ, f, g), with$$\begin{gathered} \begin{array}{ll}\hfill \chi & \equiv \frac{e\hslash \sqrt{{\left(F_{\mu \nu }{p}^{\mu }\right)}^2}}{m^3{c}^4}~\frac{a^{\prime }}{a_{\text{Sch}}}, f\equiv \frac{e^2{\hslash }^2{F}_{\mu \nu }{F}^{\mu \nu }}{4m^4{c}^6}~\frac{{{\mathit{a}}_E}^2-{{\mathit{a}}_B}^2}{4{a^2}_{\text{Sch}}},\hfill \\ \hfill & g\equiv \frac{e^2{\hslash }^2{F}_{\mu \nu }{F}^{\mu \nu *}}{4m^4{c}^6}~\frac{{{\mathit{a}}_E}^2\mathbf{·}{{\mathit{a}}_B}^2}{4{a^2}_{\text{Sch}}}, \\ \hfill \end{array} \end{gathered}$$where a_E and a_B denote the normalized field strengths of the electric and magnetic components, respectively.5,6 For most cases of weak-field (a₀ ≪ a_Sch) conditions, f, g ≪ χ², and W(χ, f, g) ∼ W(χ), i.e., the probability depends on only a single parameter χ. For electrons and positrons, nonlinear Compton scattering (NCS, e + nω_L → e′ + ω_γ) is the dominant nonlinear QED process in the strong-field regime, whereas for photons, nonlinear Breit–Wheeler (NBW) pair production (ω_γ + nω_L → e⁺ + e⁻) is the dominant process, where ω_L and ω_γ denote the laser photon and the emitted γ-photon, respectively, and n is the photon absorption number.

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 m_ec, velocity v with c, energy ɛ with m_ec², EM fields E and B with m_ecω/e, force F with m_ecω, charge q with e, charge density ρ with k³e, and current density J with k³ec, where λ and ω = 2πc/λ are the reference wavelength and frequency, respectively.

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

When radiation reaction is weak (the radiation power is much smaller than the energy gain power), the motion of charged particles is governed by the Newton–Lorentz equations:$$\frac{d\mathbf{p}}{dt}=\frac{q}{m}(\mathbf{E}+\mathit{\beta }\times \mathbf{B}),$$(1)$$\frac{d\mathbf{x}}{dt}=\frac{\mathbf{p}}{\gamma },$$(2)where p ≡ γmv, x, q, m, γ, v, and β ≡ v/c are the momentum, position, charge, mass, Lorentz factor, velocity, and normalized velocity of a particle, respectively. These coupled equations are discretized using a leapfrog algorithm as$$\frac{{\mathbf{p}}^{n+1/2}-{\mathbf{p}}^{n-1/2}}{\mathrm{\Delta }t}=\frac{q}{m}\left({\mathbf{E}}^n+\frac{{\mathbf{p}}^n}{{\gamma }^n}\times {\mathbf{B}}^n\right),$$(3)$$\frac{{\mathbf{x}}^{n+1}-{\mathbf{x}}^n}{\mathrm{\Delta }t}={\mathbf{v}}^{n+1/2},$$(4)and are solved using the standard Boris rotation:31–33$${\mathbf{p}}^{n-1/2}={\mathbf{p}}^{-}-\frac{q\mathrm{\Delta }t}{2m}{\mathbf{E}}^n,$$(5)$${\mathbf{p}}^{n+1/2}={\mathbf{p}}^{+}+\frac{q\mathrm{\Delta }t}{2m}{\mathbf{E}}^n,$$(6)$${\mathbf{p}}^{\prime }={\mathbf{p}}^{-}+{\mathbf{p}}^{-}\times \mathit{\tau },$$(7)$${\mathbf{p}}^{+}={\mathbf{p}}^{-}+{\mathbf{p}}^{\prime }\times \mathit{\varsigma },$$(8)$$\mathit{\tau }=\frac{q\mathrm{\Delta }t}{2m{\gamma }^n}{\mathbf{B}}^n,$$(9)$$\mathit{\varsigma }=\frac{2\mathit{\tau }}{1+{\mathit{\tau }}^2},$$(10)where ${\gamma }^n=\sqrt{1+{({\mathbf{p}}^{-})}^2}=\sqrt{1+{({\mathbf{p}}^{+})}^2}$. The updates in momentum and position are asynchronized by half a time step, i.e., a leapfrog algorithm is used here. This leapfrog algorithm ensures the self-consistency of the momentum and position evolutions.

In ultraintense laser–plasma interactions, the plasma particles are assumed to be distributed in vacuum and immersed in the EM field. Therefore, the field evolution is governed by the Maxwell equations in vacuum with sources. After normalization, the Maxwell equations are given in differential form as$$\nabla \cdot \mathbf{E}=\rho ,$$(11)$$\nabla \cdot \mathbf{B}=0$$(12)$$\nabla \times \mathbf{E}=-\frac{\partial \mathbf{B}}{\partial t},$$(13)$$\nabla \times \mathbf{B}=\frac{\partial \mathbf{E}}{\partial t}+\mathbf{J}.$$(14)

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 E_y and B_z components as examples here:

1D case (squeezing the y and z directions):$$\begin{array}{c}{\left.\frac{E_y^{n+1}-E_y^n}{\mathrm{\Delta }t}\right|}_{i+1/2}={\left.-\frac{B_{i+1}-B_i}{\mathrm{\Delta }x}\right|}_z^{n+1/2}-J_{y,i+1/2}^{n+1/2},\\ {\left.\frac{B_z^{n+1/2}-B_z^{n-1/2}}{\mathrm{\Delta }t}\right|}_i={\left.-\frac{E_{i+1/2}-E_{i-1/2}}{\mathrm{\Delta }x}\right|}_y^n;\end{array}$$(15)2D case (squeezing the z direction):$$\begin{array}{ll}\hfill {\left.\frac{E_y^{n+1}-E_y^n}{\mathrm{\Delta }t}\right|}_{i+1/2,j}& ={\left.-\frac{B_{i+1,j}-B_{i,j}}{\mathrm{\Delta }x}\right|}_z^{n+1/2}-J_{y,i+1/2,j}^{n+1/2},\hfill \\ \hfill {\left.\frac{B_z^{n+1/2}-B_z^{n-1/2}}{\mathrm{\Delta }t}\right|}_{i+1/2,j}& ={\left.-\frac{E_{i+1/2,j}-E_{i+1/2,j}}{\mathrm{\Delta }x}\right|}_y^n\hfill \\ \hfill & {\left.+\frac{E_{i+1/2,j+1/2}-E_{i+1/2,j-1/2}}{\mathrm{\Delta }y}\right|}_x^{n+1/2};\hfill \end{array}$$(16)

3D case:$$\begin{array}{ll}\hfill {\left.\frac{E_y^{n+1}-E_y^n}{\mathrm{\Delta }t}\right|}_{i+1/2,j,k+1/2}& ={\left.-\frac{B_{i+1,j,k+1/2}-B_{i,j,k+1/2}}{\mathrm{\Delta }x}\right|}_z^{n+1/2}\hfill \\ \hfill & {\left.+\frac{B_{i+1/2,j,k+1}-B_{i+1/2,j,k}}{\mathrm{\Delta }z}\right|}_x^{n+1/2}\hfill \\ \hfill & -J_{y,i+1/2,j,k+1/2}^{n+1/2},\hfill \\ \hfill {\left.\frac{B_z^{n+1/2}-B_z^{n-1/2}}{\mathrm{\Delta }t}\right|}_{i,j,k+1/2}& ={\left.-\frac{E_{i+1/2,j,k}-E_{i-1/2,j,k}}{\mathrm{\Delta }x}\right|}_y^n\hfill \\ \hfill & {\left.+\frac{E_{i+1/2,j+1,k}-E_{i+1/2,j,k}}{\mathrm{\Delta }y}\right|}_x^n.\hfill \end{array}$$(17)

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.

We calculate the charge current density using Esirkepov’s method, which conserves charge by satisfying Gauss’s law35$${\partial }_t\rho +\nabla \cdot \mathbf{J}=0,$$(18)and removes the need for Coulomb correction.19 This algorithm computes the charge density at time steps $t-\frac{1}{2}\mathrm{\Delta }t$ and $t+\frac{1}{2}\mathrm{\Delta }t$ on each grid cell from the particle positions and velocities, i.e.,$${\rho }_{i,j,k}^{n+1/2}=\frac{1}{\mathrm{\Delta }V}\sum _{r}W\left({\mathbf{x}}_r^n+\frac{1}{2}{\mathbf{v}}^{n+1/2}\mathrm{\Delta }t\right)q_r,$$(19)$${\rho }_{i,j,k}^{n-1/2}=\frac{1}{\mathrm{\Delta }V}\sum _{r}W\left({\mathbf{x}}_r^n-\frac{1}{2}{\mathbf{v}}^n\mathrm{\Delta }t\right)q_r,$$(20)$${\delta }^n\rho ={\rho }^{n+1/2}-{\rho }^{n-1/2},$$(21)where r denotes the particle index, |x_r − x_i,j,k| ≤ (Δx, Δy, Δz), and ΔV = ΔxΔyΔz is the cell volume. We then interpolate the charge density to the current grid to obtain the current density; see Ref. 35 for more details.

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 x_r, we find its nearest grid point ${(i,j,k)}^g=\mathrm{fl}\mathrm{o}\mathrm{o}\mathrm{r}({\mathbf{x}}_r/\mathrm{\Delta }\mathbf{x}+\frac{1}{2})$ and its nearest half grid point (i, j, k) = floor(x_r/Δ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.

The weighting function W in the current and force deposition is determined by the form factor (shape factor) of the macroparticle, which is a key concept in modern PIC algorithms. The form factor gives the macroparticle a finite size (composed of thousands of real particles) and reduces the nonphysical collisions.19 Various particle shape function models have been proposed, such as the nearest grid point (NGP) and cloud-in-cell (CIC) methods. The NGP and CIC methods use the nearest one and two grid fields as the full contribution, respectively. Higher orders of particle shape function can suppress unphysical noise and produce smoother results. We use a triangle shape function (triangular shape cloud, TSC) in each dimension:35$$W_{\text{spline}}=\left\{\begin{array}{cc}\frac{3}{4}-{\delta }^2,\hfill & \mathrm{f}\mathrm{o}\mathrm{r}j,\hfill \\ \frac{1}{2}\left(\frac{1}{2}\pm \delta \right),\hfill & \mathrm{f}\mathrm{o}\mathrm{r}j\pm 1,\hfill \end{array}\right.$$(22)where δ = (x − X_j)/Δx, x is the particle position, j is the nearest grid/half grid number, and X_j ≡ jΔx. We obtain higher-dimensional functions as products of 1D shape functions in each dimension: W_2D(i, j) = W_x(i)W_y(j) and W_3D(i, j, k) = W_x(i)W_y(j)W_z(k).

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.

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.

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

The LL equation can be employed when the radiation is relatively weak (weak radiation reaction, χ_e ≪ 10⁻²),37 and, in Gaussian units, takes the form$$\begin{array}{ll}\hfill & {\mathbf{F}}_{\text{RR,classical}}\hfill \\ \hfill & =\frac{2e^3}{3m{c}^3}\left\{\gamma \left[\left(\frac{\partial }{\partial t}+\frac{\mathbf{p}}{\gamma m}\cdot \nabla \right)\mathbf{E}+\frac{\mathbf{p}}{\gamma mc}\times \left(\frac{\partial }{\partial t}+\frac{\mathbf{p}}{\gamma m}\cdot \nabla \right)\mathbf{B}\right]\right.\hfill \\ \hfill & +\frac{e}{mc}\left[\mathbf{E}\times \mathbf{B}+\frac{1}{\gamma mc}\mathbf{B}\times (\mathbf{B}\times \mathbf{p})+\frac{1}{\gamma mc}\mathbf{E}(\mathbf{p}\cdot \mathbf{E})\right]\hfill \\ \hfill & \left.-\frac{e\gamma }{m^2{c}^2}\mathbf{p}\left[{\left(\mathbf{E}+\frac{\mathbf{p}}{\gamma mc}\times \mathbf{B}\right)}^2-\frac{1}{{\gamma }^2{m}^2{c}^2}{(\mathbf{E}\cdot \mathbf{p})}^2\right]\right\}.\hfill \end{array}$$(23)

The dimensionless form of this equation is$$\begin{array}{ll}\hfill {\mathbf{F}}_{\text{RR,classical}}& =\frac{2}{3}{\alpha }_f{\xi }_L\left\{\gamma \left[\left(\frac{\partial }{\partial t}+\frac{\mathbf{p}}{\gamma }\cdot \nabla \right)\mathbf{E}+\frac{\mathbf{p}}{\gamma }\times \left(\frac{\partial }{\partial t}+\frac{\mathbf{p}}{\gamma }\cdot \nabla \right)\mathbf{B}\right]\right.\hfill \\ \hfill & +\left[\mathbf{E}\times \mathbf{B}+\frac{1}{\gamma }\mathbf{B}\times (\mathbf{B}\times \mathbf{p})+\mathbf{E}(\mathbf{p}\cdot \mathbf{E})\right]\hfill \\ \hfill & \left.-\gamma \mathbf{p}\left[{\left(\mathbf{E}+\frac{\mathbf{p}}{\gamma }\times \mathbf{B}\right)}^2-\frac{1}{{\gamma }^2}{(\mathbf{E}\cdot \mathbf{p})}^2\right]\right\},\hfill \end{array}$$(24)where α_f = e²/cℏ is the fine structure constant and ξ_L = ℏω/m_ec² is the normalized reference photon energy. In the case of an ultraintense laser interacting with a plasma, the dominant contribution comes from the last two terms.43 In the ultrarelativistic limit, only the third term dominates the contribution, and the radiation reaction force can be simplified as$${\mathbf{F}}_{\text{RR,cl}}\simeq \frac{2}{3}{\alpha }_f\frac{{\chi }_e^2}{{\xi }_L}\mathit{\beta },$$(25)where$$\begin{array}{ll}\hfill {\chi }_e& =\frac{e\hslash }{m^3{c}^4}\sqrt{|F^{\mu \nu }{p}_{\nu }{|}^2}\hfill \\ \hfill & \equiv {\xi }_L{\gamma }_e\sqrt{{(\mathbf{E}+\mathit{\beta }\times \mathbf{B})}^2-{[\mathit{\beta }\cdot (\mathit{\beta }\cdot \mathbf{E})]}^2}\hfill \\ \hfill & \simeq {\gamma }_e{E}_{\perp }{\xi }_L(1-\cos \theta )\hfill \end{array}$$is a nonlinear quantum parameter signifying the strength of the NCS, with θ denoting the angle between the electron momentum and the EM field wavevector and E_⊥ denoting the perpendicular component of the electric field. This reduced form gives the importance of the radiation reaction when one estimates the ratio between F_RR and the Lorentz force F_L:$$\begin{array}{ll}\hfill \mathcal{R}& \equiv |F_{\text{RR}}|/|F_{\text{L}}|\hfill \\ \hfill & \sim \frac{2}{3}{\alpha }_f{\gamma }_e{\chi }_e\simeq 2\times 10^{-8}{a}_0{\gamma }_e^2(\text{for\hspace{0.17em}wavelength\hspace{0.17em}}1\mu \text{m}).\hfill \end{array}$$(26)Clearly, once ${\gamma }_e^2{a}_0\gtrsim 10^6$, the radiation reaction force should be considered.

The LL equation is only applicable when the radiation reaction force is much weaker than the Lorentz force, or the radiation per laser period is much smaller than m_ec².44 Once χ_e is larger than 10⁻², the quantum nature of the radiation dominates the process. On the one hand, the radiation spectrum will be suppressed and deviate from the radiation force in the LL equation; on the other hand, the radiation will be stochastic and discontinuous. However, when the stochasticity is not relevant for detection and one only cares about the average effect (integrated spectra), a correction to the radiation force can be made, i.e., a quantum correction45–48$${\mathbf{F}}_{\text{RR,quantum}}=q(\chi ){\mathbf{F}}_{\text{RR,classical}},$$(27)where$$q(\chi )=\frac{I_{\text{QED}}}{I_{\text{C}}},$$(28)$$I_{\text{QED}}=m{c}^2\int c(k\cdot k^{\prime })\frac{d{W}_{fi}}{d\eta d{r}_0}d{r}_0,$$(29)$$I_{\text{C}}=\frac{2e^4{E}^{\prime 2}}{3m^2{c}^3},$$(30)with W_fi being the radiation probability,49η = k₀z − ω₀t, r₀ = 2(k · k′)/3χ(k · p_i), and E′ is the electric field in the instantaneous frame of the electron. p_i is the four-momentum of the electron before radiation. k and k′ are the four-wavevectors of the local EM field and the radiated photon, respectively. See q(χ) in Fig. 4. Here, the ratio between the QED radiation power and the classical one, i.e., the re-scaling factor q(χ), is the same as the factor in Ref. 44:$$q({\chi }_e)\approx \frac{1}{{\left[1+4.8(1+{\chi }_e)\ln (1+1.7{\chi }_e)+2.44{\chi }_e^2\right]}^{2/3}},$$(31)or$$q({\chi }_e)\approx \frac{1}{{(1+8.93{\chi }_e+2.41{\chi }_e^2)}^{2/3}}.$$(32)

In the ultrarelativistic limit, the following alternative formula can be employed:23,50$${\mathbf{F}}_{\text{RR,quantum}}=q(\chi )P_{\text{cl}}{\chi }_e^2\mathit{\beta }/{\mathit{\beta }}^{2}c.$$(33)Clearly, once χ ≳ 10⁻², the quantum-corrected version should be used.

Here, we plug the radiative correction (either classical or quantum corrected version) into the standard Boris pusher as follows:43$$\frac{{\mathbf{p}}^{n+1/2}-{\mathbf{p}}^{n-1/2}}{\mathrm{\Delta }t}={\mathbf{F}}^n={\mathbf{F}}_L^n+{\mathbf{F}}_R^n.$$(34)First we use the Boris step$$\frac{{\mathbf{p}}_L^{n+1/2}-{\mathbf{p}}_L^{n-1/2}}{\mathrm{\Delta }t}={\mathbf{F}}_L^n,$$(35)and then use the radiative correction push$$\frac{{\mathbf{p}}_R^{n+1/2}-{\mathbf{p}}_R^{n-1/2}}{\mathrm{\Delta }t}={\mathbf{F}}_R^n,$$(36)where ${\mathbf{p}}_L^{n-1/2}={\mathbf{p}}_R^{n-1/2}={\mathbf{p}}^{n-1/2}$, and the final momentum is given by$${\mathbf{p}}^{n+1/2}={\mathbf{p}}_L^{n+1/2}+{\mathbf{p}}_R^{n+1/2}-{\mathbf{p}}^{n-1/2}={\mathbf{p}}_L^{n+1/2}+{\mathbf{F}}_R^n\mathrm{\Delta }t.$$(37)With this algorithm, the Boris pusher is realized.

Figure 5 presents a comparison between dynamics calculated using different solvers. For the Lorentz equation without radiation, the particle momentum and energy are given analytically by51$$\mathbf{p}(\tau )={\mathbf{p}}_0-\mathbf{A}(\tau )+\widehat{\mathbf{k}}\frac{{\mathbf{A}}^2(\tau )-2{\mathbf{p}}_0\cdot \mathbf{A}(\tau )}{2({\gamma }_0-{\mathbf{p}}_0\cdot \widehat{\mathbf{k}})}$$(38)$$\gamma (\tau )={\gamma }_0+\frac{{\mathbf{A}}^2(\tau )-2{\mathbf{p}}_0\cdot \mathbf{A}(\tau )}{2({\gamma }_0-{\mathbf{p}}_0\cdot \widehat{\mathbf{k}})}$$(39)where $\mathbf{A}\left(\tau \right)=-{\int }_{{\tau }_0}^{\tau }\mathbf{E}\left({\tau }^{\prime }\right)d{\tau }^{\prime }$ is the external field vector potential, τ is the proper time, $\widehat{\mathbf{k}}$ is the normalized wavevector of the field, and γ, p, and γ₀, p₀ are the instantaneous and initial (subscript 0) Lorentz factor and momentum of the particle, respectively. For a planewave with a temporal profile, the momentum and energy gain vanish as A(∞) = A(−∞) = 0. The planewave solution with radiation reaction can be found in Ref. 52. However, no explicit solution exists when the quantum correction term is included, as shown in Fig. 5.

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.

The nonradiative spin dynamics of an electron are given by$$\begin{array}{ll}\hfill {\left(\frac{d\mathbf{S}}{dt}\right)}_T& =\mathbf{S}\times \mathbf{\Omega }\hfill \\ \hfill & \equiv \mathbf{S}\times \left[-\left(\frac{g}{2}-1\right)\frac{{\gamma }_e}{{\gamma }_e+1}\left(\mathit{\beta }\cdot \mathbf{B}\right)\cdot \mathit{\beta }\right.\hfill \\ \hfill & \left.+\left(\frac{g}{2}-1+\frac{1}{{\gamma }_e}\right)\mathbf{B}-\left(\frac{g}{2}-\frac{{\gamma }_e}{{\gamma }_e+1}\right)\mathit{\beta }\times \mathbf{E}\right],\hfill \end{array}$$(40)where E and B are the normalized electric and magnetic fields and g is the electron Landé factor. Since the this equation is a pure rotation around the precession frequency of Ω, Boris rotation is greatly preferable to other solvers for ordinary differential equations (Runge–Kutta, etc.). Here, Ω plays the role of B/γ in Eqs. (3) and (5)–(10). For other particle species, the appropriate charge, mass, and Landé factor should be employed.

When radiation damping is no longer negligible, the radiation can also affect the spin dynamics. In the weak radiation regime, this radiation-induced modification of the spin dynamics can be handled in a similar way as in the LL equation. Thus, the modified version of the T-BMT equation, the radiative T-BMT equation, is given by$$\frac{d\mathbf{S}}{dt}={\left(\frac{d\mathbf{S}}{dt}\right)}_T+{\left(\frac{d\mathbf{S}}{dt}\right)}_R,$$(41)with the first (labeled with “T”) and second (labeled with “R”) terms corresponding to the nonradiative precession in Eq. (40) and the radiative correction, respectively. The radiative term is given by$${\left(\frac{d\mathbf{S}}{dt}\right)}_R=-P\left[{\psi }_1(\chi )\mathbf{S}+{\psi }_2(\chi )(\mathbf{S}\cdot \mathit{\beta })\mathit{\beta }+{\psi }_3(\chi ){\widehat{\mathbf{n}}}_B\right].$$(42)Here,$$\begin{array}{l}\hfill P=\frac{{\alpha }_f}{\sqrt{3}\pi {\gamma }_e{\xi }_L},{\psi }_1({\chi }_e)={\int }_0^{\infty }{u}^{″}du{\mathrm{K}}_{2/3}(u^{\prime }),\\ \hfill {\psi }_2({\chi }_e)={\int }_0^{\infty }{u}^{″}du{\int }_{u^{\prime }}^{\infty }dx{\mathrm{K}}_{1/3}(x)-{\psi }_1({\chi }_e),\\ \hfill {\psi }_3(\chi )={\int }_0^{\infty }{u}^{″}du{\mathrm{K}}_{1/3}(u^{\prime }),u^{\prime }=\frac{2u}{3{\chi }_e},u^{″}=\frac{u^2}{{(1+u)}^3},\end{array}$$where K_n is the nth-order modified Bessel function of the second kind, ${\widehat{\mathbf{n}}}_B=\mathit{\beta }\times \widehat{\mathbf{a}}$, and β and $\widehat{\mathbf{a}}$ denote the normalized velocity and acceleration vectors, respectively.58,59

The simulation algorithms for spin precession are quite similar to those for the EOM (the Lorentz equation and radiative EOM), namely, the LL/MLL equations. Therefore, the T-BMT equation is simulated via Boris rotation without the pre- and post-acceleration terms, and with only the rotation term Ω. In SLIPs, a standard Boris algorithm is used:$${\mathbf{S}}^{\prime }={\mathbf{S}}^{n-1/2}+{\mathbf{S}}^{n-1/2}\times \mathbf{t},$$(43)$${\mathbf{S}}^{n+1/2}={\mathbf{S}}^{n-1/2}+{\mathbf{S}}^{\prime }\times \mathbf{o},$$(44)$$\mathbf{t}=\frac{q\mathrm{\Delta }t}{2}{\mathbf{\Omega }}^n,$$(45)$$\mathbf{o}=\frac{2\mathbf{t}}{1+t^2}.$$(46)

For the radiative T-BMT equation, there will be an extra term (dS/dt)_R, which is equivalent to the electric field term in the Lorentz equation. Therefore, the straightforward algorithm is given by$${\mathbf{S}}_T^{n-1/2}={\mathbf{S}}^{n-1/2}+\frac{\mathrm{\Delta }t}{2}{\left(\frac{d\mathbf{S}}{dt}\right)}_R,$$(47)Boris T-BMT Eqs. (43)–(46),$${\mathbf{S}}^{n+1/2}={\mathbf{S}}_T^{n+1/2}+\frac{\mathrm{\Delta }t}{2}{\left(\frac{d\mathbf{S}}{dt}\right)}_R.$$(48)

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 S_x and S_y. The radiation reaction mainly affects S_z. In the case without radiation reaction (case A), S_x and S_y oscillate owing to precession and are conserved in Fig. 6(d). In the case with only spin radiation reaction (case B), S_x is strongly damped by the term (dS/dt)_R. S_y and S_z 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 S_x and S_z [see Fig. 6(c) for a comparison of cases B, D, and C in terms of S_z amplitude]. Simultaneously, the precession term (dS/dt)_T ∝ B/γ_e grows with decreasing γ_e, which amplifies the oscillation of S_y, as shown by the contrast between B (Lorentz), D (MLL), and C (LL) in Fig. 6(b).

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.

When the laser intensity a₀ and the electron energy γ_e are such that the locally constant cross-field approximation (LCFA) is valid, i.e., a₀ ≫ 1, χ_e ≳ 1, the polarization- and spin-resolved emission rate for the NCS is given by12,15,62$$\frac{d^2{W}_{fi}}{dudt}=\frac{W_R}{2}(F_0+{\xi }_1{F}_1+{\xi }_2{F}_2+{\xi }_3{F}_3),$$(49)where the photon polarization is represented by the Stokes parameters (ξ₁, ξ₂, ξ₃), defined with respect to the axes ${\widehat{\mathbf{P}}}_1=\widehat{\mathbf{a}}-\widehat{\mathbf{n}}(\widehat{\mathbf{n}}\cdot \widehat{\mathbf{a}})$ and ${\widehat{\mathbf{P}}}_2=\widehat{\mathbf{n}}\times {\widehat{\mathbf{P}}}_1$,63 with the photon emission direction $\widehat{\mathbf{n}}={\mathbf{p}}_e/|{\mathbf{p}}_e|$ along the momentum p_e of the ultrarelativistic electron. The variables introduced in Eq. (49) are as follows:$$\begin{array}{ll}\hfill F_0& =-{(2+u)}^2\left[{\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{K}}_{1/3}(u^{\prime })-2{\mathrm{K}}_{2/3}(u^{\prime })\right](1+S_{if})\hfill \\ \hfill & +u^2(1-S_{if})\left[{\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{K}}_{1/3}(u^{\prime })+2{\mathrm{K}}_{2/3}(u^{\prime })\right]+2u^2{S}_{if}{\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{K}}_{1/3}(u^{\prime })\hfill \\ \hfill & -(4u+2u^2)({\mathbf{S}}_f+{\mathbf{S}}_i)\cdot \left[\widehat{\mathbf{n}}\times \widehat{\mathbf{a}}\right]{\mathrm{K}}_{1/3}(u^{\prime })-2u^2({\mathbf{S}}_f-{\mathbf{S}}_i)\hfill \\ \hfill & \cdot \left[\widehat{\mathbf{n}}\times \widehat{\mathbf{a}}\right]{\mathrm{K}}_{1/3}(u^{\prime })-4u^2\left[{\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{K}}_{1/3}(u^{\prime })-{\mathrm{K}}_{2/3}(u^{\prime })\right]\hfill \\ \hfill & \times ({\mathbf{S}}_i\cdot \widehat{\mathbf{n}})({\mathbf{S}}_f\cdot \widehat{\mathbf{n}}),\hfill \end{array}$$(50)$$\begin{array}{ll}\hfill F_1& =-2u^2{\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{K}}_{1/3}(u^{\prime })\left\{({\mathbf{S}}_i\cdot \widehat{\mathbf{a}}){\mathbf{S}}_f\cdot \left[\widehat{\mathbf{n}}\times \widehat{\mathbf{a}}\right]+({\mathbf{S}}_f\cdot \widehat{\mathbf{a}}){\mathbf{S}}_i\cdot [\widehat{\mathbf{n}}\times \widehat{\mathbf{a}}]\right\}\hfill \\ \hfill & +4u\left[({\mathbf{S}}_i\cdot \widehat{\mathbf{a}})(1+u)+({\mathbf{S}}_f\cdot \widehat{\mathbf{a}})\right]{\mathrm{K}}_{1/3}(u^{\prime })\hfill \\ \hfill & +2u(2+u)\widehat{\mathbf{n}}\cdot [{\mathbf{S}}_f\times {\mathbf{S}}_i]{\mathrm{K}}_{2/3}(u^{\prime }),\hfill \end{array}$$(51)$$\begin{array}{ll}\hfill F_2& =-\left\{2u^2\left\{({\mathbf{S}}_i\cdot \widehat{\mathbf{n}}){\mathbf{S}}_f\cdot \left[\widehat{\mathbf{n}}\times \widehat{\mathbf{a}}\right]+({\mathbf{S}}_f\cdot \widehat{\mathbf{n}}){\mathbf{S}}_i\cdot \left[\widehat{\mathbf{n}}\times \widehat{\mathbf{a}}\right]\right\}\right.\hfill \\ \hfill & \left.+2u(2+u)\widehat{\mathbf{a}}\cdot [{\mathbf{S}}_f\times {\mathbf{S}}_i]\right\}{\mathrm{K}}_{1/3}(u^{\prime })-4u\left[({\mathbf{S}}_i\cdot \widehat{\mathbf{n}})\right.\hfill \\ \hfill & \left.+({\mathbf{S}}_f\cdot \widehat{\mathbf{n}})(1+u)\right]{\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{K}}_{1/3}(u^{\prime })+4u(2+u)\hfill \\ \hfill & \times \left[({\mathbf{S}}_i\cdot \widehat{\mathbf{n}})+({\mathbf{S}}_f\cdot \widehat{\mathbf{n}})\right]{\mathrm{K}}_{2/3}(u^{\prime }),\hfill \end{array}$$(52)$$\begin{array}{ll}\hfill F_3& =4\left[1+u+\left(1+u+\frac{1}{2}{u}^2\right)S_{if}-\frac{1}{2}{u}^2({\mathbf{S}}_i\cdot \widehat{\mathbf{n}})({\mathbf{S}}_f\cdot \widehat{\mathbf{n}})\right]{\mathrm{K}}_{2/3}(u^{\prime })\hfill \\ \hfill & +2u^2\left\{{\mathbf{S}}_i\cdot \left[\widehat{\mathbf{n}}\times \widehat{\mathbf{a}}\right]{\mathbf{S}}_f\cdot \left[\widehat{\mathbf{n}}\times \widehat{\mathbf{a}}\right]-({\mathbf{S}}_i\cdot \widehat{\mathbf{a}})({\mathbf{S}}_f\cdot \widehat{\mathbf{a}})\right\}{\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{K}}_{1/3}(u^{\prime })\hfill \\ \hfill & -4u\left\{(1+u){\mathbf{S}}_i\left[\widehat{\mathbf{n}}\times \widehat{\mathbf{a}}\right]+{\mathbf{S}}_f\left[\widehat{\mathbf{n}}\times \widehat{\mathbf{a}}\right]\right\}{\mathrm{K}}_{1/3}(u^{\prime }),\hfill \end{array}$$(53)where$$W_R=\frac{{\alpha }_f}{8\sqrt{3}\pi {\xi }_L{(1+u)}^3},u^{\prime }=\frac{2u}{3\chi },u=\frac{{\omega }_{\gamma }}{{\epsilon }_i-{\omega }_{\gamma }},$$ω_γ is the emitted photon energy, ɛ_i is the electron energy before radiation, $\widehat{\mathbf{a}}=\mathbf{a}/|\mathbf{a}|$ is the direction of the electron acceleration a, S_i and S_f are the electron spin vectors respectively before and after radiation (|S_i| = |S_f| = 1), and S_if ≡ S_i · S_f. The function Int K_1/3 is defined as follows:$${\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{K}}_{1/3}(u^{\prime })\equiv {\int }_{u^{\prime }}^{\infty }dz{\mathrm{K}}_{1/3}(z).$$

By summing over the photon polarizations, the electron spin-resolved emission probability can be written as12,15,64$$\begin{array}{ll}\hfill \frac{d^2{W}_{fi}}{dudt}& =W_R\left\{-{(2+u)}^2\left[{\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{K}}_{1/3}(u^{\prime })-2{\mathrm{K}}_{2/3}(u^{\prime })\right](1+S_{if})\right.\hfill \\ \hfill & +u^2\left[{\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{K}}_{1/3}(u^{\prime })+2{\mathrm{K}}_{2/3}(u^{\prime })\right](1-S_{if})\hfill \\ \hfill & +2u^2{S}_{if}{\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{K}}_{1/3}(u^{\prime })-(4u+2u^2)({\mathbf{S}}_f+{\mathbf{S}}_i)[\mathbf{n}\times \widehat{\mathbf{a}}]\hfill \\ \hfill & \times {\mathrm{K}}_{1/3}(u^{\prime })-2u^2({\mathbf{S}}_f-{\mathbf{S}}_i)[\mathbf{n}\times \widehat{\mathbf{a}}]{\mathrm{K}}_{1/3}(u^{\prime })\hfill \\ \hfill & \left.-4u^2\left[{\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{K}}_{1/3}(u^{\prime })-{\mathrm{K}}_{2/3}(u^{\prime })\right]({\mathbf{S}}_i\cdot \mathbf{n})({\mathbf{S}}_f\cdot \mathbf{n})\right\},\hfill \end{array}$$(54)and by summing over the final states S_f, the initial spin-resolved radiation probability is obtained:$$\begin{array}{ll}\hfill \frac{d^2{\bar{W}}_{fi}}{dudt}& =8W_R\left\{-(1+u){\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{K}}_{1/3}(u^{\prime })\right.\hfill \\ \hfill & \left.+(2+2u+u^2){\mathrm{K}}_{2/3}(u^{\prime })-u{\mathbf{S}}_i\cdot \left[\mathbf{n}\times \widehat{\mathbf{a}}\right]{\mathrm{K}}_{1/3}(u^{\prime })\right\}.\hfill \end{array}$$(55)By averaging the electron initial spin, one obtains the widely used radiation probability for the unpolarized initial particles.5,45,65

During the photon emission simulation, the electron/positron spin transitions to either a parallel or antiparallel orientation with respect to the spin quantized axis (SQA), depending on the occurrence of emission. Upon photon emission, the SQA is chosen to obtain the maximum transition probability, which is along the energy-resolved average polarization$${\mathbf{S}}_f^R=\frac{\mathbf{g}}{w+\mathbf{f}\cdot {\mathbf{S}}_i}.$$(56)This is obtained by summing over the photon polarization and retains the dependence on the initial and final electron spin:$$\frac{d^2{W}_{\text{rad}}}{dudt}=W_r(w+\mathbf{f}\cdot {\mathbf{S}}_i+\mathbf{g}\cdot {\mathbf{S}}_f),$$(57)where$$\begin{array}{ll}\hfill & w=-(1+u){\mathrm{K}}_{1/3}({\rho }^{\prime })+(2+2u+u^2){\mathrm{K}}_{2/3}({\rho }^{\prime }),\hfill \\ \hfill & \mathbf{f}=u{\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{K}}_{1/3}({\rho }^{\prime })\widehat{\mathbf{v}}\times \widehat{\mathbf{a}},\hfill \\ \hfill & \mathbf{g}=-(1+u)\left[{\mathrm{K}}_{1/3}({\rho }^{\prime })-2{\mathrm{K}}_{2/3}({\rho }^{\prime })\right]{\mathbf{S}}_i\hfill \\ \hfill & -(1+u)u{\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{K}}_{1/3}({\rho }^{\prime })\widehat{\mathbf{v}}\times \widehat{\mathbf{a}}\hfill \\ \hfill & -u^2\left[{\mathrm{K}}_{1/3}({\rho }^{\prime })-{\mathrm{K}}_{2/3}({\rho }^{\prime })\right]({\mathbf{S}}_i\cdot \widehat{\mathbf{v}})\widehat{\mathbf{v}}.\hfill \end{array}$$Conversely, without emission, the SQA aligns with another SQA.12,66 In both cases, the final spin is determined by assessing the probability density for alignment, either parallel or antiparallel, with the SQA. We account for the stochastic spin flip during photon emission using four random numbers r_1,2,3,4 ∈ [0, 1). The procedure is as follows. First, at each simulation time step Δt, a photon with energy ω_γ = r₁γ_e is emitted if the spin-dependent radiation probability in Eq. (55), $P\equiv d^2{\bar{W}}_{fi}({\chi }_e,r_1,{\gamma }_e,{\mathbf{S}}_i)/dudt\cdot \mathrm{\Delta }t$, meets or exceeds r₂, following the so-called von Neumann rejection method. The final momenta of the electron and photon are given by p_f = (1 − r₁)p_i and k = r₁p_i, respectively. Next, the electron spin flips either parallel (spin-up) or antiparallel (spin-down) to the SQA with probabilities $P_{\text{flip}}\equiv W_{fi}^{\uparrow }/P$ and $W_{fi}^{\downarrow }/P$, respectively, where $W_{fi}^{\uparrow ,\downarrow }\equiv d^2{W}_{fi}^{\uparrow ,\downarrow }/dudt\cdot \mathrm{\Delta }t$ from Eq. (57). In other words, the final spin S_f will flip parallel to the SQA if r₃ < P_flip, and vice versa; see the flow chart of NCS in Fig. 7. In the alternative scenario, i.e., when no photon is emitted, the average final spin is given by$${\bar{\mathbf{S}}}_f=\frac{{\mathbf{S}}_i(1-W\mathrm{\Delta }t)-\mathbf{f}\mathrm{\Delta }t}{1-(W+\mathbf{f}\cdot {\mathbf{S}}_i)\mathrm{\Delta }t},$$where$$W\equiv 16W_R[-(1+u){\mathrm{I}\mathrm{n}\mathrm{t}\mathrm{K}}_{1/3}(u^{\prime })+(2+2u+u^2){\mathrm{K}}_{2/3}(u^{\prime })]$$and $\mathbf{f}\equiv -16W_R\mathbf{n}\times \widehat{\mathbf{a}}{\mathrm{K}}_{1/3}(u^{\prime })$.12,66 Then, the SQA is given by ${\bar{\mathbf{S}}}_f/|{\bar{\mathbf{S}}}_f|$, and the probability for the aligned case is given by $|{\bar{\mathbf{S}}}_f|$ and that for the antiparallel case by $1-|{\bar{\mathbf{S}}}_f|$.

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 ${\widehat{\mathit{\xi }}}^{\pm }\equiv \pm$(${\bar{\xi }}_1$, ${\bar{\xi }}_2$, ${\bar{\xi }}_3$)/${\bar{\xi }}_0$, where ${\bar{\xi }}_0\equiv \sqrt{{({\bar{\xi }}_1)}^2+{({\bar{\xi }}_2)}^2+{({\bar{\xi }}_3)}^2}$. The probabilities of photon emission in these states, $W_{fi}^{\pm }$, are given by Eq. (49). A stochastic procedure is defined using the fourth random number r₄: if $W_{fi}^{+}/{\bar{W}}_{fi}\ge r_4$, the polarization state ${\widehat{\mathit{\xi }}}^{+}$ will be chosen; otherwise, the polarization state will be assigned as ${\widehat{\mathit{\xi }}}^{-}$. Here, ${\bar{W}}_{fi}\equiv W_R{F}_0$ and $W_{fi}^{\pm }\equiv W_R\left(F_0+{\sum }_{j=1,3}{\xi }_j^{\pm }{F}_j\right)$.

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., p_f ≈ (1 − ω_γ/|p_i|)p_i. Besides, in this simulation, interference effects between emissions in adjacent coherent lengths (l_f ≃ λ_L/a₀) are negligible when the employed laser intensity is ultrastrong, i.e., a₀ ≫ 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.

In the context of NCS and the subsequent nonlinear Breit–Wheeler pair production, the polarization state of a photon can be characterized by the polarization unit vector $\widehat{\mathbf{P}}$, which functions as the spin component of the photon wavefunction. An arbitrary polarization $\widehat{\mathbf{P}}$ can be represented as a superposition of two orthogonal basis vectors:67$$\widehat{\mathbf{P}}=\cos ({\theta }_{\alpha }){\widehat{\mathbf{P}}}_1+\sin ({\theta }_{\alpha }){\widehat{\mathbf{P}}}_2{e}^{i{\theta }_{\beta }},$$(58)where θ_α denotes the angle between P and ${\widehat{\mathbf{P}}}_1$, while θ_β represents the absolute phase. In quantum mechanics, the photon polarization state corresponding to P can be described by the density matrix$$\rho =\frac{1}{2}\left(1+\mathit{\xi }\cdot \mathit{\sigma }\right)=\frac{1}{2}\left(\begin{array}{cc}\hfill 1+{\xi }_3\hfill & \hfill {\xi }_1-i{\xi }_2\hfill \\ \hfill {\xi }_1+i{\xi }_2\hfill & \hfill 1-{\xi }_3\hfill \end{array}\right),$$(59)where σ is the Pauli matrix, and ξ = (ξ₁, ξ₂, ξ₃) denotes the Stokes parameters, with ξ₁ = sin(2θ_α)cos(θ_β), ξ₂ = sin(2θ_α)sin(θ_β), and ξ₃ = cos(2θ_α).

Calculation of the probability of pair creation requires transformation of the Stokes parameters from the initial frame of the photon (${\widehat{\mathbf{P}}}_1$, ${\widehat{\mathbf{P}}}_2$, $\widehat{\mathbf{n}}$) to the frame of pair production (${\widehat{\mathbf{P}}}_1^{\prime }$, ${\widehat{\mathbf{P}}}_2^{\prime }$, $\widehat{\mathbf{n}}$). The vector ${\widehat{\mathbf{P}}}_1^{\prime }$ is given by $[\mathbf{E}-\widehat{\mathbf{n}}\cdot (\widehat{\mathbf{n}}\cdot \mathbf{E})+\widehat{\mathbf{n}}\times \mathbf{B}]$/$|\mathbf{E}-\widehat{\mathbf{n}}\cdot (\widehat{\mathbf{n}}\cdot \mathbf{E})+\widehat{\mathbf{n}}\times \mathbf{B}|$, and the vector ${\widehat{\mathbf{P}}}_2^{\prime }$ is obtained by taking the cross product of $\widehat{\mathbf{n}}$ and ${\widehat{\mathbf{P}}}_1^{\prime }$. Here, $\widehat{\mathbf{n}}$ represents the direction of propagation of the photon, and E and B are the electric and magnetic fields. The two groups of polarization vectors are connected via rotation through an angle ψ:$${\widehat{\mathbf{P}}}_1^{\prime }={\widehat{\mathbf{P}}}_1\cos (\psi )+{\widehat{\mathbf{P}}}_2\sin (\psi ),$$(60)$${\widehat{\mathbf{P}}}_2^{\prime }=-{\widehat{\mathbf{P}}}_1\sin (\psi )+{\widehat{\mathbf{P}}}_2\cos (\psi ).$$(61)Thus, the Stokes parameters with respect to the vectors ${\widehat{\mathbf{P}}}_1^{\prime }$, ${\widehat{\mathbf{P}}}_2^{\prime }$, and $\widehat{\mathbf{n}}$ are as follows:$$\begin{array}{c}{\xi }_1^{\prime }={\xi }_1\cos (2\psi )-{\xi }_3\sin (2\psi ),\\ {\xi }_2^{\prime }={\xi }_2,\\ {\xi }_3^{\prime }={\xi }_1\sin (2\psi )+{\xi }_3\cos (2\psi ),\end{array}$$(62)which is equivalent to a rotation:68,69$$\left(\begin{array}{c}\hfill {\xi }_1^{\prime }\hfill \\ \hfill {\xi }_2^{\prime }\hfill \\ \hfill {\xi }_3^{\prime }\hfill \end{array}\right)=\left(\begin{array}{ccc}\hfill \cos 2\psi \hfill & \hfill 0\hfill & \hfill -\sin 2\psi \hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill \sin 2\psi \hfill & \hfill 0\hfill & \hfill \cos 2\psi \hfill \end{array}\right)\left(\begin{array}{c}\hfill {\xi }_1\hfill \\ \hfill {\xi }_2\hfill \\ \hfill {\xi }_3\hfill \end{array}\right)\equiv \mathrm{R}\mathrm{O}\mathrm{T}(\psi )\cdot \mathit{\xi }.$$(63)