Simulations of spin/polarization-resolved laser–plasma interactions in the nonlinear QED regime

Feng Wan; Chong Lv; Kun Xue; Zhen-Ke Dou; Qian Zhao; Mamutjan Ababekri; Wen-Qing Wei; Zhong-Peng Li; Yong-Tao Zhao; Jian-Xing Li

doi:10.1063/5.0163929

Journals >Matter and Radiation at Extremes >Volume 8 >Issue 6 >Page 064002 > Article

Matter and Radiation at Extremes
Vol. 8, Issue 6, 064002 (2023)

Simulations of spin/polarization-resolved laser–plasma interactions in the nonlinear QED regime

Feng Wan¹, Chong Lv², Kun Xue¹, Zhen-Ke Dou¹, Qian Zhao¹, Mamutjan Ababekri¹, Wen-Qing Wei¹, Zhong-Peng Li¹, Yong-Tao Zhao¹, and Jian-Xing Li^1、a)

Author Affiliations

¹Ministry of Education Key Laboratory for Nonequilibrium Synthesis and Modulation of Condensed Matter, Shaanxi Province Key Laboratory of Quantum Information and Quantum Optoelectronic Devices, School of Physics, Xi’an Jiaotong University, Xi’an 710049, China

²Department of Nuclear Physics, China Institute of Atomic Energy, P.O. Box 275(7), Beijing 102413, China

show less

DOI: 10.1063/5.0163929 Cite this Article

Feng Wan, Chong Lv, Kun Xue, Zhen-Ke Dou, Qian Zhao, Mamutjan Ababekri, Wen-Qing Wei, Zhong-Peng Li, Yong-Tao Zhao, Jian-Xing Li. Simulations of spin/polarization-resolved laser–plasma interactions in the nonlinear QED regime[J]. Matter and Radiation at Extremes, 2023, 8(6): 064002 Copy Citation Text

EndNote(RIS)

BibTex

Plain Text

show less

Abstract

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.

