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

*I*

_{0}reaches or exceeds 10

^{22}W/cm

^{2}.

^{1,2}For example, when the laser intensity is of the order of 10

^{21}–10

^{22}W/cm

^{2}, i.e., the normalized peak laser field strength parameter

*a*

_{0}≡

*eE*

_{0}/

*m*

_{e}

*cω*

_{0}∼ 10, electrons can be accelerated to GeV energies

^{3,4}(with Lorentz factor

*γ*

_{e}∼ 10

^{3}or higher) in a centimeter-long gas plasma, where −

*e*and

*m*

_{e}are the charge and mass of the electron,

*E*

_{0}and

*ω*

_{0}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

*ω*

_{0}= 2

*πc*/

*λ*

_{0}and that the wavelength of the laser is

*λ*

_{0}= 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*

_{0}∼ 10

^{4}–10

^{5}. Such a field strength is close to the QED critical field strength (Schwinger critical field strength) $ESch\u2261me2c3/e\u210f$, i.e.,

*a*

_{Sch}=

*m*

_{e}

*c*

^{2}/

*ℏω*

_{0}≃ 4.1 × 10

^{5}, 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

*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*

_{0}≪

*a*

_{Sch}) conditions,

*f*,

*g*≪

*χ*

^{2}, 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., *d***s**/*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 effect^{8} (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*_{e}*c*, velocity *v* with *c*, energy *ɛ* with *m*_{e}*c*^{2}, EM fields *E* and *B* with *m*_{e}*cω*/*e*, force *F* with *m*_{e}*cω*, charge *q* with *e*, charge density *ρ* with *k*^{3}*e*, and current density *J* with *k*^{3}*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}

### A. Particle pushing

**p**≡

*γm*

**v**,

**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

^{31–33}

### B. Field solving

The standard finite-difference method in the time domain for the Maxwell equations is to discretize field variables on a spatial grid and advance forward in time. Here, following the well-known Yee-grid approach,^{34} we put **E** and **B** as in Fig. 2(a), which automatically satisfies the two curl equations. For lower-dimensional simulations, extra dimensions are squeezed, as shown in the 2D example in Fig. 2(b). In these dimension-reduced simulations, field components in the disappeared dimensions can be seen as uniform, i.e., the gradient is 0.

Using Esirkepov’s method of current deposition,^{35} the current is calculated from the charge density via charge conservation, i.e., ∂_{t}*ρ* + ∇ · **J** = 0. Once the initial condition obeys Gauss’s law, ∇ · **E** = *ρ*, this law is automatically embedded. This can be verified by taking the gradient of Eq. (14): 0 = ∇ · (∇ × **B**) = ∂_{t}(∇ · **E**) + ∇ · **J** = ∂_{t}(∇ · **E** − *ρ*), i.e., the temporal variation in the violation of Gauss’s law is 0. Therefore, in the field solver, only the two curl equations are solved. We take the *E*_{y} and *B*_{z} components as examples here:

*y*and

*z*directions):

*z*direction):

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

^{35}

^{,}

^{19}This algorithm computes the charge density at time steps $t\u221212\Delta t$ and $t+12\Delta t$ on each grid cell from the particle positions and velocities, i.e.,

*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=floor(xr/\Delta x+12)$ and its nearest half grid point (*i*, *j*, *k*)^{h} = 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

*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}

*δ*= (

*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.

## 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

*χ*

_{e}≪ 10

^{−2}),

^{37}and, in Gaussian units, takes the form

*α*

_{f}=

*e*

^{2}/

*cℏ*is the fine structure constant and

*ξ*

_{L}=

*ℏω*/

*m*

_{e}

*c*

^{2}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

*θ*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}:

#### 2. Modified Landau–Lifshitz (MLL) equation

*m*

_{e}

*c*

^{2}.

^{44}Once

*χ*

_{e}is larger than 10

^{−2}, 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 correction

^{45–48}

^{,}

*W*

_{fi}being the radiation probability,

^{49},

*η*=

*k*

_{0}

*z*−

*ω*

_{0}

*t*,

*r*

_{0}= 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:

^{23,50}

*χ*≳ 10

^{−2}, the quantum-corrected version should be used.

#### 3. Algorithms for the radiative pusher

^{43}

^{51},

*τ*is the proper time, $k\u0302$ is the normalized wavevector of the field, and

*γ*,

**p**, and

*γ*

_{0},

**p**

_{0}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.

### 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 solver^{55} 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

**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

_{n}is the

*n*th-order modified Bessel function of the second kind, $n\u0302B=\beta \xd7a\u0302$, and

**and $a\u0302$ denote the normalized velocity and acceleration vectors, respectively.**

*β*^{58,59}