I. INTRODUCTION

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_{Sch} \equiv m_{e}^{2} c^{3} / e ℏ$ , 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{align} χ & \equiv \frac{e ℏ \sqrt{{(F_{μ ν} p^{μ})}^{2}}}{m^{3} c^{4}} ~ \frac{a^{'}}{a_{Sch}}, f \equiv \frac{e^{2} ℏ^{2} F_{μ ν} F^{μ ν}}{4 m^{4} c^{6}} ~ \frac{{a_{E}}^{2} - {a_{B}}^{2}}{4 {a^{2}}_{Sch}}, \\ g \equiv \frac{e^{2} ℏ^{2} F_{μ ν} F^{μ ν *}}{4 m^{4} c^{6}} ~ \frac{{a_{E}}^{2} \cdot {a_{B}}^{2}}{4 {a^{2}}_{Sch}}, \end{align}$ 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.

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

Figure 1.Standard particle-in-cell (PIC) loop with four kernel parts.

A. Particle pushing

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 p}{d t} = \frac{q}{m} (E + β \times B),$ (1) $\frac{d x}{d t} = \frac{p}{γ},$ (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{p^{n + 1 / 2} - p^{n - 1 / 2}}{Δ t} = \frac{q}{m} (E^{n} + \frac{p^{n}}{γ^{n}} \times B^{n}),$ (3) $\frac{x^{n + 1} - x^{n}}{Δ t} = v^{n + 1 / 2},$ (4)and are solved using the standard Boris rotation:31–33 $p^{n - 1 / 2} = p^{-} - \frac{q Δ t}{2 m} E^{n},$ (5) $p^{n + 1 / 2} = p^{+} + \frac{q Δ t}{2 m} E^{n},$ (6) $p^{'} = p^{-} + p^{-} \times τ,$ (7) $p^{+} = p^{-} + p^{'} \times ς,$ (8) $τ = \frac{q Δ t}{2 m γ^{n}} B^{n},$ (9) $ς = \frac{2 τ}{1 + τ^{2}},$ (10)where $γ^{n} = \sqrt{1 + {(p^{-})}^{2}} = \sqrt{1 + {(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.

B. Field solving

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 E = ρ,$ (11) $\nabla \cdot B = 0$ (12) $\nabla \times E = - \frac{\partial B}{\partial t},$ (13) $\nabla \times B = \frac{\partial E}{\partial t} + 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.

Figure 2.(a) and (b) Yee grid and position of each field component in 3D and 2D cases, respectively. In (b), the z direction is squeezed.

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{gathered} {\frac{E_{y}^{n + 1} - E_{y}^{n}}{Δ t}|}_{i + 1 / 2} = {- \frac{B_{i + 1} - B_{i}}{Δ x}|}_{z}^{n + 1 / 2} - J_{y, i + 1 / 2}^{n + 1 / 2}, \\ {\frac{B_{z}^{n + 1 / 2} - B_{z}^{n - 1 / 2}}{Δ t}|}_{i} = {- \frac{E_{i + 1 / 2} - E_{i - 1 / 2}}{Δ x}|}_{y}^{n}; \end{gathered}$ (15)2D case (squeezing the z direction): $\begin{align} {\frac{E_{y}^{n + 1} - E_{y}^{n}}{Δ t}|}_{i + 1 / 2, j} & = {- \frac{B_{i + 1, j} - B_{i, j}}{Δ x}|}_{z}^{n + 1 / 2} - J_{y, i + 1 / 2, j}^{n + 1 / 2}, \\ {\frac{B_{z}^{n + 1 / 2} - B_{z}^{n - 1 / 2}}{Δ t}|}_{i + 1 / 2, j} & = {- \frac{E_{i + 1 / 2, j} - E_{i + 1 / 2, j}}{Δ x}|}_{y}^{n} \\ {+ \frac{E_{i + 1 / 2, j + 1 / 2} - E_{i + 1 / 2, j - 1 / 2}}{Δ y}|}_{x}^{n + 1 / 2}; \end{align}$ (16)

3D case: $\begin{align} {\frac{E_{y}^{n + 1} - E_{y}^{n}}{Δ t}|}_{i + 1 / 2, j, k + 1 / 2} & = {- \frac{B_{i + 1, j, k + 1 / 2} - B_{i, j, k + 1 / 2}}{Δ x}|}_{z}^{n + 1 / 2} \\ {+ \frac{B_{i + 1 / 2, j, k + 1} - B_{i + 1 / 2, j, k}}{Δ z}|}_{x}^{n + 1 / 2} \\ - J_{y, i + 1 / 2, j, k + 1 / 2}^{n + 1 / 2}, \\ {\frac{B_{z}^{n + 1 / 2} - B_{z}^{n - 1 / 2}}{Δ t}|}_{i, j, k + 1 / 2} & = {- \frac{E_{i + 1 / 2, j, k} - E_{i - 1 / 2, j, k}}{Δ x}|}_{y}^{n} \\ {+ \frac{E_{i + 1 / 2, j + 1, k} - E_{i + 1 / 2, j, k}}{Δ y}|}_{x}^{n} . \end{align}$ (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.

C. Current deposition

We calculate the charge current density using Esirkepov’s method, which conserves charge by satisfying Gauss’s law35 $\partial_{t} ρ + \nabla \cdot J = 0,$ (18)and removes the need for Coulomb correction.19 This algorithm computes the charge density at time steps $t - \frac{1}{2} Δ t$ and $t + \frac{1}{2} Δ t$ on each grid cell from the particle positions and velocities, i.e., $ρ_{i, j, k}^{n + 1 / 2} = \frac{1}{Δ V} \sum_{r} W (x_{r}^{n} + \frac{1}{2} v^{n + 1 / 2} Δ t) q_{r},$ (19) $ρ_{i, j, k}^{n - 1 / 2} = \frac{1}{Δ V} \sum_{r} W (x_{r}^{n} - \frac{1}{2} v^{n} Δ t) q_{r},$ (20) $δ^{n} ρ = ρ^{n + 1 / 2} - ρ^{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.

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 x_r, we find its nearest grid point ${(i, j, k)}^{g} = fl o o r (x_{r} / Δ 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.

E. Particle shape function

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_{spline} = \{\begin{cases} \frac{3}{4} - δ^{2}, & f o r j, \\ \frac{1}{2} (\frac{1}{2} \pm δ), & f o r j \pm 1, \end{cases}$ (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).

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.

Figure 3.Leapfrog algorithm for particle pushing and field advancing.

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

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{align} F_{RR,classical} \\ = \frac{2 e^{3}}{3 m c^{3}} \{γ [(\frac{\partial}{\partial t} + \frac{p}{γ m} \cdot \nabla) E + \frac{p}{γ m c} \times (\frac{\partial}{\partial t} + \frac{p}{γ m} \cdot \nabla) B] \\ + \frac{e}{m c} [E \times B + \frac{1}{γ m c} B \times (B \times p) + \frac{1}{γ m c} E (p \cdot E)] \\ - \frac{e γ}{m^{2} c^{2}} p [{(E + \frac{p}{γ m c} \times B)}^{2} - \frac{1}{γ^{2} m^{2} c^{2}} {(E \cdot p)}^{2}]\} . \end{align}$ (23)

The dimensionless form of this equation is $\begin{align} F_{RR,classical} & = \frac{2}{3} α_{f} ξ_{L} \{γ [(\frac{\partial}{\partial t} + \frac{p}{γ} \cdot \nabla) E + \frac{p}{γ} \times (\frac{\partial}{\partial t} + \frac{p}{γ} \cdot \nabla) B] \\ + [E \times B + \frac{1}{γ} B \times (B \times p) + E (p \cdot E)] \\ - γ p [{(E + \frac{p}{γ} \times B)}^{2} - \frac{1}{γ^{2}} {(E \cdot p)}^{2}]\}, \end{align}$ (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 $F_{RR,cl} ≃ \frac{2}{3} α_{f} \frac{χ_{e}^{2}}{ξ_{L}} β,$ (25)where $\begin{array}{l} χ_{e} & = \frac{e ℏ}{m^{3} c^{4}} \sqrt{| F^{μ ν} p_{ν} |^{2}} \\ \equiv ξ_{L} γ_{e} \sqrt{{(E + β \times B)}^{2} - {[β \cdot (β \cdot E)]}^{2}} \\ ≃ γ_{e} E_{⊥} ξ_{L} (1 - \cos θ) \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{align} R & \equiv | F_{RR} | / | F_{L} | \\ \sim \frac{2}{3} α_{f} γ_{e} χ_{e} ≃ 2 \times 1 0^{- 8} a_{0} γ_{e}^{2} (for wavelength 1 μ m) . \end{align}$ (26)Clearly, once $γ_{e}^{2} a_{0} ≳ 1 0^{6}$ , the radiation reaction force should be considered.

2. Modified Landau–Lifshitz (MLL) equation

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 $F_{RR,quantum} = q (χ) F_{RR,classical},$ (27)where $q (χ) = \frac{I_{QED}}{I_{C}},$ (28) $I_{QED} = m c^{2} \int c (k \cdot k^{'}) \frac{d W_{f i}}{d η d r_{0}} d r_{0},$ (29) $I_{C} = \frac{2 e^{4} E^{' 2}}{3 m^{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 (χ_{e}) \approx \frac{1}{{[1 + 4.8 (1 + χ_{e}) \ln (1 + 1.7 χ_{e}) + 2.44 χ_{e}^{2}]}^{2 / 3}},$ (31)or $q (χ_{e}) \approx \frac{1}{{(1 + 8.93 χ_{e} + 2.41 χ_{e}^{2})}^{2 / 3}} .$ (32)

Figure 4.q(χ) vs χ.

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

3. Algorithms for the radiative pusher

Here, we plug the radiative correction (either classical or quantum corrected version) into the standard Boris pusher as follows:43 $\frac{p^{n + 1 / 2} - p^{n - 1 / 2}}{Δ t} = F^{n} = F_{L}^{n} + F_{R}^{n} .$ (34)First we use the Boris step $\frac{p_{L}^{n + 1 / 2} - p_{L}^{n - 1 / 2}}{Δ t} = F_{L}^{n},$ (35)and then use the radiative correction push $\frac{p_{R}^{n + 1 / 2} - p_{R}^{n - 1 / 2}}{Δ t} = F_{R}^{n},$ (36)where $p_{L}^{n - 1 / 2} = p_{R}^{n - 1 / 2} = p^{n - 1 / 2}$ , and the final momentum is given by $p^{n + 1 / 2} = p_{L}^{n + 1 / 2} + p_{R}^{n + 1 / 2} - p^{n - 1 / 2} = p_{L}^{n + 1 / 2} + F_{R}^{n} Δ 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 $p (τ) = p_{0} - A (τ) + \hat{k} \frac{A^{2} (τ) - 2 p_{0} \cdot A (τ)}{2 (γ_{0} - p_{0} \cdot \hat{k})}$ (38) $γ (τ) = γ_{0} + \frac{A^{2} (τ) - 2 p_{0} \cdot A (τ)}{2 (γ_{0} - p_{0} \cdot \hat{k})}$ (39)where $A (τ) = - \int_{τ_{0}}^{τ} E (τ^{'}) d τ^{'}$ is the external field vector potential, τ is the proper time, $\hat{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.

Figure 5.Dynamics of an electron [p₀ = (4000, 0, 0)] scattering with an ultraintense linearly polarized laser pulse of $E_{y} = 100 \exp [- {(\frac{ϕ - 100}{10 π})}^{2}] \cos ϕ$ , with ϕ ≡ t + x. Here, Lo., LL., and MLL. denote results calculated from the Lorentz, LL, and modified LL equations, respectively.

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

The nonradiative spin dynamics of an electron are given by $\begin{align} {(\frac{d S}{d t})}_{T} & = S \times Ω \\ \equiv S \times [- (\frac{g}{2} - 1) \frac{γ_{e}}{γ_{e} + 1} (β \cdot B) \cdot β \\ + (\frac{g}{2} - 1 + \frac{1}{γ_{e}}) B - (\frac{g}{2} - \frac{γ_{e}}{γ_{e} + 1}) β \times E], \end{align}$ (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.

2. Radiative T-BMT equation

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 S}{d t} = {(\frac{d S}{d t})}_{T} + {(\frac{d S}{d t})}_{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 ${(\frac{d S}{d t})}_{R} = - P [ψ_{1} (χ) S + ψ_{2} (χ) (S \cdot β) β + ψ_{3} (χ) {\hat{n}}_{B}] .$ (42)Here, $\begin{array}{l} P = \frac{α_{f}}{\sqrt{3} π γ_{e} ξ_{L}}, ψ_{1} (χ_{e}) = \int_{0}^{\infty} u^{″} d u K_{2 / 3} (u^{'}), \\ ψ_{2} (χ_{e}) = \int_{0}^{\infty} u^{″} d u \int_{u^{'}}^{\infty} d x K_{1 / 3} (x) - ψ_{1} (χ_{e}), \\ ψ_{3} (χ) = \int_{0}^{\infty} u^{″} d u K_{1 / 3} (u^{'}), u^{'} = \frac{2 u}{3 χ_{e}}, u^{″} = \frac{u^{2}}{{(1 + u)}^{3}}, \end{array}$ where K_n is the nth-order modified Bessel function of the second kind, ${\hat{n}}_{B} = β \times \hat{a}$ , and β and $\hat{a}$ denote the normalized velocity and acceleration vectors, respectively.58,59

3. Algorithms for simulating spin precession

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: $S^{'} = S^{n - 1 / 2} + S^{n - 1 / 2} \times t,$ (43) $S^{n + 1 / 2} = S^{n - 1 / 2} + S^{'} \times o,$ (44) $t = \frac{q Δ t}{2} Ω^{n},$ (45) $o = \frac{2 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 $S_{T}^{n - 1 / 2} = S^{n - 1 / 2} + \frac{Δ t}{2} {(\frac{d S}{d t})}_{R},$ (47)Boris T-BMT Eqs. (43)–(46), $S^{n + 1 / 2} = S_{T}^{n + 1 / 2} + \frac{Δ t}{2} {(\frac{d S}{d t})}_{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).

Figure 6.Spin dynamics of an electron [p₀ = (4000, 0, 0), s₀ = (1, 0, 0)] scattering with an ultraintense linearly polarized laser pulse of $E_{y} = 100 \exp [- {(\frac{ϕ - 100}{10 π})}^{2}] \cos ϕ$ , with ϕ ≡ t + x. Here, A, B, C, and D denote results calculated using the Lorentz + T-BMT, Lorentz + radiative T-BMT, LL + radiative T-BMT, and MLL + radiative T-BMT equations, respectively.

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

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_{f i}}{d u d t} = \frac{W_{R}}{2} (F_{0} + ξ_{1} F_{1} + ξ_{2} F_{2} + ξ_{3} F_{3}),$ (49)where the photon polarization is represented by the Stokes parameters (ξ₁, ξ₂, ξ₃), defined with respect to the axes ${\hat{P}}_{1} = \hat{a} - \hat{n} (\hat{n} \cdot \hat{a})$ and ${\hat{P}}_{2} = \hat{n} \times {\hat{P}}_{1}$ ,63 with the photon emission direction $\hat{n} = p_{e} / | p_{e} |$ along the momentum p_e of the ultrarelativistic electron. The variables introduced in Eq. (49) are as follows: $\begin{align} F_{0} & = - {(2 + u)}^{2} [{I n t K}_{1 / 3} (u^{'}) - 2 K_{2 / 3} (u^{'})] (1 + S_{i f}) \\ + u^{2} (1 - S_{i f}) [{I n t K}_{1 / 3} (u^{'}) + 2 K_{2 / 3} (u^{'})] + 2 u^{2} S_{i f} {I n t K}_{1 / 3} (u^{'}) \\ - (4 u + 2 u^{2}) (S_{f} + S_{i}) \cdot [\hat{n} \times \hat{a}] K_{1 / 3} (u^{'}) - 2 u^{2} (S_{f} - S_{i}) \\ \cdot [\hat{n} \times \hat{a}] K_{1 / 3} (u^{'}) - 4 u^{2} [{I n t K}_{1 / 3} (u^{'}) - K_{2 / 3} (u^{'})] \\ \times (S_{i} \cdot \hat{n}) (S_{f} \cdot \hat{n}), \end{align}$ (50) $\begin{align} F_{1} & = - 2 u^{2} {I n t K}_{1 / 3} (u^{'}) \{(S_{i} \cdot \hat{a}) S_{f} \cdot [\hat{n} \times \hat{a}] + (S_{f} \cdot \hat{a}) S_{i} \cdot [\hat{n} \times \hat{a}]\} \\ + 4 u [(S_{i} \cdot \hat{a}) (1 + u) + (S_{f} \cdot \hat{a})] K_{1 / 3} (u^{'}) \\ + 2 u (2 + u) \hat{n} \cdot [S_{f} \times S_{i}] K_{2 / 3} (u^{'}), \end{align}$ (51) $\begin{align} F_{2} & = - \{2 u^{2} \{(S_{i} \cdot \hat{n}) S_{f} \cdot [\hat{n} \times \hat{a}] + (S_{f} \cdot \hat{n}) S_{i} \cdot [\hat{n} \times \hat{a}]\} \\ + 2 u (2 + u) \hat{a} \cdot [S_{f} \times S_{i}]\} K_{1 / 3} (u^{'}) - 4 u [(S_{i} \cdot \hat{n}) \\ + (S_{f} \cdot \hat{n}) (1 + u)] {I n t K}_{1 / 3} (u^{'}) + 4 u (2 + u) \\ \times [(S_{i} \cdot \hat{n}) + (S_{f} \cdot \hat{n})] K_{2 / 3} (u^{'}), \end{align}$ (52) $\begin{align} F_{3} & = 4 [1 + u + (1 + u + \frac{1}{2} u^{2}) S_{i f} - \frac{1}{2} u^{2} (S_{i} \cdot \hat{n}) (S_{f} \cdot \hat{n})] K_{2 / 3} (u^{'}) \\ + 2 u^{2} \{S_{i} \cdot [\hat{n} \times \hat{a}] S_{f} \cdot [\hat{n} \times \hat{a}] - (S_{i} \cdot \hat{a}) (S_{f} \cdot \hat{a})\} {I n t K}_{1 / 3} (u^{'}) \\ - 4 u \{(1 + u) S_{i} [\hat{n} \times \hat{a}] + S_{f} [\hat{n} \times \hat{a}]\} K_{1 / 3} (u^{'}), \end{align}$ (53)where $W_{R} = \frac{α_{f}}{8 \sqrt{3} π ξ_{L} {(1 + u)}^{3}}, u^{'} = \frac{2 u}{3 χ}, u = \frac{ω_{γ}}{ε_{i} - ω_{γ}},$ ω_γ is the emitted photon energy, ɛ_i is the electron energy before radiation, $\hat{a} = a / | 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: ${I n t K}_{1 / 3} (u^{'}) \equiv \int_{u^{'}}^{\infty} d z K_{1 / 3} (z) .$

By summing over the photon polarizations, the electron spin-resolved emission probability can be written as12,15,64 $\begin{align} \frac{d^{2} W_{f i}}{d u d t} & = W_{R} \{- {(2 + u)}^{2} [{I n t K}_{1 / 3} (u^{'}) - 2 K_{2 / 3} (u^{'})] (1 + S_{i f}) \\ + u^{2} [{I n t K}_{1 / 3} (u^{'}) + 2 K_{2 / 3} (u^{'})] (1 - S_{i f}) \\ + 2 u^{2} S_{i f} {I n t K}_{1 / 3} (u^{'}) - (4 u + 2 u^{2}) (S_{f} + S_{i}) [n \times \hat{a}] \\ \times K_{1 / 3} (u^{'}) - 2 u^{2} (S_{f} - S_{i}) [n \times \hat{a}] K_{1 / 3} (u^{'}) \\ - 4 u^{2} [{I n t K}_{1 / 3} (u^{'}) - K_{2 / 3} (u^{'})] (S_{i} \cdot n) (S_{f} \cdot n)\}, \end{align}$ (54)and by summing over the final states S_f, the initial spin-resolved radiation probability is obtained: $\begin{align} \frac{d^{2} {\bar{W}}_{f i}}{d u d t} & = 8 W_{R} \{- (1 + u) {I n t K}_{1 / 3} (u^{'}) \\ + (2 + 2 u + u^{2}) K_{2 / 3} (u^{'}) - u S_{i} \cdot [n \times \hat{a}] K_{1 / 3} (u^{'})\} . \end{align}$ (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 $S_{f}^{R} = \frac{g}{w + f \cdot 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_{rad}}{d u d t} = W_{r} (w + f \cdot S_{i} + g \cdot S_{f}),$ (57)where $\begin{array}{l} w = - (1 + u) K_{1 / 3} (ρ^{'}) + (2 + 2 u + u^{2}) K_{2 / 3} (ρ^{'}), \\ f = u {I n t K}_{1 / 3} (ρ^{'}) \hat{v} \times \hat{a}, \\ g = - (1 + u) [K_{1 / 3} (ρ^{'}) - 2 K_{2 / 3} (ρ^{'})] S_{i} \\ - (1 + u) u {I n t K}_{1 / 3} (ρ^{'}) \hat{v} \times \hat{a} \\ - u^{2} [K_{1 / 3} (ρ^{'}) - K_{2 / 3} (ρ^{'})] (S_{i} \cdot \hat{v}) \hat{v} . \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}}_{f i} (χ_{e}, r_{1}, γ_{e}, S_{i}) / d u d t \cdot Δ 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_{flip} \equiv W_{f i}^{↑} / P$ and $W_{f i}^{↓} / P$ , respectively, where $W_{f i}^{↑, ↓} \equiv d^{2} W_{f i}^{↑, ↓} / d u d t \cdot Δ 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{S}}_{f} = \frac{S_{i} (1 - W Δ t) - f Δ t}{1 - (W + f \cdot S_{i}) Δ t},$ where $W \equiv 16 W_{R} [- (1 + u) {I n t K}_{1 / 3} (u^{'}) + (2 + 2 u + u^{2}) K_{2 / 3} (u^{'})]$ and $f \equiv - 16 W_{R} n \times \hat{a} K_{1 / 3} (u^{'})$ .12,66 Then, the SQA is given by ${\bar{S}}_{f} / | {\bar{S}}_{f} |$ , and the probability for the aligned case is given by $| {\bar{S}}_{f} |$ and that for the antiparallel case by $1 - | {\bar{S}}_{f} |$ .

Figure 7.Flowchart of spin- and polarization-resolved 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 ${\hat{ξ}}^{\pm} \equiv \pm$ ( ${\bar{ξ}}_{1}$ , ${\bar{ξ}}_{2}$ , ${\bar{ξ}}_{3}$ )/ ${\bar{ξ}}_{0}$ , where ${\bar{ξ}}_{0} \equiv \sqrt{{({\bar{ξ}}_{1})}^{2} + {({\bar{ξ}}_{2})}^{2} + {({\bar{ξ}}_{3})}^{2}}$ . The probabilities of photon emission in these states, $W_{f i}^{\pm}$ , are given by Eq. (49). A stochastic procedure is defined using the fourth random number r₄: if $W_{f i}^{+} / {\bar{W}}_{f i} \geq r_{4}$ , the polarization state ${\hat{ξ}}^{+}$ will be chosen; otherwise, the polarization state will be assigned as ${\hat{ξ}}^{-}$ . Here, ${\bar{W}}_{f i} \equiv W_{R} F_{0}$ and $W_{f i}^{\pm} \equiv W_{R} (F_{0} + \sum_{j = 1, 3} ξ_{j}^{\pm} F_{j})$ .

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.

Figure 8.Dynamics of 1000 electrons via stochastic NCS, with the simulation parameters the same as those in Fig. 6. Blue lines are for ten sampled electrons, and black ones are the average value over 1000 sample particles.

Figure 9.Dynamics of an electron beam (particle number N_e = 10⁴), with colors denoting the number density in arbitrary units and a logarithmic scale (a.u.); other parameters are the same as those in Fig. 6.

Figure 10.(a) Energy spectra of scattered electrons (black curve) and generated photons (red curve). (b) Energy-dependent Stokes parameters ${\bar{ξ}}_{2}$ and $\bar{ξ_{3}}$ , i.e., circular and linear polarization with respect to the y and z axes. The simulation parameters are the same as those in Fig. 6.

2. Definition and transformation of Stokes parameters

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 $\hat{P}$ , which functions as the spin component of the photon wavefunction. An arbitrary polarization $\hat{P}$ can be represented as a superposition of two orthogonal basis vectors:67 $\hat{P} = \cos (θ_{α}) {\hat{P}}_{1} + \sin (θ_{α}) {\hat{P}}_{2} e^{i θ_{β}},$ (58)where θ_α denotes the angle between P and ${\hat{P}}_{1}$ , while θ_β represents the absolute phase. In quantum mechanics, the photon polarization state corresponding to P can be described by the density matrix $ρ = \frac{1}{2} (1 + ξ \cdot σ) = \frac{1}{2} (\begin{matrix} 1 + ξ_{3} & ξ_{1} - i ξ_{2} \\ ξ_{1} + i ξ_{2} & 1 - ξ_{3} \end{matrix}),$ (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 ( ${\hat{P}}_{1}$ , ${\hat{P}}_{2}$ , $\hat{n}$ ) to the frame of pair production ( ${\hat{P}}_{1}^{'}$ , ${\hat{P}}_{2}^{'}$ , $\hat{n}$ ). The vector ${\hat{P}}_{1}^{'}$ is given by $[E - \hat{n} \cdot (\hat{n} \cdot E) + \hat{n} \times B]$ / $| E - \hat{n} \cdot (\hat{n} \cdot E) + \hat{n} \times B |$ , and the vector ${\hat{P}}_{2}^{'}$ is obtained by taking the cross product of $\hat{n}$ and ${\hat{P}}_{1}^{'}$ . Here, $\hat{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 ψ: ${\hat{P}}_{1}^{'} = {\hat{P}}_{1} \cos (ψ) + {\hat{P}}_{2} \sin (ψ),$ (60) ${\hat{P}}_{2}^{'} = - {\hat{P}}_{1} \sin (ψ) + {\hat{P}}_{2} \cos (ψ) .$ (61)Thus, the Stokes parameters with respect to the vectors ${\hat{P}}_{1}^{'}$ , ${\hat{P}}_{2}^{'}$ , and $\hat{n}$ are as follows: $\begin{gathered} ξ_{1}^{'} = ξ_{1} \cos (2 ψ) - ξ_{3} \sin (2 ψ), \\ ξ_{2}^{'} = ξ_{2}, \\ ξ_{3}^{'} = ξ_{1} \sin (2 ψ) + ξ_{3} \cos (2 ψ), \end{gathered}$ (62)which is equivalent to a rotation:68,69 $(\begin{matrix} ξ_{1}^{'} \\ ξ_{2}^{'} \\ ξ_{3}^{'} \end{matrix}) = (\begin{matrix} \cos 2 ψ & 0 & - \sin 2 ψ \\ 0 & 1 & 0 \\ \sin 2 ψ & 0 & \cos 2 ψ \end{matrix}) (\begin{matrix} ξ_{1} \\ ξ_{2} \\ ξ_{3} \end{matrix}) \equiv R O T (ψ) \cdot ξ .$ (63)

D. Nonlinear Breit–Wheeler (NBW) pair production

When the energy of a photon exceeds the rest mass of an electron–positron pair, i.e., ω_γ ≥ 2m_ec², and the photon is subjected to an ultraintense field a₀ ≫ 1, the related nonlinear quantum parameter χ_γ can reach unity. Here, $χ_{γ} \equiv (e ℏ^{2} / m^{3} c^{4}) \sqrt{| F^{μ ν} k_{ν} |^{2}}$ and is approximately equal to $2 a_{0} ω_{γ} ξ_{L}$ 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

The polarization-resolved NBW probability rate with dependence on the positron energy is given by $\frac{d^{2} W_{pair}^{\pm}}{d ε_{+} d t} = \frac{1}{2} (G_{0} + ξ_{1} G_{1} + ξ_{2} G_{2} + ξ_{3} G_{3}),$ (64)where the polarization-independent term G₀ and polarization-related terms G_1,2,3 are given by $\begin{align} G_{0} & = \frac{W_{0}}{2} \{{I n t K}_{1 / 3} (ρ) + \frac{ε_{-}^{2} + ε_{+}^{2}}{ε_{-} ε_{+}} K_{2 / 3} (ρ) \\ + [{I n t K}_{1 / 3} (ρ) - 2 K_{2 / 3} (ρ)] (S_{-} \cdot S_{+}) \\ + K_{1 / 3} (ρ) [- \frac{ε_{γ}}{ε_{+}} (S_{+} \cdot {\hat{b}}_{+}) + \frac{ε_{γ}}{ε_{-}} (S_{-} \cdot {\hat{b}}_{+})] \\ + [\frac{ε_{-}^{2} + ε_{+}^{2}}{ε_{-} ε_{+}} {I n t K}_{1 / 3} (ρ) - \frac{{(ε_{+} - ε_{-})}^{2}}{ε_{-} ε_{+}} K_{2 / 3} (ρ)] \\ \times (S_{+} \cdot {\hat{v}}_{+}) (S_{-} \cdot {\hat{v}}_{+})\}, \end{align}$ (65) $\begin{align} G_{1} & = \frac{W_{0}}{2} \{K_{1 / 3} (ρ) [- \frac{ε_{γ}}{ε_{-}} (S_{+} \cdot {\hat{a}}_{+}) + \frac{ε_{γ}}{ε_{+}} (S_{-} \cdot {\hat{a}}_{+})] \\ + \frac{ε_{+}^{2} - ε_{-}^{2}}{2 ε_{-} ε_{+}} K_{2 / 3} (ρ) (S_{-} \times S_{+}) \cdot {\hat{v}}_{+} - \frac{ε_{γ}^{2}}{2 ε_{-} ε_{+}} {I n t K}_{1 / 3} (ρ) \\ \times [(S_{+} \cdot \hat{a}) (S_{-} \cdot \hat{b}) + (S_{-} \cdot {\hat{a}}_{+}) (S_{+} \cdot {\hat{b}}_{+})]\}, \end{align}$ (66) $\begin{align} G_{2} & = \frac{W_{0}}{2} \{\frac{ε_{γ}^{2}}{2 ε_{-} ε_{+}} K_{1 / 3} (ρ) (S_{-} \times S_{+}) \cdot {\hat{a}}_{+} + \frac{ε_{+}^{2} - ε_{-}^{2}}{2 ε_{-} ε_{+}} K_{1 / 3} (ρ) \\ \times [(S_{-} \cdot {\hat{v}}_{+}) (S_{+} \cdot {\hat{b}}_{+}) + (S_{+} \cdot {\hat{v}}_{+}) (S_{-} \cdot {\hat{b}}_{+})] \\ + [\frac{ε_{γ}}{ε_{-}} {I n t K}_{1 / 3} (ρ) - \frac{ε_{+}^{2} - ε_{-}^{2}}{ε_{-} ε_{+}} K_{2 / 3} (ρ)] (S_{-} \cdot {\hat{v}}_{+}) \\ + [\frac{ε_{γ}}{ε_{+}} {I n t K}_{1 / 3} (ρ) + \frac{ε_{+}^{2} - ε_{-}^{2}}{ε_{-} ε_{+}} K_{2 / 3} (ρ)] (S_{+} \cdot {\hat{v}}_{+})\}, \end{align}$ (67) $\begin{align} G_{3} & = \frac{W_{0}}{2} \{- K_{2 / 3} (ρ) + \frac{ε_{-}^{2} + ε_{+}^{2}}{2 ε_{-} ε_{+}} K_{2 / 3} (ρ) (S_{-} \cdot S_{+}) - K_{1 / 3} (ρ) \\ \times [\frac{ε_{γ}}{ε_{+}} (S_{-} \cdot {\hat{b}}_{+}) - \frac{ε_{γ}}{ε_{-}} (S_{+} \cdot {\hat{b}}_{+})] + \frac{ε_{γ}^{2}}{2 ε_{-} ε_{+}} {I n t K}_{1 / 3} (ρ) \\ \times [(S_{+} \cdot {\hat{b}}_{+}) (S_{-} \cdot {\hat{b}}_{+}) - (S_{+} \cdot {\hat{a}}_{+}) (S_{-} \cdot {\hat{a}}_{+})] \\ - \frac{{(ε_{+} - ε_{-})}^{2}}{2 ε_{-} ε_{+}} K_{2 / 3} (ρ) (S_{+} \cdot {\hat{v}}_{+}) (S_{-} \cdot {\hat{v}}_{+})\}, \end{align}$ (68)where $\begin{array}{l} W_{0} = \frac{α}{\sqrt{3} π ω_{γ}^{' 2}}, ω_{γ}^{'} = \frac{ε_{γ}}{m_{e} c^{2}}, \\ ρ = \frac{2 ε_{γ}^{2}}{3 χ_{γ} ε_{-} ε_{+}} = \frac{2}{3 δ (1 - δ)}, δ = \frac{ε_{+}}{ε_{γ}}, \end{array}$ α is the fine structure constant, ɛ_γ, ɛ₋, and ɛ₊ are the energies of the parent photon and the created electron and positron, respectively, ${\hat{v}}_{+} = v_{+} / | v_{+} |$ (with v₊ the positron velocity), ${\hat{a}}_{+} = a_{+} / | a_{+} |$ (with a₊ the positron acceleration in the rest frame of the positron), ${\hat{b}}_{+} = v_{+} \times a_{+} / | v_{+} \times a_{+} |$ , ξ₁, ξ₂, and ξ₃ are the Stokes parameters of the γ-photon, and S₊ and S₋ are the positron and electron spin vectors, respectively. K_n is again the nth-order modified Bessel function of the second kind, and the function Int K_1/3 is defined after Eq. (53). Note that the Stokes parameters must be transformed from the photon initial frame $({\hat{P}}_{1}, {\hat{P}}_{2}, \hat{n})$ to the pair production frame $({\hat{P}}_{1}^{'}, {\hat{P}}_{2}^{'}, \hat{n})$ ; see the transformations of the Stokes parameters in Sec. III C 2.

By summing over the electron spin, the pair production probability depending on the positron spin S₊ and the photon polarization ξ is obtained as $\begin{align} \frac{d^{2} W_{pair}^{+}}{d ε_{+} d t} & = W_{0} \{{I n t K}_{1 / 3} (ρ) + \frac{ε_{-}^{2} + ε_{+}^{2}}{ε_{-} ε_{+}} K_{2 / 3} (ρ) \\ - \frac{ε_{γ}}{ε_{+}} K_{1 / 3} (ρ) (S_{+} \cdot {\hat{b}}_{+}) - ξ_{1} [\frac{ε_{γ}}{ε_{-}} K_{1 / 3} (ρ) (S_{+} \cdot {\hat{a}}_{+})] \\ + ξ_{2} [\frac{ε_{+}^{2} - ε_{-}^{2}}{ε_{-} ε_{+}} K_{2 / 3} (ρ) + \frac{ε_{γ}}{ε_{+}} {I n t K}_{1 / 3} (ρ)] (S_{+} \cdot {\hat{v}}_{+}) \\ - ξ_{3} [K_{2 / 3} (ρ) - \frac{ε_{γ}}{ε_{-}} K_{1 / 3} (ρ) (S_{+} \cdot {\hat{b}}_{+})]\} . \end{align}$ (69)This can be rewritten as $\frac{d^{2} W_{pair}^{+}}{d ε_{+} d t} = W_{0} (C + S_{+} \cdot D),$ (70)where $C = {I n t K}_{1 / 3} (ρ) + \frac{ε_{-}^{2} + ε_{+}^{2}}{ε_{-} ε_{+}} K_{2 / 3} (ρ) - ξ_{3} K_{2 / 3} (ρ),$ (71) $\begin{align} D & = - (\frac{ε_{γ}}{ε_{+}} - ξ_{3} \frac{ε_{γ}}{ε_{-}}) K_{1 / 3} (ρ) {\hat{b}}_{+} - ξ_{1} \frac{ε_{γ}}{ε_{-}} K_{1 / 3} (ρ) {\hat{a}}_{+} \\ + ξ_{2} [\frac{ε_{+}^{2} - ε_{-}^{2}}{ε_{-} ε_{+}} K_{2 / 3} (ρ) + \frac{ε_{γ}}{ε_{+}} {I n t K}_{1 / 3} (ρ)] {\hat{v}}_{+} . \end{align}$ (72)When a photon decays to a pair, the positron spin state is instantaneously collapsed into one of its basis states defined by the instantaneous SQA, along the energy-resolved average polarization $S_{+}^{(ε_{+})} = D / C$ .

Similarly, by summing over the positron spin, the pair production probability depending on the electron spin S₋ and the photon polarization is obtained as $\frac{d^{2} W_{pair}^{-}}{d ε_{+} d t} = W_{0} (C + S_{-} \cdot D^{'}),$ (73) $\begin{align} D^{'} & = (\frac{ε_{γ}}{ε_{-}} - ξ_{3} \frac{ε_{γ}}{ε_{+}}) K_{1 / 3} (ρ) {\hat{b}}_{+} + ξ_{1} \frac{ε_{γ}}{ε_{+}} K_{1 / 3} (ρ) {\hat{a}}_{+} \\ - ξ_{2} [\frac{ε_{+}^{2} - ε_{-}^{2}}{ε_{-} ε_{+}} K_{2 / 3} (ρ) - \frac{ε_{γ}}{ε_{-}} {I n t K}_{1 / 3} (ρ)] {\hat{v}}_{+} . \end{align}$ (74)

The pair production probability, depending solely on the photon polarization, is determined by summing over both positron and electron spins: $\frac{d^{2} W_{pair}}{d ε_{+} d t} = 2 W_{0} \{{I n t K}_{1 / 3} (ρ) + \frac{ε_{-}^{2} + ε_{+}^{2}}{ε_{-} ε_{+}} K_{2 / 3} (ρ) - ξ_{3} K_{2 / 3} (ρ)\} .$ (75)

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 ɛ₊ = r₁ɛ_γ when the probability density P ≡ d²W_pair/dɛ₊dt · Δt of pair production is greater than or equal to a random number r₂ within the range [0, 1). Here, d²W_pair/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 P_1,2,3,4, each representing spin parallel or antiparallel to the SQA, where P_1,2,3,4 is computed from Eq. (64). Finally, a random number r₃ 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.

Figure 11.Flowchart of the spin- and polarization-resolved nonlinear Breit–Wheeler (NBW) pair production process.

Figure 12.(a) Normalized energy spectrum (black solid curve) and energy-resolved longitudinal spin polarization (red solid curve) of positrons. (b) Statistics of the longitudinal spin components of generated positrons. The laser and electron beam parameters are consistent with those in Fig. 9.

E. High-energy bremsstrahlung

High-energy bremsstrahlung is another important emission mechanism, and it can also be modeled using an MC collision model.73 The MC collision model was tested using the Geant4 code,74 and the results are presented here. The bremsstrahlung emission is described by the cross-section from Ref. 75: $\begin{align} \frac{d σ_{e Z}}{d ω} (ω, y) \\ = \frac{α r_{0}^{2}}{ω} \{(\frac{4}{3} - \frac{4}{3} y + y^{2}) [Z^{2} (ϕ_{1} - \frac{4}{3} \ln Z - 4 f) + Z (ψ_{1} - \frac{8}{3} \ln Z)] \\ + \frac{2}{3} (1 - y) [Z^{2} (ϕ_{1} - ϕ_{2}) + Z (ψ_{1} - ψ_{2})]\}, \end{align}$ (76)where y = ℏω/E_e is the ratio of the energy of the emitted photon to that of the incident electron, r₀ is the classical electron radius, the functions ϕ_1,2 and ψ_1,2 depend on the screening potential by atomic electrons, and f is the Coulomb correction term. When the atomic number of the target is greater than 5, we use Eqs. (3.38)–(3.41) from Ref. 75 to calculate these functions. However, for targets with Z < 5, the approximate screening functions are unsuitable and require modification.

The PENELOPE code76 utilizes another method, which involves tabulated data from Ref. 77. This method transforms the “scaled” bremsstrahlung differential cross-section (DCS) to a differential cross-section as follows:76 $\frac{d σ_{br}}{d ω} = \frac{Z^{2}}{β^{2}} \frac{1}{ω} χ (Z, E_{e}, y),$ (77)where β = v/c is the normalized electron velocity. Integrating this expression over the photon frequencies yields a tabulated total cross-section σ_br(E_e, y) for MC simulation, i.e., the direct sampling method can be used.

The electron and positron DCS are related by $\frac{d σ_{br}^{+}}{d ω} = F_{p} (Z, E_{e}) \frac{d σ_{br}^{-}}{d ω},$ (78)where F_p(Z, E_e) is an analytical approximation factor that can be found in Ref. 76. A high level of accuracy was demonstrated in Ref. 76, with a difference of only ∼0.5% compared with the results reported in Ref. 78.

The bremsstrahlung implementation is based on direct MC sampling. Given an incident electron with energy E_e and velocity v, the probability of triggering a bremsstrahlung event is calculated as P_br = 1 − e^Δ/, where Δs = vΔt, v = |v| is the incident particle velocity, Δt is the time interval, λ = 1/nσ(E_e), n is the target particle density, and σ(E_e) is the total cross-section. A random number r₁ is then generated and compared with P_br. If r₁ < P_br, then a bremsstrahlung event is triggered. The energy of the resulting photon is determined by generating another random number r₂, which is then multiplied by σ_br(E_e) to obtain the energy ratio y through σ(y, E_e) = σ(E_e)r₂. Finally, a photon with energy ℏω = E_ey 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 10⁵ primaries, colliding with a 5 mm Au target with Z = 79 and ρ = 19.3 g/cm³ and a 5 mm Al target with Z = 13 and ρ = 2.7 g/cm³. 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.

Figure 13.Bremsstrahlung of 100 MeV electrons: (a) scattered electron spectra; (b) yield photon spectra. Solid curves represent PIC results and dashed curves Geant4 results. Reproduced with permission from F. Wan et al., Eur. Phys. J. D 71, 236 (2017). Copyright 2017, EDP Sciences, SIF, Springer-Verlag GmbH Germany.

Figure 14.Bremsstrahlung of 1 GeV electrons: (a) scattered electron spectra; (b) yield photon spectra. Solid curves represent PIC results and dashed curves Geant4 results. Reproduced with permission from F. Wan et al., Eur. Phys. J. D 71, 236 (2017). Copyright 2017, EDP Sciences, SIF, Springer-Verlag GmbH Germany.

F. Vacuum birefringence

In addition to the NBW processes, another important process for polarized photons in ultraintense laser–matter interactions is vacuum birefringence (VB). In this paper, we utilize Eq. (4.26) from Ref. 79 to calculate the refractive index n for a photon with arbitrary energy ω (wavelength λ) in a constant weak EM field [|E|(|B|) ≪ E_cr]. We include the electric field and assume relativistic units c = ℏ = 1. The resulting expression is $\begin{align} n & \approx 1 - \frac{α χ_{γ}^{2} m^{2}}{16 π ω^{2}} \int_{- 1}^{1} d υ (1 - υ^{2}) \{\begin{matrix} \frac{1}{2} (1 + \frac{1}{3} υ^{2}) \\ 1 - \frac{1}{3} υ^{2} \end{matrix}\} \\ \times [π x^{4 / 3} {G i}^{'} (x^{2 / 3}) - i \frac{x^{2}}{\sqrt{3}} K_{2 / 3} (\frac{2}{3} x)], \end{align}$ (79)where α is the fine structure constant, m is the electron mass, χ_γ is the nonlinear quantum parameter as defined earlier, x = 4/[(1 − υ²)χ_γ], and Gi′(x) is the derivative of the Scorer function. $E_{red, ⊥} = E_{⊥} + \hat{k} \times B_{⊥}$ is the transverse reduced field (acceleration field for electrons). The first and second rows in the ${\begin{matrix} \dots \\ \dots \end{matrix}}$ correspond to the eigenmodes parallel and perpendicular to the reduced field, respectively. After extraction of a factor $D = \frac{α}{90 π} {(\frac{e | E_{red ⊥} |}{m^{2}})}^{2} \equiv \frac{α}{90 π} \frac{χ_{γ}^{2}}{ω^{2} / m^{2}} .$ and separation into real and imaginary parts, Eq. (79) becomes $\begin{align} R e (n) & = 1 - \frac{45}{4} D \int_{0}^{1} d υ (1 - υ^{2}) \{\begin{matrix} \frac{1}{2} (1 + \frac{1}{3} υ^{2}) \\ 1 - \frac{1}{3} υ^{2} \end{matrix}\} \\ \times [π x^{4 / 3} {G i}^{'} (x^{2 / 3})], \end{align}$ (80) $I m (n) = \frac{45}{4} D \int_{0}^{1} d υ (1 - υ^{2}) \{\begin{matrix} \frac{1}{2} (1 + \frac{1}{3} υ^{2}) \\ 1 - \frac{1}{3} υ^{2} \end{matrix}\} [\frac{x^{2}}{\sqrt{3}} K_{2 / 3} (\frac{2}{3} x)] .$ (81)

In the weak-field limit of χ_γ ≪ 1, the imaginary part associated with pair production is negligible. We now define $M (χ_{γ}) = - \frac{45}{4} \int_{0}^{1} d υ (1 - υ^{2}) \{\begin{matrix} \frac{1}{2} (1 + \frac{1}{3} υ^{2}) \\ 1 - \frac{1}{3} υ^{2} \end{matrix}\} π x^{4 / 3} {G i}^{'} (x^{2 / 3}),$ (82)yielding $R e (n) = 1 + M (χ_{γ}) D \equiv 1 + M (χ_{γ}) \frac{α}{90 π} \frac{χ_{γ}^{2}}{ω^{2} / m^{2}} .$ (83)The numerical results for M(χ_γ) and comparisons with the low-energy-limit (ω_γ ≪ m) constants are given in Fig. 15.

Figure 15.(a) M(χ_γ) (red and blue solid curves) and the corresponding low-energy-limit constants, with red and blue dash-dotted lines equal to 4 and 7, respectively. (b) Relative error between M(χ_γ) and the low-energy-limit constant.

In the limit of χ_γ ≪ 1, the real part simplifies to $R e (n) = 1 + D \{\begin{matrix} 4_{+} \\ 7_{-} \end{matrix}\}$ (84)and can be used to simulate the VB effect with good accuracy for χ_γ ≪ 1. Note that these results are identical to those in Refs. 79–81. For large χ_γ, two interpolated refractive indices are used.

The phase retardation between two orthogonal components is given by $δ ϕ = ϕ_{+} - ϕ_{-} = Δ n 2 π l / λ = - 3 D 2 π l / λ$ , where l is the propagation length, and the VB effect is equivalent to a rotation of the Stokes parameters: $(\begin{matrix} ξ_{1}^{'} \\ ξ_{2}^{'} \\ ξ_{3}^{'} \end{matrix}) = (\begin{matrix} \cos δ ϕ & - \sin δ ϕ & 0 \\ \sin δ ϕ & \cos δ ϕ & 0 \\ 0 & 0 & 1 \end{matrix}) (\begin{matrix} ξ_{1} \\ ξ_{2} \\ ξ_{3} \end{matrix}) \equiv Q E D (δ ϕ) \cdot ξ .$ (85)

The VB effect of the probe photons in the PIC code is simulated with Algorithm 1.82ALGORITHM 1.

VB effect in SLIPs.

For an example of the VB effect, see Fig. 16.

Figure 16.VB effect of a γ-photon [ɛ_γ = 1 GeV, ξ = (1, 0, 0)] propagating through (a) static crossed fields with E_y = −B_z = 100 and (b) a laser field (the same as in Fig. 6).

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.

Figure 17.Data structure of SLIPs.

Figure 18.Framework of SLIPs.

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: $L o r e n t z and T - B M T \Rightarrow R a d i a t i v e c o r r e c t i o n \Rightarrow V B \Rightarrow N B W \Rightarrow N C S$ and Bremss.

V. POLARIZED PARTICLE SIMULATIONS

In this section, we present known results that were calculated from the single-particle mode using SLIPs. The spin-resolved NCS/NBW are evaluated by generating spin-polarized electrons/positrons. The simulation setups used in this study are identical to those described in Refs. 10 and 64.

A. Polarized electron/positron simulation

To simulate the generation of spin-polarized electrons, we utilized an elliptically polarized laser with an intensity a₀ = 30, a wavelength λ₀ = 1 μm, and an ellipticity a_y,0/a_x,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.

$Generation of polarized electrons: (a) number density log_{10}(d2N/dθxdθy) (a.u.); (b) spin polarization Sx.$

Figure 19.Generation of polarized electrons: (a) number density log_{10}(d²N/dθ_xdθ_y) (a.u.); (b) spin polarization S_x.

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.

Figure 20.Generation of LP γ-photons: (a) number density log₁₀(d²N/dθ_xdθ_y) (a.u.); (b) linear polarization ξ₃.

Figure 21.Generation of CP γ-photons with longitudinally polarized electrons: (a) number density log₁₀(d²N/dθ_xdθ_y) (a.u.); (b) circular polarization |ξ₂|.

C. Laser–plasma interactions

Finally, we present a simulation result demonstrating the interaction between an ultraintense laser with a normalized intensity a₀ = 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 $e_{1} \equiv β \times \dot{β}$ and $e_{2} \equiv \hat{n} \times e_{1}$ , where $\hat{n}$ 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% $(\bar{ξ_{3}} ≃ 0.6)$ , 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.

Figure 22.Laser–plasma interaction via 2D simulation: (a)–(c) spatial distributions of E_x, E_y, and B_z, respectively; (d)–(f) number densities (in logrithm) of target electrons, generated NBW positrons, and NCS γ-photons, respectively.

Figure 23.Photons generated by laser–plasma interaction: (a) number density with respect to energy and angle, i.e., log₁₀(dN²/dγ_γdθ) (a.u.) with $γ_{γ} {\equiv E}_{γ} / m_{e} c^{2}$ and θ ≡ p_y/p_x; (b) energy- and angle-resolved linear polarization degree ${\bar{ξ}}_{3}$ ; (c) energy-resolved number and polarization distributions; (d) angle-resolved number and polarization distributions.

Figure 24.Positrons generated by laser–plasma interaction: (a) number density with respect to energy and angle, i.e., dN²/dγ₊dθ (a.u.), with γ₊ ≡ ɛ₊/m_ec² and θ ≡ arctan(p_y/p_x); (b) energy- and angle-resolved spin component ${\bar{S}}_{z}$ ; (c) normalized angular distribution n(θ) ≡ dN/dθ (a.u.); (d) angular distribution of ${\bar{S}}_{z}$ (i.e., energy-averaged); (e) normalized energy distribution n(ɛ₊) ≡ dN/dɛ₊(a.u.).

VI. OUTLOOK

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.

ACKNOWLEDGMENTS

Acknowledgment. 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).

References

[1] A.Di Piazza, Piazza A.Di, C.Müller, Z.Hatsagortsyan K., C. and, C.Müller, Piazza A.Di, C.Müller, Z.Hatsagortsyan K., C. and, K. Z.Hatsagortsyan, Piazza A.Di, C.Müller, Z.Hatsagortsyan K., C. and, C. H.Keitel. Extremely high-intensity laser interactions with fundamental quantum systems. Rev. Mod. Phys., 84, 1177-1228(2012).

[2] A. R.Bell, R.Bell A., J. G.Kirk. Possibility of prolific pair production with high-power lasers. Phys. Rev. Lett., 101, 200403(2008).

[4] E.Esarey, E.Esarey, B.Schroeder C., W. and, C. B.Schroeder, E.Esarey, B.Schroeder C., W. and, W. P.Leemans. Physics of laser-driven plasma-based electron accelerators. Rev. Mod. Phys., 81, 1229-1285(2009).

[5] V. I.Ritus. Quantum effects of the interaction of elementary particles with an intense electromagnetic field. J. Russ. Laser Res., 6, 497-617(1985).

[6] V. N.Baier, N.Baier V., M.Katkov V., V. and, V. M.Katkov, N.Baier V., M.Katkov V., V. and, V. M.Strakhovenko. Electromagnetic Processes at High Energies in Oriented Single Crystals(1998).

[7] J. D.Jackson. Classical Electrodynamics(2021).

[8] I. M.Ternov, M.Ternov I., A. A.Sokolov. Radiation from Relativistic Electrons(1986).

[9] D.Del Sorbo, Sorbo D.Del, D.Seipt, G.Blackburn T., G. A., D.Murphy C., G.Kirk J., C. and, D.Seipt, Sorbo D.Del, D.Seipt, G.Blackburn T., G. A., D.Murphy C., G.Kirk J., C. and, T. G.Blackburn, Sorbo D.Del, D.Seipt, G.Blackburn T., G. A., D.Murphy C., G.Kirk J., C. and, A. G. R.Thomas, Sorbo D.Del, D.Seipt, G.Blackburn T., G. A., D.Murphy C., G.Kirk J., C. and, C. D.Murphy, Sorbo D.Del, D.Seipt, G.Blackburn T., G. A., D.Murphy C., G.Kirk J., C. and, J. G.Kirk, Sorbo D.Del, D.Seipt, G.Blackburn T., G. A., D.Murphy C., G.Kirk J., C. and, C. P.Ridgers. Spin polarization of electrons by ultraintense lasers. Phys. Rev. A, 96, 043407(2017).

[10] Y.-F.Li, Y.-F.Li, R.Shaisultanov, Z.Hatsagortsyan K., F.Wan, H.Keitel C., J.-X.Li and, R.Shaisultanov, Y.-F.Li, R.Shaisultanov, Z.Hatsagortsyan K., F.Wan, H.Keitel C., J.-X.Li and, K. Z.Hatsagortsyan, Y.-F.Li, R.Shaisultanov, Z.Hatsagortsyan K., F.Wan, H.Keitel C., J.-X.Li and, F.Wan, Y.-F.Li, R.Shaisultanov, Z.Hatsagortsyan K., F.Wan, H.Keitel C., J.-X.Li and, C. H.Keitel, Y.-F.Li, R.Shaisultanov, Z.Hatsagortsyan K., F.Wan, H.Keitel C., J.-X.Li and, J.-X.Li. Ultrarelativistic electron-beam polarization in single-shot interaction with an ultraintense laser pulse. Phys. Rev. Lett., 122, 154801(2019).

[11] B.King, and B.King, S.Tang. Nonlinear Compton scattering of polarized photons in plane-wave backgrounds. Phys. Rev. A, 102, 022809(2020).

[12] Y.-F.Li, Y.-F.Li, R.Shaisultanov, Y.-Y.Chen, F.Wan, Z.Hatsagortsyan K., H.Keitel C., J.-X.Li and, R.Shaisultanov, Y.-F.Li, R.Shaisultanov, Y.-Y.Chen, F.Wan, Z.Hatsagortsyan K., H.Keitel C., J.-X.Li and, Y.-Y.Chen, Y.-F.Li, R.Shaisultanov, Y.-Y.Chen, F.Wan, Z.Hatsagortsyan K., H.Keitel C., J.-X.Li and, F.Wan, Y.-F.Li, R.Shaisultanov, Y.-Y.Chen, F.Wan, Z.Hatsagortsyan K., H.Keitel C., J.-X.Li and, K. Z.Hatsagortsyan, Y.-F.Li, R.Shaisultanov, Y.-Y.Chen, F.Wan, Z.Hatsagortsyan K., H.Keitel C., J.-X.Li and, C. H.Keitel, Y.-F.Li, R.Shaisultanov, Y.-Y.Chen, F.Wan, Z.Hatsagortsyan K., H.Keitel C., J.-X.Li and, J.-X.Li. Polarized ultrashort brilliant multi-GeV γ-rays via single-shot laser–electron interaction. Phys. Rev. Lett., 124, 014801(2020).

[13] Z.Gong, Z.Gong, Z.Hatsagortsyan K., C. and, K. Z.Hatsagortsyan, Z.Gong, Z.Hatsagortsyan K., C. and, C. H.Keitel. Retrieving transient magnetic fields of ultrarelativistic laser plasma via ejected electron polarization. Phys. Rev. Lett., 127, 165002(2021).

[16] S.Tang, S.Tang, B.King, H.Hu and, B.King, S.Tang, B.King, H.Hu and, H.Hu. Highly polarised gamma photons from electron-laser collisions. Phys. Lett. B, 809, 135701(2020).

[17] H.-H.Song, H.-H.Song, W.-M.Wang, Y.-T.Li and, W.-M.Wang, H.-H.Song, W.-M.Wang, Y.-T.Li and, Y.-T.Li. Dense polarized positrons from laser-irradiated foil targets in the QED regime. Phys. Rev. Lett., 129, 035001(2022).

[19] C.Birdsall, and C.Birdsall, A.Langdon. Plasma Physics via Computer Simulation(2018).

[21] R. A.Fonseca, P. M. A.Sloot, A.Fonseca R., M. P., O.Silva L., G.Hoekstra A., J. C., S.Tsung F., K.Decyk V., J. and, A. G.Hoekstra, W.Lu, C.Ren, M. P., G.Hoekstra A., B.Mori W., S.Deng, J. C., J. and, S.Lee, T.Katsouleas, C. J. K.Tan, M. P., J. and, L. O.Silva, G.Hoekstra A., A.Fonseca R., J. C., J. and, O.Silva L., J. J.Dongarra, S.Tsung F., K.Decyk V., W.Lu, C.Ren, B.Mori W., S.Deng, S.Lee, T.Katsouleas, J. and, F. S.Tsung, A.Fonseca R., O.Silva L., S.Tsung F., K.Decyk V., W.Lu, C.Ren, B.Mori W., S.Deng, S.Lee, T.Katsouleas, J. and, V. K.Decyk, A.Fonseca R., O.Silva L., S.Tsung F., K.Decyk V., W.Lu, C.Ren, B.Mori W., S.Deng, S.Lee, T.Katsouleas, J. and, W.Lu, A.Fonseca R., O.Silva L., S.Tsung F., K.Decyk V., W.Lu, C.Ren, B.Mori W., S.Deng, S.Lee, T.Katsouleas, J. and, C.Ren, A.Fonseca R., O.Silva L., S.Tsung F., K.Decyk V., W.Lu, C.Ren, B.Mori W., S.Deng, S.Lee, T.Katsouleas, J. and, W. B.Mori, A.Fonseca R., O.Silva L., S.Tsung F., K.Decyk V., W.Lu, C.Ren, B.Mori W., S.Deng, S.Lee, T.Katsouleas, J. and, S.Deng, A.Fonseca R., O.Silva L., S.Tsung F., K.Decyk V., W.Lu, C.Ren, B.Mori W., S.Deng, S.Lee, T.Katsouleas, J. and, S.Lee, A.Fonseca R., O.Silva L., S.Tsung F., K.Decyk V., W.Lu, C.Ren, B.Mori W., S.Deng, S.Lee, T.Katsouleas, J. and, T.Katsouleas, A.Fonseca R., O.Silva L., S.Tsung F., K.Decyk V., W.Lu, C.Ren, B.Mori W., S.Deng, S.Lee, T.Katsouleas, J. and, J. C.Adam. OSIRIS: A three-dimensional, fully relativistic particle in cell code for modeling plasma based accelerators. Computational Science—ICCS 2002, 342-351(2002).

[24] D.Wu, D.Wu, W.Yu, S.Fritzsche, X. and, W.Yu, D.Wu, W.Yu, S.Fritzsche, X. and, S.Fritzsche, D.Wu, W.Yu, S.Fritzsche, X. and, X. T.He. Particle-in-cell simulation method for macroscopic degenerate plasmas. Phys. Rev. E, 102, 033312(2020).

[25] Y.-F.Li, Y.-F.Li, Y.-Y.Chen, Z.Hatsagortsyan K., C. and, Y.-Y.Chen, Y.-F.Li, Y.-Y.Chen, Z.Hatsagortsyan K., C. and, K. Z.Hatsagortsyan, Y.-F.Li, Y.-Y.Chen, Z.Hatsagortsyan K., C. and, C. H.Keitel. Helicity transfer in strong laser fields via the electron anomalous magnetic moment. Phys. Rev. Lett., 128, 174801(2022).

[26] C.Danson, C.Danson, D.Hillier, N.Hopps, D.Neely and, D.Hillier, C.Danson, D.Hillier, N.Hopps, D.Neely and, N.Hopps, C.Danson, D.Hillier, N.Hopps, D.Neely and, D.Neely. Petawatt class lasers worldwide. High Power Laser Sci. Eng., 3, e3(2015).

[27]

[29]

[30] Z.Gan, K.Yamanouchi, Z.Gan, K.Yamanouchi, L.Yu, K.Midorikawa, C.Wang, L.Roso and, Y.Liu, K.Midorikawa, Y.Xu, K.Yamanouchi, K.Midorikawa, W.Li, L.Roso and, S.Li, L.Roso, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, L.Yu, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, C.Wang, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, Y.Liu, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, Y.Xu, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, W.Li, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, S.Li, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, L.Yu, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, X.Wang, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, X.Liu, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, J.Chen, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, Y.Peng, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, L.Xu, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, B.Yao, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, X.Zhang, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, L.Chen, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, Y.Tang, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, X.Wang, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, D.Yin, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, X.Liang, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, Y.Leng, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, R.Li, Z.Gan, L.Yu, C.Wang, Y.Liu, Y.Xu, W.Li, S.Li, L.Yu, X.Wang, X.Liu, J.Chen, Y.Peng, L.Xu, B.Yao, X.Zhang, L.Chen, Y.Tang, X.Wang, D.Yin, X.Liang, Y.Leng, R.Li, Z.Xu and, Z.Xu. The Shanghai superintense ultrafast laser facility (SULF) project. Progress in Ultrafast Intense Laser Science XVI, 199-217(2021).

[31] O.Buneman. Time-reversible difference procedures. J. Comput. Phys., 1, 517-535(1967).

[32] J. P.Boris, P.Boris J., R. A.Shanny(1970).

[33] H.Qin, H.Qin, S.Zhang, J.Xiao, J.Liu, Y.Sun, W. and, S.Zhang, H.Qin, S.Zhang, J.Xiao, J.Liu, Y.Sun, W. and, J.Xiao, H.Qin, S.Zhang, J.Xiao, J.Liu, Y.Sun, W. and, J.Liu, H.Qin, S.Zhang, J.Xiao, J.Liu, Y.Sun, W. and, Y.Sun, H.Qin, S.Zhang, J.Xiao, J.Liu, Y.Sun, W. and, W. M.Tang. Why is Boris algorithm so good?. Phys. Plasmas, 20, 084503(2013).

[34] Y.Kane. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media. IEEE Trans. Antennas Propag., 14, 302-307(1966).

[35] T.Esirkepov. Exact charge conservation scheme for particle-in-cell simulation with an arbitrary form-factor. Comput. Phys. Commun., 135, 144-153(2001).

[36] P. A. M.Dirac. Classical theory of radiating electrons. Proc. R. Soc. London, Ser. A, 167, 148-169(1938).

[37] L. D.Landau, D.Landau L., E.Lifshitz. The Classical Theory of Fields(1999).

[38] R.Ekman, R.Ekman, T.Heinzl, A.Ilderton and, T.Heinzl, R.Ekman, T.Heinzl, A.Ilderton and, A.Ilderton. Reduction of order, resummation, and radiation reaction. Phys. Rev. D, 104, 036002(2021).

[39] R.Ekman. Reduction of order and transseries structure of radiation reaction. Phys. Rev. D, 105, 056016(2022).

[40] A.Ilderton, and A.Ilderton, G.Torgrimsson. Radiation reaction in strong field QED. Phys. Lett. B, 725, 481-486(2013).

[41] D.Seipt, and D.Seipt, A. G. R.Thomas. Kinetic theory for spin-polarized relativistic plasmas(2023).

[42] N.Neitz, and N.Neitz, A.Di Piazza. Electron-beam dynamics in a strong laser field including quantum radiation reaction. Phys. Rev. A, 90, 022102(2014).

[43] M.Tamburini, M.Tamburini, F.Pegoraro, Piazza A.Di, H.Keitel C., A.Macchi and, F.Pegoraro, M.Tamburini, F.Pegoraro, Piazza A.Di, H.Keitel C., A.Macchi and, A.Di Piazza, M.Tamburini, F.Pegoraro, Piazza A.Di, H.Keitel C., A.Macchi and, C. H.Keitel, M.Tamburini, F.Pegoraro, Piazza A.Di, H.Keitel C., A.Macchi and, A.Macchi. Radiation reaction effects on radiation pressure acceleration. New J. Phys., 12, 123005(2010).

[44] J.Hein, S. V.Bulanov, J.Hein, V.Bulanov S., G.Korn, Z.Esirkepov T., L. and, M.Kando, G.Korn, K.Koga J., T.Nakamura, J.Hein, S.Bulanov S., G.Korn, L. and, G.Zhidkov A., L. O.Silva, Y.Kato, G.Korn and, T. Z.Esirkepov, V.Bulanov S., Z.Esirkepov T., M.Kando, K.Koga J., T.Nakamura, S.Bulanov S., G.Zhidkov A., Y.Kato, G.Korn and, M.Kando, V.Bulanov S., Z.Esirkepov T., M.Kando, K.Koga J., T.Nakamura, S.Bulanov S., G.Zhidkov A., Y.Kato, G.Korn and, J. K.Koga, V.Bulanov S., Z.Esirkepov T., M.Kando, K.Koga J., T.Nakamura, S.Bulanov S., G.Zhidkov A., Y.Kato, G.Korn and, T.Nakamura, V.Bulanov S., Z.Esirkepov T., M.Kando, K.Koga J., T.Nakamura, S.Bulanov S., G.Zhidkov A., Y.Kato, G.Korn and, S. S.Bulanov, V.Bulanov S., Z.Esirkepov T., M.Kando, K.Koga J., T.Nakamura, S.Bulanov S., G.Zhidkov A., Y.Kato, G.Korn and, A. G.Zhidkov, V.Bulanov S., Z.Esirkepov T., M.Kando, K.Koga J., T.Nakamura, S.Bulanov S., G.Zhidkov A., Y.Kato, G.Korn and, Y.Kato, V.Bulanov S., Z.Esirkepov T., M.Kando, K.Koga J., T.Nakamura, S.Bulanov S., G.Zhidkov A., Y.Kato, G.Korn and, G.Korn. On extreme field limits in high power laser matter interactions: Radiation dominant regimes in high intensity electromagnetic wave interaction with electrons. High-Power, High-Energy, and High-Intensity Laser Technology; and Research Using Extreme Light: Entering New Frontiers with Petawatt-Class Lasers(2013).

[45] A. I.Nikishov, I.Nikishov A., V. I.Ritus. Quantum processes in the field of a plane electromagnetic wave and in a constant field. Sov. Phys. JETP, 19, 529-541(1964).

[46] I. V.Sokolov, V.Sokolov I., M.Naumova N., A.Nees J., A.Mourou G., V. and, N. M.Naumova, V.Sokolov I., M.Naumova N., A.Nees J., A.Mourou G., V. and, J. A.Nees, V.Sokolov I., M.Naumova N., A.Nees J., A.Mourou G., V. and, G. A.Mourou, V.Sokolov I., M.Naumova N., A.Nees J., A.Mourou G., V. and, V. P.Yanovsky. Dynamics of emitting electrons in strong laser fields. Phys. Plasmas, 16, 093115(2009).

[47] A.Di Piazza, Piazza A.Di, Z.Hatsagortsyan K., C. and, K. Z.Hatsagortsyan, Piazza A.Di, Z.Hatsagortsyan K., C. and, C. H.Keitel. Quantum radiation reaction effects in multiphoton Compton scattering. Phys. Rev. Lett., 105, 220403(2010).

[48] A. G. R.Thomas, G. A., P.Ridgers C., S.Bulanov S., J.Griffin B., S. and, C. P.Ridgers, G. A., P.Ridgers C., S.Bulanov S., J.Griffin B., S. and, S. S.Bulanov, G. A., P.Ridgers C., S.Bulanov S., J.Griffin B., S. and, B. J.Griffin, G. A., P.Ridgers C., S.Bulanov S., J.Griffin B., S. and, S. P. D.Mangles. Strong radiation-damping effects in a gamma-ray source generated by the interaction of a high-intensity laser with a wakefield-accelerated electron beam. Phys. Rev. X, 2, 041004(2012).

[49] I. V.Sokolov, V.Sokolov I., A.Nees J., P.Yanovsky V., M.Naumova N., G. and, J. A.Nees, V.Sokolov I., A.Nees J., P.Yanovsky V., M.Naumova N., G. and, V. P.Yanovsky, V.Sokolov I., A.Nees J., P.Yanovsky V., M.Naumova N., G. and, N. M.Naumova, V.Sokolov I., A.Nees J., P.Yanovsky V., M.Naumova N., G. and, G. A.Mourou. Emission and its back-reaction accompanying electron motion in relativistically strong and QED-strong pulsed laser fields. Phys. Rev. E, 81, 036412(2010).

[50] F.Niel, F.Niel, C.Riconda, F.Amiranoff, R.Duclous, M.Grech and, C.Riconda, F.Niel, C.Riconda, F.Amiranoff, R.Duclous, M.Grech and, F.Amiranoff, F.Niel, C.Riconda, F.Amiranoff, R.Duclous, M.Grech and, R.Duclous, F.Niel, C.Riconda, F.Amiranoff, R.Duclous, M.Grech and, M.Grech. From quantum to classical modeling of radiation reaction: A focus on stochasticity effects. Phys. Rev. E, 97, 043209(2018).

[51] H. K.Avetissian. Relativistic Nonlinear Electrodynamics: The QED Vacuum and Matter in Super-Strong Radiation Fields, 88(2015).

[52] A. D.Piazza. Exact solution of the Landau–Lifshitz equation in a plane wave. Lett. Math. Phys., 83, 305-313(2008).

[53] S. X.Hu, X.Hu S., C. H.Keitel. Spin signatures in intense laser–ion interaction. Phys. Rev. Lett., 83, 4709-4712(1999).

[54] M. W.Walser, W.Walser M., J.Urbach D., Z.Hatsagortsyan K., X.Hu S., C. and, D. J.Urbach, W.Walser M., J.Urbach D., Z.Hatsagortsyan K., X.Hu S., C. and, K. Z.Hatsagortsyan, W.Walser M., J.Urbach D., Z.Hatsagortsyan K., X.Hu S., C. and, S. X.Hu, W.Walser M., J.Urbach D., Z.Hatsagortsyan K., X.Hu S., C. and, C. H.Keitel. Spin and radiation in intense laser fields. Phys. Rev. A, 65, 043410(2002).

[55] G. R.Mocken, R.Mocken G., C. H.Keitel. FFT-split-operator code for solving the Dirac equation in 2 + 1 dimensions. Comput. Phys. Commun., 178, 868-882(2008).

[56] H.Bauke, H.Bauke, S.Ahrens, H.Keitel C., R.Grobe and, S.Ahrens, H.Bauke, S.Ahrens, H.Keitel C., R.Grobe and, C. H.Keitel, H.Bauke, S.Ahrens, H.Keitel C., R.Grobe and, R.Grobe. Relativistic spin operators in various electromagnetic environments. Phys. Rev. A, 89, 052101(2014).

[57] M.Wen, M.Wen, H.Keitel C., H.Bauke and, C. H.Keitel, M.Wen, H.Keitel C., H.Bauke and, H.Bauke. Spin-one-half particles in strong electromagnetic fields: Spin effects and radiation reaction. Phys. Rev. A, 95, 042102(2017).

[58] V. N.Ba?er. Radiative polarization of electron in storage rings. Sov. Phys.-Usp., 14, 695-714(1972).

[59] R.-T.Guo, R.-T.Guo, Y.Wang, R.Shaisultanov, F.Wan, Z.-F.Xu, Y.-Y.Chen, Z.Hatsagortsyan K., J.-X.Li and, Y.Wang, R.-T.Guo, Y.Wang, R.Shaisultanov, F.Wan, Z.-F.Xu, Y.-Y.Chen, Z.Hatsagortsyan K., J.-X.Li and, R.Shaisultanov, R.-T.Guo, Y.Wang, R.Shaisultanov, F.Wan, Z.-F.Xu, Y.-Y.Chen, Z.Hatsagortsyan K., J.-X.Li and, F.Wan, R.-T.Guo, Y.Wang, R.Shaisultanov, F.Wan, Z.-F.Xu, Y.-Y.Chen, Z.Hatsagortsyan K., J.-X.Li and, Z.-F.Xu, R.-T.Guo, Y.Wang, R.Shaisultanov, F.Wan, Z.-F.Xu, Y.-Y.Chen, Z.Hatsagortsyan K., J.-X.Li and, Y.-Y.Chen, R.-T.Guo, Y.Wang, R.Shaisultanov, F.Wan, Z.-F.Xu, Y.-Y.Chen, Z.Hatsagortsyan K., J.-X.Li and, K. Z.Hatsagortsyan, R.-T.Guo, Y.Wang, R.Shaisultanov, F.Wan, Z.-F.Xu, Y.-Y.Chen, Z.Hatsagortsyan K., J.-X.Li and, J.-X.Li. Stochasticity in radiative polarization of ultrarelativistic electrons in an ultrastrong laser pulse. Phys. Rev. Res., 2, 033483(2020).

[60] J. G.Kirk, G.Kirk J., R.Bell A., I.Arka and, A. R.Bell, G.Kirk J., R.Bell A., I.Arka and, I.Arka. Pair production in counter-propagating laser beams. Plasma Phys. Controlled Fusion, 51, 085008(2009).

[61] C.Ridgers, C.Ridgers, J.Kirk, R.Duclous, T.Blackburn, C.Brady, K.Bennett, T.Arber, A.Bell and, J.Kirk, C.Ridgers, J.Kirk, R.Duclous, T.Blackburn, C.Brady, K.Bennett, T.Arber, A.Bell and, R.Duclous, C.Ridgers, J.Kirk, R.Duclous, T.Blackburn, C.Brady, K.Bennett, T.Arber, A.Bell and, T.Blackburn, C.Ridgers, J.Kirk, R.Duclous, T.Blackburn, C.Brady, K.Bennett, T.Arber, A.Bell and, C.Brady, C.Ridgers, J.Kirk, R.Duclous, T.Blackburn, C.Brady, K.Bennett, T.Arber, A.Bell and, K.Bennett, C.Ridgers, J.Kirk, R.Duclous, T.Blackburn, C.Brady, K.Bennett, T.Arber, A.Bell and, T.Arber, C.Ridgers, J.Kirk, R.Duclous, T.Blackburn, C.Brady, K.Bennett, T.Arber, A.Bell and, A.Bell. Modelling gamma-ray photon emission and pair production in high-intensity laser–matter interactions. J. Comput. Phys., 260, 273-285(2014).

[63] D.Green, and D.Green, C.Harvey. SIMLA: Simulating particle dynamics in intense laser and other electromagnetic fields via classical and quantum electrodynamics. Comput. Phys. Commun., 192, 313-321(2015).

[64] F.Wan, F.Wan, R.Shaisultanov, Y.-F.Li, Z.Hatsagortsyan K., H.Keitel C., J.-X.Li and, R.Shaisultanov, F.Wan, R.Shaisultanov, Y.-F.Li, Z.Hatsagortsyan K., H.Keitel C., J.-X.Li and, Y.-F.Li, F.Wan, R.Shaisultanov, Y.-F.Li, Z.Hatsagortsyan K., H.Keitel C., J.-X.Li and, K. Z.Hatsagortsyan, F.Wan, R.Shaisultanov, Y.-F.Li, Z.Hatsagortsyan K., H.Keitel C., J.-X.Li and, C. H.Keitel, F.Wan, R.Shaisultanov, Y.-F.Li, Z.Hatsagortsyan K., H.Keitel C., J.-X.Li and, J.-X.Li. Ultrarelativistic polarized positron jets via collision of electron and ultraintense laser beams. Phys. Lett. B, 800, 135120(2020).

[65] V. I., and V.I.Ritus, A. I.. Quantum Electrodynamics of Phenomena in a Strong Field. Trudy Fiz. Inst. Akad. Nauk SSSR, 111(1979).

[66] K.Yokoyaet?al.. Cain, version 2.42(2011).

[67] V. B.Berestetskii, B.Berestetskii V., M.Lifshitz E., L. and, E. M.Lifshitz, B.Berestetskii V., M.Lifshitz E., L. and, L. P.Pitaevskii. Quantum Electrodynamics: Volume 4, 4(1982).

[68] T. N.Wistisen, N.Wistisen T., U. I.Uggerh?j. Vacuum birefringence by Compton backscattering through a strong field. Phys. Rev. D, 88, 053009(2013).

[69] S.Bragin, S.Bragin, S.Meuren, H.Keitel C., A.Di and, S.Meuren, S.Bragin, S.Meuren, H.Keitel C., A.Di and, C. H.Keitel, S.Bragin, S.Meuren, H.Keitel C., A.Di and, A.Di Piazza. High-energy vacuum birefringence and dichroism in an ultrastrong laser field. Phys. Rev. Lett., 119, 250403(2017).

[70] Y.-Y.Chen, Y.-Y.Chen, P.-L.He, R.Shaisultanov, Z.Hatsagortsyan K., C. and, P.-L.He, Y.-Y.Chen, P.-L.He, R.Shaisultanov, Z.Hatsagortsyan K., C. and, R.Shaisultanov, Y.-Y.Chen, P.-L.He, R.Shaisultanov, Z.Hatsagortsyan K., C. and, K. Z.Hatsagortsyan, Y.-Y.Chen, P.-L.He, R.Shaisultanov, Z.Hatsagortsyan K., C. and, C. H.Keitel. Polarized positron beams via intense two-color laser pulses. Phys. Rev. Lett., 123, 174801(2019).

[71] Y.-Y.Chen, Y.-Y.Chen, Z.Hatsagortsyan K., H.Keitel C., R.Shaisultanov and, K. Z.Hatsagortsyan, Y.-Y.Chen, Z.Hatsagortsyan K., H.Keitel C., R.Shaisultanov and, C. H.Keitel, Y.-Y.Chen, Z.Hatsagortsyan K., H.Keitel C., R.Shaisultanov and, R.Shaisultanov. Electron spin- and photon polarization-resolved probabilities of strong-field QED processes. Phys. Rev. D, 105, 116013(2022).

[73] F.Wan, F.Wan, C.Lv, M.Jia, H.Sang, B.Xie and, C.Lv, F.Wan, C.Lv, M.Jia, H.Sang, B.Xie and, M.Jia, F.Wan, C.Lv, M.Jia, H.Sang, B.Xie and, H.Sang, F.Wan, C.Lv, M.Jia, H.Sang, B.Xie and, B.Xie. Photon emission by bremsstrahlung and nonlinear Compton scattering in the interaction of ultraintense laser with plasmas. Eur. Phys. J. D, 71, 236(2017).

[74] S.Agostinelli, S.Agostinelli, J.Allison, K.Amako, J.Apostolakiset?al., J.Allison, S.Agostinelli, J.Allison, K.Amako, J.Apostolakiset?al., K.Amako, S.Agostinelli, J.Allison, K.Amako, J.Apostolakiset?al., J.Apostolakiset?al.. Geant4—A simulation toolkit. Nucl. Instrum. Methods Phys. Res., Sect. A, 506, 250-303(2003).

[75] Y.-S.Tsai. Pair production and bremsstrahlung of charged leptons. Rev. Mod. Phys., 46, 815-851(1974).

[76] NEA, “PENELOPE. A code system for Monte Carlo simulation of electron and photon transport(2018).

[77] S. M.Seltzer, M.Seltzer S., M. J.Berger. Bremsstrahlung energy spectra from electrons with kinetic energy 1 keV–10 GeV incident on screened nuclei and orbital electrons of neutral atoms with Z = 1–100. At. Data Nucl. Data Tables, 35, 345-418(1986).

[78] L.Kim, L.Kim, H.Pratt R., M.Seltzer S., M. and, R. H.Pratt, L.Kim, H.Pratt R., M.Seltzer S., M. and, S. M.Seltzer, L.Kim, H.Pratt R., M.Seltzer S., M. and, M. J.Berger. Ratio of positron to electron bremsstrahlung energy loss: An approximate scaling law. Phys. Rev. A, 33, 3002-3009(1986).

[79] G. M.Shore. Superluminality and UV completion. Nucl. Phys. B, 778, 219-258(2007).

[80] F.Karbstein. Photon polarization tensor in a homogeneous magnetic or electric field. Phys. Rev. D, 88, 085033(2013).