#### 3. Algorithms for simulating spin precession

**Ω**. In SLIPs, a standard Boris algorithm is used:

*d*

**S**/

*dt*)

_{R}, which is equivalent to the electric field term in the Lorentz equation. Therefore, the straightforward algorithm is given by

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 (*d***S**/*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 (*d***S**/*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 (*d***S**/*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).

### 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

*a*

_{0}and the electron energy

*γ*

_{e}are such that the locally constant cross-field approximation (LCFA) is valid, i.e.,

*a*

_{0}≫ 1,

*χ*

_{e}≳ 1, the polarization- and spin-resolved emission rate for the NCS is given by

^{12,15,62}

*ξ*

_{1},

*ξ*

_{2},

*ξ*

_{3}), defined with respect to the axes $P\u03021=a\u0302\u2212n\u0302(n\u0302\u22c5a\u0302)$ and $P\u03022=n\u0302\xd7P\u03021$,

^{63}with the photon emission direction $n\u0302=pe/|pe|$ along the momentum

**p**

_{e}of the ultrarelativistic electron. The variables introduced in Eq. (49) are as follows:

*ω*

_{γ}is the emitted photon energy,

*ɛ*

_{i}is the electron energy before radiation, $a\u0302=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:

^{12,15,64}

**S**

_{f}, the initial spin-resolved radiation probability is obtained:

^{5,45,65}

^{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*

_{1}

*γ*

_{e}is emitted if the spin-dependent radiation probability in Eq. (55), $P\u2261d2W\u0304fi(\chi e,r1,\gamma e,Si)/dudt\u22c5\Delta t$, meets or exceeds

*r*

_{2}, following the so-called von Neumann rejection method. The final momenta of the electron and photon are given by

**p**

_{f}= (1 −

*r*

_{1})

**p**

_{i}and

**k**=

*r*

_{1}

**p**

_{i}, respectively. Next, the electron spin flips either parallel (spin-up) or antiparallel (spin-down) to the SQA with probabilities $Pflip\u2261Wfi\u2191/P$ and $Wfi\u2193/P$, respectively, where $Wfi\u2191,\u2193\u2261d2Wfi\u2191,\u2193/dudt\u22c5\Delta t$ from Eq. (57). In other words, the final spin

**S**

_{f}will flip parallel to the SQA if

*r*

_{3}<

*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

^{12,66}Then, the SQA is given by $S\u0304f/|S\u0304f|$, and the probability for the aligned case is given by $|S\u0304f|$ and that for the antiparallel case by $1\u2212|S\u0304f|$.

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 $\xi \u0302\xb1\u2261\xb1$($\xi \u03041$, $\xi \u03042$, $\xi \u03043$)/$\xi \u03040$, where $\xi \u03040\u2261(\xi \u03041)2+(\xi \u03042)2+(\xi \u03043)2$. The probabilities of photon emission in these states, $Wfi\xb1$, are given by Eq. (49). A stochastic procedure is defined using the fourth random number *r*_{4}: if $Wfi+/W\u0304fi\u2265r4$, the polarization state $\xi \u0302+$ will be chosen; otherwise, the polarization state will be assigned as $\xi \u0302\u2212$. Here, $W\u0304fi\u2261WRF0$ and $Wfi\xb1\u2261WRF0+\u2211j=1,3\xi j\xb1Fj$.

Between photon emissions, the electron dynamics in the external laser field are described by the Lorentz equation *d***p**/*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**

*. Besides, in this simulation, interference effects between emissions in adjacent coherent lengths (*

_{i}*l*

_{f}≃

*λ*

_{L}/

*a*

_{0}) are negligible when the employed laser intensity is ultrastrong, i.e.,

*a*

_{0}≫ 1. Therefore, the photon emissions occurring in each coherent length are independent of each other.

Examples of the electron dynamics and spin can be seen in Fig. 8: clearly, the average value matches the MLL equations for dynamics and the MLL + radiative T-BMT equations for spins. The beam evolution is also shown in Fig. 9. The energy spectra of electrons and photons, as well as the photon polarization, can be seen in Fig. 10.

#### 2. Definition and transformation of Stokes parameters

^{67}

*θ*

_{α}denotes the angle between

**P**and $P\u03021$, while

*θ*

_{β}represents the absolute phase. In quantum mechanics, the photon polarization state corresponding to

**P**can be described by the density matrix

**is the Pauli matrix, and**

*σ***= (**

*ξ**ξ*

_{1},

*ξ*

_{2},

*ξ*

_{3}) denotes the Stokes parameters, with

*ξ*

_{1}= sin(2

*θ*

_{α})cos(

*θ*

_{β}),

*ξ*

_{2}= sin(2

*θ*

_{α})sin(

*θ*

_{β}), and

*ξ*

_{3}= cos(2

*θ*

_{α}).

**E**and

**B**are the electric and magnetic fields. The two groups of polarization vectors are connected via rotation through an angle

*ψ*:

^{68,69}

### D. Nonlinear Breit–Wheeler (NBW) pair production

When the energy of a photon exceeds the rest mass of an electron–positron pair, i.e., *ω*_{γ} ≥ 2*m*_{e}*c*^{2}, and the photon is subjected to an ultraintense field *a*_{0} ≫ 1, the related nonlinear quantum parameter *χ*_{γ} can reach unity. Here, $\chi \gamma \u2261(e\u210f2/m3c4)|F\mu \nu k\nu |2$ and is approximately equal to $2a0\omega \gamma \xi 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

*G*

_{0}and polarization-related terms

*G*

_{1,2,3}are given by

*α*is the fine structure constant,

*ɛ*

_{γ},

*ɛ*

_{−}, and

*ɛ*

_{+}are the energies of the parent photon and the created electron and positron, respectively, $v\u0302+=v+/|v+|$ (with

**v**

_{+}the positron velocity), $a\u0302+=a+/|a+|$ (with

**a**

_{+}the positron acceleration in the rest frame of the positron), $b\u0302+=v+\xd7a+/|v+\xd7a+|$,

*ξ*

_{1},

*ξ*

_{2}, and

*ξ*

_{3}are the Stokes parameters of the

*γ*-photon, and

**S**

_{+}and

**S**

_{−}are the positron and electron spin vectors, respectively. K

_{n}is again the

*n*th-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 $(P\u03021,P\u03022,n\u0302)$ to the pair production frame $(P\u03021\u2032,P\u03022\u2032,n\u0302)$; see the transformations of the Stokes parameters in Sec. III C 2.

**S**

_{+}and the photon polarization

**is obtained as**

*ξ***S**

_{−}and the photon polarization is obtained as

#### 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*_{1}*ɛ*_{γ} when the probability density *P* ≡ *d*^{2}*W*_{pair}/*dɛ*_{+}*dt* · Δ*t* of pair production is greater than or equal to a random number *r*_{2} within the range [0, 1). Here, *d*^{2}*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*_{3} is used to sample the final spin states for the electron and positron. Note that here all random numbers are sampled uniformly from [0, 1), as in the NCS algorithm. An example of the production of secondary electrons and positrons resulting from a collision between a laser and an electron beam is illustrated in Fig. 12.

### E. High-energy bremsstrahlung

^{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:

*y*=

*ℏω*/

*E*

_{e}is the ratio of the energy of the emitted photon to that of the incident electron,

*r*

_{0}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.

^{76}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}

*β*=

*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.

*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*^{Δs/λ}, 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*_{1} is then generated and compared with *P*_{br}. If *r*_{1} < *P*_{br}, then a bremsstrahlung event is triggered. The energy of the resulting photon is determined by generating another random number *r*_{2}, which is then multiplied by *σ*_{br}(*E*_{e}) to obtain the energy ratio *y* through *σ*(*y*, *E*_{e}) = *σ*(*E*_{e})*r*_{2}. Finally, a photon with energy *ℏω* = *E*_{e}*y* 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^{5} primaries, colliding with a 5 mm Au target with *Z* = 79 and *ρ* = 19.3 g/cm^{3} and a 5 mm Al target with *Z* = 13 and *ρ* = 2.7 g/cm^{3}. We disabled the field updater and weighting procedure in the PIC code, and enabled only the particle pusher and bremsstrahlung MC module. The electron and photon spectra were found to be in good agreement with the Geant4 results, except for a slightly higher photon emission in the high-energy tail (which is due to the difference in the cross-section data). Figure 13 displays the spectra of electrons and photons from a 100 MeV electron bunch normally incident onto the Al and Au slabs, and similar distributions for a 1 GeV electron bunch are shown in Fig. 14.

### F. Vacuum birefringence

*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

*α*is the fine structure constant,

*m*is the electron mass,

*χ*

_{γ}is the nonlinear quantum parameter as defined earlier,

*x*= 4/[(1 −

*υ*

^{2})

*χ*

_{γ}], and Gi′(

*x*) is the derivative of the Scorer function. $Ered,\u22a5=E\u22a5+k\u0302\xd7B\u22a5$ is the transverse reduced field (acceleration field for electrons). The first and second rows in the ${\cdots \cdots}$ correspond to the eigenmodes parallel and perpendicular to the reduced field, respectively. After extraction of a factor

*χ*

_{γ}≪ 1, the imaginary part associated with pair production is negligible. We now define

*M*(

*χ*

_{γ}) and comparisons with the low-energy-limit (

*ω*

_{γ}≪

*m*) constants are given in Fig. 15.

*χ*

_{γ}≪ 1, the real part simplifies to

*χ*

_{γ}≪ 1. Note that these results are identical to those in Refs.

^{79–81}. For large

*χ*

_{γ}, two interpolated refractive indices are used.

*l*is the propagation length, and the VB effect is equivalent to a rotation of the Stokes parameters:

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

## IV. FRAMEWORK OF SLIPS

These physical processes have been incorporated into a spin-resolved laser–plasma interaction simulation code, known as SLIPs. The data structure and framework layout are illustrated in Figs. 17 and 18.

As depicted in Fig. 17, SLIPs utilizes a toml file to store simulation information, which is then parsed into a SimInfo structure that includes domainInfo, speciesInfo, boundaryInfo, laserInfo, pusherInfo, and other metadata. Subsequently, this metadata are employed to generate a SimBox that comprises all ParticleList and Fields, and to define the FieldSolver and EOMSolver and initialize QED processes.

The internal data structure of SLIPs is constructed using the open-source numerical library, Armadillo C++.^{83,84} String expressions are parsed using the ExprTk library.^{85} The data are then dumped using serial-hdf5 and merged with external Python scripts to remove ghost cells.

The spin-resolved processes, i.e., those tagged as Spin-QED in the diagram in Fig. 18, are implemented in conjunction with the Lorentz equation. In the coding, the Spin-QED part is arranged as a sequential series of processes. For example, Lorentz and T-BMT are followed by radiative correction, VB, NBW, and NCS with bremsstrahlung: $LorentzandT\u2212BMT\u21d2Radiativecorrection\u21d2VB\u21d2NBW\u21d2NCS$ and Bremss.

## V. POLARIZED PARTICLE SIMULATIONS

### A. Polarized electron/positron simulation

To simulate the generation of spin-polarized electrons, we utilized an elliptically polarized laser with an intensity *a*_{0} = 30, a wavelength *λ*_{0} = 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.

### B. Polarized *γ*-photons via NCS

The polarization state of emitted photons can be determined in spin/polarization-resolved NCS. Here, following Ref. 25, we utilized a linearly polarized (LP) laser to collide with an unpolarized electron bunch to generate LP *γ*-photons. Additionally, we used an LP laser to collide with a longitudinally polarized electron bunch to generate circularly polarized (CP) *γ*-photons, which were also observed in a previous study.^{12} The final polarization states of LP and CP *γ*-photons are presented in Figs. 20 and 21, respectively.

### C. Laser–plasma interactions

Finally, we present a simulation result demonstrating the interaction between an ultraintense laser with a normalized intensity *a*_{0} = 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 $e1\u2261\beta \xd7\beta \u0307$ and $e2\u2261n\u0302\xd7e1$, where $n\u0302$ 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% $(\xi 3\u0304\u22430.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.

## 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

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Feng Wan**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Chong Lv**: Formal analysis (equal); Resources (equal); Validation (equal); Writing – review & editing (equal). **Kun Xue**: Software (equal); Writing – review & editing (equal). **Zhen-Ke Dou**: Software (equal); Writing – review & editing (equal). **Qian Zhao**: Writing – review & editing (equal). **Mamutjan Ababekri**: Conceptualization (equal); Writing – review & editing (equal). **Wen-Qing Wei**: Writing – review & editing (equal). **Zhong-Peng Li**: Writing – review & editing (equal). **Yong-Tao Zhao**: Conceptualization (equal); Writing – review & editing (equal). **Jian-Xing Li**: Conceptualization (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data supporting this study's findings are available from the corresponding author upon reasonable request.

## REFERENCES

*Electromagnetic Processes at High Energies in Oriented Single Crystals*

*Radiation from Relativistic Electrons*

*γ*-rays via single-shot laser–electron interaction

*γ*-photon polarization in nonlinear Breit–Wheeler pair production and

*γ*polarimetry

*γ*-rays via laser–plasma interaction

*Plasma Physics via Computer Simulation*

*Computational Science—ICCS 2002*

*Progress in Ultrafast Intense Laser Science XVI*

*The Classical Theory of Fields*

*High-Power, High-Energy, and High-Intensity Laser Technology; and Research Using Extreme Light: Entering New Frontiers with Petawatt-Class Lasers*

*Relativistic Nonlinear Electrodynamics: The QED Vacuum and Matter in Super-Strong Radiation Fields*

*γ*photon splitting in wakefields

*Quantum Electrodynamics: Volume 4*

*Z*= 1–100

*γ*rays from laser-irradiated solids

Workshop Proceedings