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

Feng Wan 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 Chong Lv Department of Nuclear Physics, China Institute of Atomic Energy, P. O. Box 275(7), Beijing 102413, China Kun Xue 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 Zhen-Ke Dou 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 Qian Zhao 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 Mamutjan Ababekri 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 Wen-Qing Wei 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 Zhong-Peng Li 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 Yong-Tao Zhao 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 Jian-Xing Li [email protected] 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

(February 9, 2025)

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

^†^†preprint: AIP/123-QED

I Introduction

Laser-matter interactions can trigger strong-field quantum-electrodynamics (SF-QED) processes when the laser intensity $I_{0}$ reaches or exceeds $10^{22}\leavevmode\nobreak\ \mathrm{W/cm^{2}}$ Piazza et al. (2012); Bell and Kirk (2008). For example, when the laser intensity is in the order of $10^{21}$ - $10^{22}\leavevmode\nobreak\ \mathrm{W/cm^{2}}$ , i.e., the normalized peak laser field strength parameter $a_{0}\equiv eE_{0}/m_{e}c\omega_{0}\sim 10$ , electrons can be accelerated to GeV energies Gonsalves et al. (2019); Esarey, Schroeder, and Leemans (2009) (with Lorentz factor $\gamma_{e}\sim 10^{3}$ or higher) in a centimeter-long gas plasma, where $-e,m_{e}$ are the charge and mass of the electron, respectively, $E_{0},\omega_{0}$ are the field strength and angular frequency of the laser, respectively, and $c$ is the light speed in vacuum (here, for convenience, $\omega_{0}=\frac{2\pi c}{\lambda_{0}}$ and the wavelength of the laser $\lambda_{0}=1\mu\mathrm{m}$ are assumed). 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^{\prime}\simeq 2\gamma a_{0}\sim 10^{4}$ - $10^{5}$ . Such a field strength is close to the QED critical field strength (Schwinger critical field strength) $E_{\mathrm{Sch}}\equiv\frac{m_{e}^{2}c^{3}}{e\hbar}$ , i.e., $a_{\mathrm{Sch}}=\frac{m_{e}c^{2}}{\hbar\omega_{0}}\simeq 4.1\times 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(\chi,f,g)$ , with $\chi\equiv\frac{e\sqrt{(F_{\mu\nu}p^{\mu})^{2}}}{m^{3}}\sim a^{\prime}/a_{\mathrm{Sch}}$ , $f\equiv\frac{e^{2}F_{\mu\nu}F^{\mu\nu}}{4m^{4}}\sim\frac{(\bm{a}_{E}^{2}-\bm{a}_{B}^{2})}{4a^{2}_{\mathrm{Sch}}}$ and $g\equiv\frac{e^{2}F_{\mu\nu}F^{\mu\nu*}}{4m^{4}}\sim\frac{\bm{a}_{E}\cdot\bm{a}_{B}}{4a^{2}_{\mathrm{Sch}}}$ (here, $\bm{a}_{E,B}$ denote the normalized field strength of electric and magnetic component, respectively) Ritus (1985); Baier, Katkov, and Strakhovenko (1998). For most cases of weak field ( $a_{0}\ll a_{\mathrm{Sch}}$ ) condition, $f,g\ll\chi^{2}$ , and $W(\chi,f,g)\sim W(\chi)$ , i.e., the probability only depends on a single parameter $\chi$ . For electrons/positrons, the nonlinear Compton scattering (NCS, $e+n\omega_{L}\rightarrow e^{\prime}+\omega_{\gamma}$ ) is the dominant nonlinear QED process in the strong field regime, and for photons, the nonlinear Breit-Wheeler pair production (NBW, $\omega_{\gamma}+n\omega_{L}\rightarrow e^{+}+e^{-}$ ) is the dominant one, where $\omega_{L},\omega_{\gamma}$ denote the laser photon and the emitted $\gamma$ photon, respectively, and $n$ is the photon absorption number.

Apart from these kinetic effects, the 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. $\mathrm{d}{\bf s}/\mathrm{d}t\propto{\bf B}\times{\bf s}$ , where $\bf{s}$ denotes the classical spin vector Jackson (2021). In storage rings, due to the radiation reaction, the spin of an electron/positron will flip to the direction parallel/antiparallel to the external magnetic field, i.e., the Sokolov-Ternov effect (an unpolarized electron beam will be polarized to a degree of $\sim 92.5\%$ ) A. A. Sokolov , and a similar process also occurs in the NCS Del Sorbo et al. (2017); Li et al. (2019); King and Tang (2020). 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 the transient fields in plasmas Li et al. (2020); Gong, Hatsagortsyan, and Keitel (2021). 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 $\sim 30\%$ Wan et al. (2020a), and will be inherited by the subsequent QED cascade. For example, in the laser-plasma/beam interactions, the polarization degree for linearly polarized (LP) photons is about $60\%$ or higher and for circular polarized (CP) $\gamma$ photons can reach 59% when employing longitudinally polarized primariesXue et al. (2020); King and Tang (2020); Tang, King, and Hu (2020); Song, Wang, and Li (2022).

Analytical solutions in ultraintense laser-matter interactions are scarce due to the high nonlinearity and complexity of the problem. Moreover, the micro-level processes such as ionization, recombination and Coulomb collisions, etc., coupled with the complicated configurations of lasers and plasmas make the explicit derivation almost impossible. Fortunately, computer simulation methods provide alternative and more robust tools to study those unsolvable processes even in more realistic situations Arber et al. (2015). In general, simulation methods for laser-plasma (ionized matter) interactions can be categorized as kinetic or fluid simulations, and specifically kinetic methods include the Fokker-Planck equation (F-P) (or the Vlasov equation for the collisionless case) and the particle-in-cell (PIC) method, and the fluid method mainly uses the magnetohydrodynamic equations (MHD) Birdsall and Langdon (2018). Among these methods, both F-P and MHD discretize the momentum space of particles and are prone to the nonphysical multi-stream 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 Arber et al. (2015); Birdsall and Langdon (2018); Gonoskov et al. (2015).

Thanks to the emerging PIC simulation methods, the development of parallelism and large-scale cluster deployment, simulations of laser-plasma wakefield acceleration, laser-ion acceleration, THz radiation and also SF-QED, etc., have become accessible for general laser-plasma scientists Fonseca et al. (2002); Burau et al. (2010); Arber et al. (2015); Derouillat et al. (2018a); Wu et al. (2020). However, the spin and polarization properties of the plasma particles and QED products are not widely introduced to the mainstream due to the lack of appropriate algorithms. In some recent studies, the spin and polarization resolved SF-QED processes have been investigated in the laser-beam colliding configurations and have shown that these processes are prominent in generating polarized beams Wan et al. (2020a); King and Tang (2020); Tang, King, and Hu (2020); Song, Wang, and Li (2022); Li et al. (2019, 2022). And these local-constant-approximated version of these 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 the spin/polarization-dependent SF-QED processes are given in detail and have been coded into our PIC code SLIPs (which stands for “Spin-resolved Laser Interaction with Plasma simulation code”). These 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 available or will be achievable in the near future multi-petawatt (PW) to exawatt (EW) laser facilities Danson et al. (2015), such as Apollo Apo ; Zou et al. (2015), ELI ELI , SULF Gan et al. (2021) and SEL, etc. Throughout the paper, Gaussian units will be adopted, and all quantities are normalized as follows: time $t$ with $1/\omega$ (i.e., $t^{\prime}\equiv t/(1/\omega)=\omega t$ ), position $x$ with $1/k=\frac{\lambda}{2\pi}$ , momentum $p$ with $m_{e}c$ , velocity $v$ with $c$ , energy $\varepsilon$ with $m_{e}c^{2}$ , EM field $E,B$ with $\frac{m_{e}c\omega}{e}$ , force $F$ with $m_{e}c\omega$ , charge $q$ with $e$ , charge density $\rho$ with $k^{3}e$ , current density $J$ to $k^{3}ec$ , where $\lambda$ and $\omega=\frac{2\pi c}{\lambda}$ are the reference wavelength and frequency, respectively.

II PIC Algorithm

The simulation of laser-plasma interactions requires two essential components: the evolution of the EM field and the motion of particles. The corresponding governing equations are the Maxwell equations (either with $\bf{A}$ - $\phi$ or $\bf{E}$ - $\bf{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 Maxwell equations; see Figure. 1. Here, we review each part briefly (these algorithms are used in the SLIPs) and refer to the standard literature or textbooks for more details Arber et al. (2015); Birdsall and Langdon (2018).

Refer to caption — Figure 1: Standard particle-in-cell (PIC) loop with four kernel parts.

II.1 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 equation:

	$\displaystyle\frac{d{\mathbf{p}}}{dt}$	$\displaystyle=$	$\displaystyle\frac{q}{m}(\mathbf{E}+\bm{\beta}\times\mathbf{B}),$		(1)
	$\displaystyle\frac{d\mathbf{x}}{dt}$	$\displaystyle=$	$\displaystyle\frac{\mathbf{p}}{\gamma},$		(2)

where $\mathbf{p}\equiv\gamma m\mathbf{v}$ , $\mathbf{x}$ , $q$ , $m$ , $\gamma$ , $\mathbf{v}$ , and $\bm{\beta}\equiv\mathbf{v}/c$ are the momentum, position, charge, mass, Lorentz factor, velocity, and normalized velocity of the particle, respectively. These coupled equations are discretized using a leapfrog algorithm as

	$\displaystyle\frac{\mathbf{p}^{n+1/2}-\mathbf{p}^{n-1/2}}{\Delta t}$	$\displaystyle=$	$\displaystyle\frac{q}{m}\left(\mathbf{E}^{n}+\frac{\mathbf{p}^{n}}{\gamma^{n}}\times\mathbf{B}^{n}\right),$		(3)
	$\displaystyle\frac{\mathbf{x}^{n+1}-\mathbf{x}^{n}}{\Delta t}$	$\displaystyle=$	$\displaystyle\mathbf{v}^{n+1/2},$		(4)

and solved using the standard Boris rotation Buneman (1967); Boris and Shanny (1971); Qin et al. (2013):

$\displaystyle\mathbf{p}^{n-1/2}$	$\displaystyle=$	$\displaystyle\mathbf{p}^{-}-\frac{q\Delta t}{2m}\mathbf{E}^{n},$	(5)
$\displaystyle\mathbf{p}^{n+1/2}$	$\displaystyle=$	$\displaystyle\mathbf{p}^{+}+\frac{q\Delta t}{2m}\mathbf{E}^{n},$	(6)
$\displaystyle\mathbf{p}^{\prime}$	$\displaystyle=$	$\displaystyle\mathbf{p}^{-}+\mathbf{p}^{-}\times\bm{\tau},$	(7)
$\displaystyle\mathbf{p}^{+}$	$\displaystyle=$	$\displaystyle\mathbf{p}^{-}+\mathbf{p}^{\prime}\times\bm{\varsigma},$	(8)
$\displaystyle\bm{\tau}$	$\displaystyle=$	$\displaystyle\frac{q\Delta t}{2m\gamma^{n}}\mathbf{B}^{n},$	(9)
$\displaystyle\bm{\varsigma}$	$\displaystyle=$	$\displaystyle\frac{2\bm{\tau}}{1+\bm{\tau}^{2}},$	(10)

where $\gamma^{n}=\sqrt{1+(\mathbf{p}^{-})^{2}}=\sqrt{1+(\mathbf{p}^{+})^{2}}$ . The update 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 evolution.

II.2 Field solving

In the ultraintense laser-plasma interactions, the plasma particles are assumed to be distributed in the 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 in differential form are given by

$\displaystyle\nabla\cdot\mathbf{E}$	$\displaystyle=$	$\displaystyle\rho$	(11)
$\displaystyle\nabla\cdot\mathbf{B}$	$\displaystyle=$	$\displaystyle 0$	(12)
$\displaystyle\nabla\times\mathbf{E}$	$\displaystyle=$	$\displaystyle-\frac{\partial\mathbf{B}}{\partial t}$	(13)
$\displaystyle\nabla\times\mathbf{B}$	$\displaystyle=$	$\displaystyle\frac{\partial\mathbf{E}}{\partial t}+\mathbf{J}.$	(14)

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

By using Esierkepov’s method of current deposition Esirkepov (2001a), the current is calculated from the charge density via charge conservation, i.e., $\partial_{t}\rho+\nabla\cdot{\bf J}=0$ . Once the initial condition obeys Gauss’s law, $\nabla\cdot{\bf E}=\rho$ , Gauss’s law is automatically embedded. This can be verified with a gradient on Eq. (14) $0=\nabla\cdot(\nabla\times{\bf B})=\partial_{t}(\nabla\cdot{\bf E})+\nabla\cdot{\bf J}=\partial_{t}(\nabla\cdot{\bf E}-\rho)$ , 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. Here, we take the $E_{y}$ and $B_{z}$ components as examples:
1D case (squeezing the $y,z$ directions):

\displaystyle\begin{split}\frac{E_{y}^{n+1}-E_{y}^{n}}{\Delta t}\bigg{|}_{i+1/2}=&-\frac{B_{i+1}-B_{i}}{\Delta x}\bigg{|}_{z}^{n+1/2}-J_{y,i+1/2}^{n+1/2}\\ \frac{B_{z}^{n+1/2}-B_{z}^{n-1/2}}{\Delta t}\bigg{|}_{i}=&-\frac{E_{i+1/2}-E_{i-1/2}}{\Delta x}\bigg{|}^{n}_{y}\end{split}

(15)

2D case (squeezing the $z$ direction):

\displaystyle\begin{split}\frac{E_{y}^{n+1}-E_{y}^{n}}{\Delta t}\bigg{|}_{i+1/2,j}=&-\frac{B_{i+1,j}-B_{i,j}}{\Delta x}\bigg{|}_{z}^{n+1/2}-J^{n+1/2}_{y,i+1/2,j}\\ \frac{B_{z}^{n+1/2}-B_{z}^{n-1/2}}{\Delta t}\bigg{|}_{i+1/2,j}=&-\frac{E_{i+1/2,j}-E_{i+1/2,j}}{\Delta x}\bigg{|}_{y}^{n}+\\ &\frac{E_{i+1/2,j+1/2}-E_{i+1/2,j-1/2}}{\Delta y}\bigg{|}_{x}^{n+1/2}\end{split}

(16)

3D case:

\displaystyle\begin{split}\frac{E_{y}^{n+1}-E_{y}^{n}}{\Delta t}\bigg{|}_{i+1/2,j,k+1/2}=&-\frac{B_{i+1,j,k+1/2}-B_{i,j,k+1/2}}{\Delta x}\bigg{|}_{z}^{n+1/2}+\\ &\frac{B_{i+1/2,j,k+1}-B_{i+1/2,j,k}}{\Delta z}\bigg{|}_{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}}{\Delta t}\bigg{|}_{i,j,k+1/2}=&-\frac{E_{i+1/2,j,k}-E_{i-1/2,j,k}}{\Delta x}\bigg{|}_{y}^{n}+\\ &\frac{E_{i+1/2,j+1,k}-E_{i+1/2,j,k}}{\Delta y}\bigg{|}_{x}^{n}\end{split},

(17)

where 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.6.

II.3 Current deposition

We calculate the charge current density using Esirkepov’s method, which conserves charge by satisfying the Gauss law Esirkepov (2001b)

\partial_{t}\rho+\nabla\cdot\mathbf{J}=0,

(18)

and removes the need for Coulomb correction Birdsall and Langdon (2018). This algorithm computes the charge density at time step $t-\frac{1}{2}\Delta t$ and $t+\frac{1}{2}\Delta t$ on each grid cell from the particle positions and velocities, i.e.,

$\displaystyle\rho^{n+1/2}_{i,j,k}$	$\displaystyle=$	$\displaystyle\frac{1}{\Delta V}\sum_{r}W(\mathbf{x}^{n}_{r}+\frac{1}{2}\mathbf{v}^{n+1/2}\Delta t)q_{r},$	(19)
$\displaystyle\rho^{n-1/2}_{i,j,k}$	$\displaystyle=$	$\displaystyle\frac{1}{\Delta V}\sum_{r}W(\mathbf{x}^{n}_{r}-\frac{1}{2}\mathbf{v}^{n}\Delta t)q_{r},$	(20)
$\displaystyle\delta^{n}\rho$	$\displaystyle=$	$\displaystyle\rho^{n+1/2}-\rho^{n-1/2}$	(21)

where $r$ denotes the particle index, $|{\bf x}_{r}-{\bf x}_{i,j,k}|\leq(\Delta x,\Delta y,\Delta z)$ , and $\Delta V=\Delta x\Delta y\Delta z$ is the cell volume. We then interpolate the charge density to the current grid to obtain the current density; see Ref. Esirkepov (2001b) for more details.

II.4 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 $\mathbf{x}_{r}$ , we find its nearest grid point $(i,j,k)^{g}=\mathrm{floor}\left(\frac{\mathbf{x}_{r}}{\Delta\mathbf{x}}+\frac{1}{2}\right)$ and its nearest half grid point $(i,j,k)^{h}=\mathrm{floor}\left(\frac{\mathbf{x}_{r}}{\Delta\mathbf{x}}\right)$ , where $\Delta\mathbf{x}=(\Delta x,\Delta y,\Delta 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)\cdot F(i,j,k)$ , where $W(i,j,k)$ is the particle weighting function (see Sec. II.5 for more details) on the grid (half grid) $(i,j,k)$ and $F(i,j,k)$ is the field component of $\bf{E}$ or $\bf{B}$ on the spatial grid with proper staggering according to Figure. 2.

II.5 Particle shape function

The weighting function $W$ in the current and force deposition is determined by the form factor (shape factor) of the macro-particle, which is a key concept in modern PIC code algorithms. The form factor gives the macro-particle a finite size (composed of thousands of real particles) and reduces the nonphysical collisions Birdsall and Langdon (2018). 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 the unphysical noise and produce smoother results. We use the triangle shape function (triangular shape cloud, TSC) in each dimension Esirkepov (2001b)

\displaystyle W_{\rm spline}=\Bigg{\{}\begin{array}[]{ll}\frac{3}{4}-{\delta}^{2},&\mathrm{for}j\\ \frac{1}{2}\left(\frac{1}{2}\pm\delta\right),&\mathrm{for}j\pm 1\end{array},

(24)

where $\delta=\frac{x-X_{j}}{\Delta x}$ , $x$ is the particle position, $j$ the nearest grid/half grid number and $X_{j}\equiv j\Delta x$ . We obtain higher dimension functions by multiplying 1D shape function in each dimension: $W_{\mathrm{2D}}(i,j)=W_{x}(i)W_{y}(j)$ and $W_{\mathrm{3D}}(i,j,k)=W_{x}(i)W_{y}(j)W_{z}(k)$ .

II.6 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 Arber et al. (2015). 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 $\bf{E}$ and $\bf{B}$ are updated alternately at integer and half-integer time steps; and the position $\bf{x}$ and momentum $\bf{p}$ of particles are updated alternately as well; see Figure. 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/polarized version) that are relevant for laser-plasma interactions. The classical and quantum radiation correction to the Newton-Lorentz equations, namely the Laudau-Lifshitz (LL) equation and the modified Landau-Lifshitz (MLL) 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 (T-BMT) equation and the Radiative T-BMT equation; and their discretized algorithms are reviewed next. NCS with unpolarized and polarized version and their Monte-Carlo (MC) algorithms are reviewed. NBW with unpolarized and polarized version and their MC implementations are presented as well. Finally, the implementations of high-energy bremsstralung and vacuum birefringence in the conditions of weak pair productions ( $\chi_{\gamma}\lesssim 0.1$ ) are briefly discussed.

III.1 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 equation of motion (EOM) for charged particles is the Lorentz-Abraham-Dirac (LAD) equation Dirac (1938). However, this equation suffers from the runaway problem, as 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 Landau and Lifshitz (1999). The LL equation can be obtained from the LAD equation by applying iterative and order-reduction procedures Ekman, Heinzl, and Ilderton (2021); Ekman (2022), which are valid when the radiation force is much smaller than the Lorentz force. More importantly, in the limit of $\hbar\rightarrow 0$ , the QED results in the planewave background field are consistent with both the LAD equation and LL equation Ilderton and Torgrimsson (2013); Seipt and Thomas (2023). Depending on the value of the quantum nonlinear parameter $\chi_{e}$ (defined in Sec. III.1.1), the particle dynamics can be governed by either the LL equation or its quantum-corrected version Piazza et al. (2012); Landau and Lifshitz (1999); Neitz and Di Piazza (2014); Derouillat et al. (2018a).

III.1.1 Landau-Lifshitz equation

The Landau-Lifshitz equation can be employed when the radiation is relatively weak (weak radiation reaction, $\chi_{e}\ll 10^{-2}$ ) Landau and Lifshitz (1999):

\begin{split}\mathbf{F}_{\mathrm{RR,classical}}&=\frac{2e^{3}}{3mc^{3}}\bigg{\{}\gamma\bigg{[}\bigg{(}\frac{\partial}{\partial t}+\frac{\mathbf{p}}{\gamma m}\cdot\nabla\bigg{)}\mathbf{E}+\frac{\mathbf{p}}{\gamma mc}\times\bigg{(}\frac{\partial}{\partial t}+\frac{\mathbf{p}}{\gamma m}\cdot\nabla\bigg{)}\mathbf{B}\bigg{]}\\ &+\frac{e}{mc}\bigg{[}\mathbf{E}\times\mathbf{B}+\frac{1}{\gamma mc}\mathbf{B}\times(\mathbf{B}\times\mathbf{p})+\frac{1}{\gamma mc}\mathbf{E}(\mathbf{p\cdot E})\bigg{]}\\ &-\frac{e\gamma}{m^{2}c^{2}}\mathbf{p}\bigg{[}\bigg{(}\mathbf{E}+\frac{\mathbf{p}}{\gamma mc}\times\mathbf{B}\bigg{)}^{2}-\frac{1}{\gamma^{2}m^{2}c^{2}}(\mathbf{E\cdot p})^{2}\bigg{]}\bigg{\}},\end{split}

(25)

where all quantities are given in Gaussian units, and the dimensionless one is

\begin{split}\mathbf{F}_{\mathrm{RR,classical}}&=\frac{2}{3}\alpha_{f}\xi_{L}\bigg{\{}\gamma\bigg{[}\bigg{(}\frac{\partial}{\partial t}+\frac{\mathbf{p}}{\gamma}\cdot\nabla\bigg{)}\mathbf{E}+\frac{\mathbf{p}}{\gamma}\times\bigg{(}\frac{\partial}{\partial t}+\frac{\mathbf{p}}{\gamma}\cdot\nabla\bigg{)}\mathbf{B}\bigg{]}\\ &+\bigg{[}\mathbf{E}\times\mathbf{B}+\frac{1}{\gamma}\mathbf{B}\times(\mathbf{B}\times\mathbf{p})+\mathbf{E}(\mathbf{p\cdot E})\bigg{]}\\ &-\gamma\mathbf{p}\bigg{[}\bigg{(}\mathbf{E}+\frac{\mathbf{p}}{\gamma}\times\mathbf{B}\bigg{)}^{2}-\frac{1}{\gamma^{2}}(\mathbf{E\cdot p})^{2}\bigg{]}\bigg{\}},\end{split}

(26)

with $\alpha_{f}=\frac{e^{2}}{c\hbar}$ is the fine structure constant and $\xi_{L}=\frac{\hbar\omega}{m_{e}c^{2}}$ is the normalized reference photon energy. In the ultra-intense laser interacting with plasmas, the dominant contribution comes from the last two terms Tamburini et al. (2010). In the ultrarelativisitc limits, only the third term dominates the contribution, and the radiation reaction force can be simplified as

\mathbf{F}_{\mathrm{RR,cl}}\simeq\frac{2}{3}\alpha_{f}\frac{\chi^{2}_{e}}{\xi_{L}}\bm{\beta},

(27)

where $\chi_{e}=\frac{e\hbar}{m^{3}c^{4}}\sqrt{|F^{\mu\nu}p_{\nu}|^{2}}\equiv\xi_{L}\gamma_{e}\sqrt{(\mathbf{E}+\bm{\beta}\times\mathbf{B})^{2}-(\bm{\beta}\cdot(\bm{\beta}\cdot\mathbf{E}))^{2}}\simeq\gamma_{e}E_{\perp}\xi_{L}(1-\cos\theta)$ is a nonlinear quantum parameter signifies the strength of the NCS, with $\theta$ denote the angle between electron momentum and EM field wavevector, and $E_{\perp}$ denotes the perpendicular component of electric field. This reduced form gives the importance of the radiation reaction when one estimates the ratio between $F_{\mathrm{RR}}$ and Lorentz force $F_{L}$ :

\mathcal{R}\equiv|F_{\mathrm{RR}}|/|F_{\mathrm{L}}|\sim\frac{2}{3}\alpha_{f}\gamma_{e}\chi_{e}\simeq 2\times 10^{-8}a_{0}\gamma_{e}^{2}\leavevmode\nobreak\ (\mathrm{for\leavevmode\nobreak\ wavelength=1\leavevmode\nobreak\ \mu m}),

(28)

apparently, once $\gamma_{e}^{2}a_{0}\gtrsim 10^{6}$ , the radiation reaction force should be considered.

III.1.2 Modified Landau-Lifshitz 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_{e}c^{2}$ Bulanov et al. (2013a). Once $\chi_{e}$ is larger than $10^{-2}$ and above, the quantum nature of the radiation dominates the process. On one hand, the radiation spectra will be suppressed and deviate from the radiation force in LL equation; on the other hand, the radiation will be stochastic and discontinuous. However, when the stochasticity is not relevance for the detection and one only cares about the average effect (integrated spectra), a correction to the radiation force can be made, i.e., quantum correction Nikishov and Ritus (1964); Sokolov et al. (2009); Piazza, Hatsagortsyan, and Keitel (2010); Thomas et al. (2012)

\mathbf{F}_{\mathrm{RR,quantum}}=q(\chi)\mathbf{F}_{\mathrm{RR,classical}},

(29)

where

$\displaystyle q(\chi)$	$\displaystyle=$	$\displaystyle\frac{I_{\mathrm{QED}}}{I_{\mathrm{C}}},$	(30)
$\displaystyle I_{\mathrm{QED}}$	$\displaystyle=$	$\displaystyle mc^{2}\int c(k\cdot k^{\prime})\frac{dW_{fi}}{d\eta dr_{0}}dr_{0},$	(31)
$\displaystyle I_{\mathrm{C}}$	$\displaystyle=$	$\displaystyle\frac{2e^{4}E^{\prime 2}}{3m^{2}c^{3}},$	(32)

with $W_{fi}$ being the radiation probability Sokolov et al. (2010), $\eta=k_{0}z-\omega_{0}t$ , $r_{0}=\frac{2\left(k\cdot k^{\prime}\right)}{3\chi\left(k\cdot p_{i}\right)}$ , and $E^{\prime}$ the electric fields in the instantaneous frame of the electron. $p_{i}$ is the four-momentum of the electron before radiation. $k$ and $k^{\prime}$ are the four-wavevector of local EM field, and the radiated photon, respectively. See $q(\chi)$ in Figure. 4. Here, the ratio between the QED radiation power and the classical one, i.e., the re-scaling factor $q(\chi)$ , is the same with the factor in Ref. Bulanov et al. (2013b):

q(\chi_{e})\approx\frac{1}{\left[1+4.8(1+\chi_{e})\ln(1+1.7\chi_{e})+2.44\chi_{e}^{2}\right]^{2/3}},

(33)

q(\chi_{e})\approx\frac{1}{(1+8.93\chi_{e}+2.41\chi_{e}^{2})^{2/3}}.

(34)

In the ultrarelativistic limit, an alternative formula can be employed as Niel et al. (2018); Derouillat et al. (2018b)

\mathbf{F}_{\mathrm{RR,quantum}}=q(\chi)P_{\mathrm{cl}}\chi^{2}_{e}\bm{\beta}/\bm{\beta}^{2}c.

(35)

Apparently, once $\chi\gtrsim 10^{-2}$ , the quantum corrected version should be used.

III.1.3 Algorithms of the Radiative Pusher

Here, we plug the radiative correction (either classical or quantum corrected version) into the standard Boris pusher as follows Tamburini et al. (2010):

\frac{\mathbf{p}^{n+1/2}-\mathbf{p}^{n-1/2}}{\Delta t}=\mathbf{F}^{n}=\mathbf{F}^{n}_{L}+\mathbf{F}^{n}_{R}.

(36)

First we use the Boris step

\frac{\mathbf{p}_{L}^{n+1/2}-\mathbf{p}_{L}^{n-1/2}}{\Delta t}=\mathbf{F}_{L}^{n},

(37)

and then use the radiative correction push

\frac{\mathbf{p}^{n+1/2}_{R}-\mathbf{p}^{n-1/2}_{R}}{\Delta t}=\mathbf{F}_{R}^{n},

(38)

where $\mathbf{p}^{n-1/2}_{L}=\mathbf{p}^{n-1/2}_{R}=\mathbf{p}^{n-1/2}$ , and the final momentum is given by

\mathbf{p}^{n+1/2}=\mathbf{p}^{n+1/2}_{L}+\mathbf{p}^{n+1/2}_{R}-\mathbf{p}^{n-1/2}=\mathbf{p}^{n+1/2}_{L}+\mathbf{F}_{R}^{n}\Delta t.

(39)

With this algorithm, the Boris pusher is attained.

See a comparison between different solver calculated dynamics in Figure. 5. For the Lorentz equation without radiation, the particle momentum and energy are analytically given by Avetissian (2015)

	$\displaystyle{\bf p}(\tau)$	$\displaystyle=$	$\displaystyle{\bf p}_{0}-{\bf A}(\tau)+\hat{k}\frac{{\bf A}^{2}(\tau)-2{\bf p}_{0}\cdot{\bf A}(\tau)}{2(\gamma_{0}-{\bf p}_{0}\cdot\hat{k})}$		(40)
	$\displaystyle\gamma(\tau)$	$\displaystyle=$	$\displaystyle\gamma_{0}+\frac{{\bf A}^{2}(\tau)-2{\bf p}_{0}\cdot{\bf A}(\tau)}{2(\gamma_{0}-{\bf p}_{0}\cdot\hat{k})}$		(41)

where ${\bf A}(\tau)=-\int{\tau_{0}}^{\tau}{\bf E(\tau)}d\tau$ is the external field vector potential, $\tau$ is the proper time, $\hat{k}$ is the normalized wavevector of the field, $\gamma,{\bf p}$ , and $\gamma_{0},{\bf p}_{0}$ the instantaneous and initial (with notation of $0$ ) Lorentz factor and momentum of the particle, respectively. For a planewave with a temporal profile, the momentum and energy gain vanish as ${\bf A}(\infty)={\bf A}(-\infty)=0$ . The planewave solution with radiation reaction can be found in Ref. Piazza (2008). However, no explicit solution exists when the quantum correction term is included, as shown in Fig. 5.

III.2 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 electron/positron and $\gamma$ photons. The significance of this aspect has been highlighted in recent literature, particularly in the context of relativistic charged particles in EM waves and laser-matter interactions Hu and Keitel (1999); Walser et al. (2002). This issue can be addressed either by employing the computational Dirac solver Mocken and Keitel (2008) or by utilizing the Foldy-Wouthuysen transformation and the quantum operator formalism, such as the reduction of the Heisenberg equation to a classical precession equation Bauke et al. (2014); Wen, Keitel, and Bauke (2017). However, these approaches are not directly applicable to many-particle systems. Here and throughout this paper, the spin is defined as a unit vector $\mathbf{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 Bauke et al. (2014); Wen, Keitel, and Bauke (2017); Jackson (2021). 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 Landau-Lifshitz equation.

III.2.1 T-BMT equation

The non-radiative spin dynamics of an electron is given by

\begin{split}\left(\frac{{\rm d}{\bf S}}{{\rm d}t}\right)_{T}=\leavevmode\nobreak\ &\mathbf{S}\times\mathbf{\Omega}\equiv{\bf S}\times\left[-\left(\frac{g}{2}-1\right)\frac{\gamma_{e}}{\gamma_{e}+1}\left({\bm{\beta}}\cdot{\bf B}\right)\cdot{\bm{\beta}}\right.\\ &\left.+\left(\frac{g}{2}-1+\frac{1}{\gamma_{e}}\right){\bf B}-\left(\frac{g}{2}-\frac{\gamma_{e}}{\gamma_{e}+1}\right)\bm{\beta}\times{\bf E}\right],\end{split}

(42)

where ${\bf E}$ and ${\bf B}$ are the normalized electric and magnetic fields, respectively, and $g$ is the electron Landé factor, respectively. Since this equation is a pure rotation around the precession frequency of $\mathbf{\Omega}$ , the Boris rotation is much more preferable than other solver for ordinary differential equations (for instance, the Runge-Kutta, etc.). Here, $\mathbf{\Omega}$ plays the role of $\mathbf{B}/\gamma$ in the Eqs. (3) and (5-10). For other particle species, both the charge, mass and Landé factor for that species should be employed.

III.2.2 Radiative T-BMT equation

When the radiation damping is no longer negligible, the radiation can also affect the spin dynamics. In the weakly radiation regime, this radiation induced modification of the spin dynamics can be handled with a similar way as in the Landau-Lifshitz equation. Here, the modified version of T-BMT equation or the radiative T-BMT equation is thus given by

\frac{{\rm d}{\bf S}}{{\rm d}t}=\left(\frac{{\rm d}{\bf S}}{{\rm d}t}\right)_{T}+\left(\frac{{\rm d}{\bf S}}{{\rm d}t}\right)_{R},

(43)

with the first (labeled with “T”) and second (labeled with “R”) terms denote the non-radiative precession in Eq. (42), and radiative correction, respectively. The radiative term is given by

\left(\frac{{\rm d}{\bf S}}{{\rm d}t}\right)_{R}=-P\left[\psi_{1}(\chi){\bf S}+\psi_{2}(\chi)({\bf S}\cdot{\bm{\beta}}){\bm{\beta}}+\psi_{3}(\chi){\hat{\bf n}}_{B}\right],

(44)

where, $P=\frac{\alpha_{f}}{\sqrt{3}\pi\gamma_{e}\xi_{L}}$ , $\psi_{1}(\chi_{e})$ = $\int_{0}^{\infty}u^{\prime\prime}{\rm d}u{\rm K}_{\frac{2}{3}}(u^{\prime})$ , $\psi_{2}(\chi_{e})$ = $\int_{0}^{\infty}u^{\prime\prime}{\rm d}u\int_{u^{\prime}}^{\infty}{\rm d}x{\rm K}_{\frac{1}{3}}(x)$ - $\psi_{1}(\chi_{e})$ , $\psi_{3}(\chi)$ = $\int_{0}^{\infty}u^{\prime\prime}{\rm d}u{\rm K}_{\frac{1}{3}}(u^{\prime})$ , $u^{\prime}=2u/3\chi_{e}$ , $u^{\prime\prime}=u^{2}/(1+u)^{3}$ , and ${\rm K}_{n}$ is the $n$ th-order modified Bessel function of the second kind, $\hat{{\bf n}}_{B}={\bm{\beta}}\times\hat{{\bf a}}$ , ${\bm{\beta}}$ and $\hat{{\bf a}}$ denote the normalized velocity and acceleration vector Baĭer (1972); Guo et al. (2020).

III.2.3 Algorithms of simulating the spin precession

The simulation algorithms of spin precession are quite similar to the cases of EOM, i.e., Lorentz equation and radiative EOM, i.e., LL/MLL equations. Therefore, the T-BMT equation is simulated via the Boris rotation without the pre- and post- acceleration term, but with only the rotation term $\mathbf{\Omega}$ . In SLIPs, a standard Boris algorithm is used

$\displaystyle\mathbf{S}^{\prime}$	$\displaystyle=$	$\displaystyle\mathbf{S}^{n-1/2}+\mathbf{S}^{n-1/2}\times\mathbf{t},$	(45)
$\displaystyle\mathbf{S}^{n+1/2}$	$\displaystyle=$	$\displaystyle\mathbf{S}^{n-1/2}+\mathbf{S}^{\prime}\times\mathbf{o},$	(46)
$\displaystyle\mathbf{t}$	$\displaystyle=$	$\displaystyle\frac{q\Delta t}{2}\mathbf{\Omega}^{n},$	(47)
$\displaystyle\mathbf{o}$	$\displaystyle=$	$\displaystyle\frac{2\mathbf{t}}{1+t^{2}}.$	(48)

For the radiative T-BMT equation, there will be an extra term $\left(\frac{d\mathbf{S}}{dt}\right)_{R}$ which is equivalent to the electric field term in the Lorentz equation. Therefore, the straightforward algorithm is given by

\mathbf{S}^{n-1/2}_{T}=\mathbf{S}^{n-1/2}+\frac{\Delta t}{2}\left(\frac{d\mathbf{S}}{dt}\right)_{R},

(49)

Boris T-BMT Eqs. (45-48),

\mathbf{S}^{n+1/2}=\mathbf{S}^{n+1/2}_{T}+\frac{\Delta t}{2}\left(\frac{d\mathbf{S}}{dt}\right)_{R}.

(50)

Figure. 6 shows the comparison between the T-BMT and radiative T-BMT equations for different cases: Lorentz equation + BMT equation (“A”), Lorentz equation + M-BMT equation (“B”), LL equation + M-BMT equation (“C”), and MLL equation + M-BMT equation (“D”). The evolution of each spin component depends on different terms. In our setup, the magnetic field is along the $z$ direction, 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 due to the precession and are conserved in Fig. 6(d). In the case with only the spin radiation reaction, case “B”, $S_{x}$ is strongly damped by the term $\left(\frac{d{\bf S}}{dt}\right)_{R}$ . $S_{y}$ and $S_{z}$ oscillate due to the combined effects of precession and radiation reaction, as shown in Figs. 6(a) and (b). When both spin and momentum radiation reactions are included, case “C” (LL equation), the particle’s momentum and energy decrease, i.e., $\gamma_{e}$ decreases, which lowers the spin radiation reaction term $\left(\frac{d{\bf S}}{dt}\right)_{R}(\chi_{e})$ and the damping of $S_{x}$ and $S_{z}$ (see Fig. 6(c) for the comparison of “B”, “D” and “C” in terms of $S_{z}$ amplitude). Simultaneously, the precession term $\left(\frac{d{\bf S}}{dt}\right)_{T}\propto B/\gamma_{e}$ grows with decreasing $\gamma_{e}$ , which amplifies the oscillation of $S_{y}$ , as shown by the contrast of “B” (Lorentz), “D” (MLL) and “C” (LL) in Fig. 6(b).

III.3 Nonlinear Compton scattering

When the radiation is strong ( $\chi_{e}\gtrsim 0.1$ ), the stochastic nature of the radiation can no longer be neglected in the laser-beam/plasma interactions. And 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., the nonlinear Compton scattering (NCS) Bell and Kirk (2008); Kirk, Bell, and Arka (2009); Ridgers et al. (2014). Therefore, the radiation reaction and photon emission process will be calculated via the MC simulation based on the NCS probabilities. Besides, the spin of electron/positron and polarization of the NCS photons will be also included in the MC simulations.

III.3.1 Spin-resolved/summed nonlinear Compton scattering

When the laser intensity $a_{0}$ and the electron energy $\gamma_{e}$ permits the local-constant-cross field approximate (LCFA), i.e., $a_{0}\gg 1$ , $\chi_{e}\gtrsim 1$ , the polarization- and spin-resolved emission rate for the NCS is given by Li et al. (2020); Xue et al. (2020); Liu et al. (2022)

\frac{{\rm d^{2}}W_{fi}}{{\rm d}u{\rm d}t}=\frac{W_{R}}{2}\left(F_{0}+\xi_{1}F_{1}+\xi_{2}F_{2}+\xi_{3}F_{3}\right),

(51)

where, the photon polarization is represented by the Stokes parameters ( $\xi_{1}$ , $\xi_{2}$ , $\xi_{3}$ ), defined with respect to the axes $\hat{{\bf P}}_{1}=\hat{\bf a}-\hat{\bf n}(\hat{\bf n}\cdot\hat{\bf a})$ and $\hat{{\bf P}}_{2}=\hat{\bf n}\times\hat{\bf P}_{1}$ Green and Harvey (2015), with the photon emission direction $\hat{\bf n}={\bf p}_{e}/|{\bf p}_{e}|$ along the momentum ${\bf p}_{e}$ of the ultrarelativistic electron. The variables introduced in Eq. (51) read:

$\displaystyle F_{0}$	$\displaystyle=$	$\displaystyle-(2+u)^{2}\left[{\rm IntK}_{\frac{1}{3}}(u^{\prime})-2{\rm K}_{\frac{2}{3}}(u^{\prime})\right](1+{S}_{if})+u^{2}(1-{S}_{if})$	(52)
		$\displaystyle\left[{\rm IntK}_{\frac{1}{3}}(u^{\prime})+2{\rm K}_{\frac{2}{3}}(u^{\prime})\right]+2u^{2}{S}_{if}{\rm IntK}_{\frac{1}{3}}(u^{\prime})-(4u+2u^{2})$
		$\displaystyle({\bf S}_{f}+{\bf S}_{i})\cdot\left[\hat{\bf n}\times\hat{{\bf a}}\right]{\rm K}_{\frac{1}{3}}(u^{\prime})-2u^{2}({\bf S}_{f}-{\bf S}_{i})\cdot\left[\hat{\bf n}\times\hat{{\bf a}}\right]{\rm K}_{\frac{1}{3}}(u^{\prime})$
		$\displaystyle-4u^{2}\left[{\rm IntK}_{\frac{1}{3}}(u^{\prime})-{\rm K}_{\frac{2}{3}}(u^{\prime})\right]({\bf S}_{i}\cdot\hat{\bf n})({\bf S}_{f}\cdot\hat{\bf n}),$

$\displaystyle F_{1}$	$\displaystyle=$	$\displaystyle-2u^{2}{\rm IntK}_{\frac{1}{3}}(u^{\prime})\left\{({\bf S}_{i}\cdot\hat{\bf a}){\bf S}_{f}\cdot\left[\hat{\bf n}\times\hat{\bf a}\right]+({\bf S}_{f}\cdot\hat{\bf a}){\bf S}_{i}\cdot\left[\hat{\bf n}\times\hat{\bf a}\right]\right\}+$	(53)
		$\displaystyle 4u\left[({\bf S}_{i}\cdot\hat{\bf a})(1+u)+({\bf S}_{f}\cdot\hat{\bf a})\right]{\rm K}_{\frac{1}{3}}(u^{\prime})+$
		$\displaystyle 2u(2+u)\hat{\bf n}\cdot[{\bf S}_{f}\times{\bf S}_{i}]{\rm K}_{\frac{2}{3}}(u^{\prime}),$

$\displaystyle F_{2}$	$\displaystyle=$	$\displaystyle-\left\{2u^{2}\left\{({\bf S}_{i}\cdot\hat{\bf n}){\bf S}_{f}\cdot\left[\hat{\bf n}\times\hat{\bf a}\right]+({\bf S}_{f}\cdot\hat{\bf n}){\bf S}_{i}\cdot\left[\hat{\bf n}\times\hat{\bf a}\right]\right\}+2u(2+u)\right.$	(54)
		$\displaystyle\left.\hat{\bf a}\cdot[{\bf S}_{f}\times{\bf S}_{i}]\right\}{\rm K}_{\frac{1}{3}}(u^{\prime})-4u\left[({\bf S}_{i}\cdot\hat{\bf n})+({\bf S}_{f}\cdot\hat{\bf n})(1+u)\right]$
		$\displaystyle{\rm IntK}_{\frac{1}{3}}(u^{\prime})+4u(2+u)\left[({\bf S}_{i}\cdot\hat{\bf n})+({\bf S}_{f}\cdot\hat{\bf n})\right]{\rm K}_{\frac{2}{3}}(u^{\prime}),$

$\displaystyle F_{3}$	$\displaystyle=$	$\displaystyle 4\left[1+u+(1+u+\frac{u^{2}}{2}){S}_{if}-\frac{u^{2}}{2}({\bf S}_{i}\cdot\hat{\bf n})({\bf S}_{f}\cdot\hat{\bf n})\right]{\rm K}_{\frac{2}{3}}(u^{\prime})$	(55)
		$\displaystyle+2u^{2}\left\{{\bf S}_{i}\cdot\left[\hat{\bf n}\times\hat{\bf a}\right]{\bf S}_{f}\cdot\left[\hat{\bf n}\times\hat{\bf a}\right]-({\bf S}_{i}\cdot\hat{\bf a})({\bf S}_{f}\cdot\hat{\bf a})\right\}{\rm IntK}_{\frac{1}{3}}(u^{\prime})$
		$\displaystyle-4u\left[(1+u){\bf S}_{i}\left[\hat{\bf n}\times\hat{\bf a}\right]+{\bf S}_{f}\left[\hat{\bf n}\times\hat{\bf a}\right]\right]{\rm K}_{\frac{1}{3}}(u^{\prime}),$

where $W_{R}=\alpha_{f}/\left[8\sqrt{3}\pi\xi_{L}{\left(1+u\right)^{3}}\right]$ , $u^{\prime}=2u/3\chi$ , $u=\omega_{\gamma}/\left(\varepsilon_{i}-\omega_{\gamma}\right)$ , ${\rm IntK}_{\frac{1}{3}}(u^{\prime})\equiv\int_{u^{\prime}}^{\infty}{\rm d}z{\rm K}_{\frac{1}{3}}(z)$ , $\omega_{\gamma}$ the emitted photon energy, $\varepsilon_{i}$ the electron energy before radiation, $\hat{{\bf a}}={\bf a}/|{\bf a}|$ the direction of the electron acceleration ${\bf a}$ , ${\bf S}_{i}$ and ${\bf S}_{f}$ denote the electron spin vectors before and after radiation, respectively, $|{\bf S}_{i,f}|=1$ , and ${S}_{if}\equiv{\bf S}_{i}\cdot{\bf S}_{f}$ .

Summing over the photon polarization, the electron spin-resolved emission probability can be written as Li et al. (2020); Wan et al. (2020b); Xue et al. (2020):

	$\displaystyle\frac{{\rm d}^{2}W_{fi}}{{\rm d}u{\rm d}t}=W_{R}\left\{-(2+u)^{2}\left[{\rm IntK}_{\frac{1}{3}}(u^{\prime})-2{\rm K}_{\frac{2}{3}}(u^{\prime})\right](1+{S}_{if})+\right.$
	$\displaystyle\left.u^{2}\left[{\rm IntK}_{\frac{1}{3}}(u^{\prime})+2{\rm K}_{\frac{2}{3}}(u^{\prime})\right](1-{\bf S}_{if})+2u^{2}{S}_{if}{\rm IntK}_{\frac{1}{3}}(u^{\prime})-\right.$
	$\displaystyle\left.(4u+2u^{2})({\bf S}_{f}+{\bf S}_{i})\left[{\bf n}\times\hat{{\bf a}}\right]{\rm K}_{\frac{1}{3}}(u^{\prime})-2u^{2}({\bf S}_{f}-{\bf S}_{i})\left[{\bf n}\times\hat{{\bf a}}\right]\right.$
	$\displaystyle\left.{\rm K}_{\frac{1}{3}}(u^{\prime})-4u^{2}\left[{\rm IntK}_{\frac{1}{3}}(u^{\prime})-{\rm K}_{\frac{2}{3}}(u^{\prime})\right]({\bf S}_{i}\cdot{\bf n})({\bf S}_{f}\cdot{\bf n})\right\}.$		(56)

Summing over the final states ${\bf S}_{f}$ , the initial spin-resolved radiation probability is obtained:

	$\displaystyle\frac{{\rm d}^{2}{\overline{W}}_{fi}}{{\rm d}u{\rm d}t}$	$\displaystyle=$	$\displaystyle 8W_{R}\left\{-(1+u){\rm IntK}_{\frac{1}{3}}(u^{\prime})+(2+2u+u^{2}){\rm K}_{\frac{2}{3}}(u^{\prime})\right.$		(57)
			$\displaystyle\left.-u{\bf S}_{i}\cdot\left[{\bf n}\times\hat{{\bf a}}\right]{\rm K}_{\frac{1}{3}}(u^{\prime})\right\}.$		(57)

And by averaging the electron initial spin, one obtains the widely used radiation probability for the unpolarized initial particles Ritus (1979); Nikishov and Ritus (1964); Ritus (1985).

During the photon emission simulation, 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 choosen to obtains the maximum transition probabily, which is along the energy-resolved average polarization

\displaystyle\mathbf{S}_{\rm f}^{\rm R}=\frac{\mathbf{g}}{(w+\mathbf{f}\cdot\mathbf{S}_{i})}.

(58)

These are obtained by summing over the photon polarization and keeps the dependence on initial and final spin of electrons:

\displaystyle\frac{{\rm d}^{2}W_{\rm rad}}{{\rm d}u{\rm d}t}

\displaystyle=

\displaystyle{W_{\rm r}}(w+\mathbf{f}\cdot\mathbf{S}_{i}+\mathbf{g}\cdot\mathbf{S}_{f}),

(59)

where

$\displaystyle w$	$\displaystyle=$	$\displaystyle-(1+u){\rm K}_{\frac{1}{3}}(\rho^{\prime})+(2+2u+u^{2}){\rm K}_{\frac{2}{3}}(\rho^{\prime}),$
$\displaystyle\mathbf{f}$	$\displaystyle=$	$\displaystyle u{\rm IntK}_{\frac{1}{3}}(\rho^{\prime})\hat{{\mathbf{v}}}\times\hat{{\mathbf{a}}},$
$\displaystyle\mathbf{g}$	$\displaystyle=$	$\displaystyle-(1+u)\left[{\rm K}_{\frac{1}{3}}(\rho^{\prime})-2{\rm K}_{\frac{2}{3}}(\rho^{\prime})\right]\mathbf{S}_{i}-(1+u)u$
		$\displaystyle\times{\rm IntK}_{\frac{1}{3}}(\rho^{\prime})\hat{{\mathbf{v}}}\times\hat{{\mathbf{a}}}-u^{2}\left[{\rm K}_{\frac{1}{3}}(\rho^{\prime})-{\rm K}_{\frac{2}{3}}(\rho^{\prime})\right](\mathbf{S}_{i}\cdot\hat{{\mathbf{v}}})\hat{{\mathbf{v}}}.$

Conversely, without emission, the SQA aligns with another SQA Li et al. (2020); Yokoya et al. (2011). 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}\in[0,1)$ . The procedure is as follows: First, at each simulation time step $\Delta t$ , a photon with energy $\omega_{\gamma}=r_{1}\gamma_{e}$ is emitted if the spin-dependent radiation probability in Eq.(57), $P\equiv\mathrm{d}^{2}\overline{W}_{fi}(\chi_{e},r_{1},\gamma_{e},{\bf S}_{i})/\mathrm{d}u\mathrm{d}t\cdot\Delta t$ , meets or exceeds $r_{2}$ , following the “von Neumann’s rejection method”. The final momentum for the electron and photon is given by ${\bf p}_{f}=(1-r_{1}){\bf p}_{i}$ and $\hbar{\bf k}=r_{1}{\bf p}_{i}$ , respectively. Next, the electron spin flips either parallel (spin-up) or antiparallel (spin-down) to the SQA with probabilities of $P_{\mathrm{flip}}\equiv W_{fi}^{\uparrow}/P$ and $W_{fi}^{\downarrow}/P$ , respectively, where $W_{fi}^{\uparrow,\downarrow}\equiv\mathrm{d^{2}}W_{fi}^{\uparrow,\downarrow}/\mathrm{d}u\mathrm{d}t\cdot\Delta t$ from Eq. (59). In other words, the final spin ${\bf S}_{f}$ will flip parallel to the SQA if $r_{3}<P_{\mathrm{flip}}$ , and vice versa; see the flow chart of the NCS in Figure. 7. In the alternative scenario, i.e., no photon is emitted, the average final spin is given by $\overline{\bf S}_{f}=\frac{{\bf S}_{i}(1-W\Delta t)-{\bf f}\Delta t}{1-(W+{\bf f}\cdot{\bf S}_{i})\Delta t}$ , where, $W\equiv 16W_{R}(-(1+u)\mathrm{IntK}_{1/3}(u^{\prime})+(2+2u+u^{2})\mathrm{K}_{2/3}(u^{\prime}))$ , and ${\bf f}\equiv-16W_{R}({\bf n}\times{\bf\hat{a}}\mathrm{K}_{1/3}(u^{\prime}))$ Yokoya et al. (2011); Li et al. (2020). Then the SQA is given by $\overline{\bf S}_{f}/|\overline{\bf S}_{f}|$ , and the probability for the aligned case is given by $|\overline{\bf S}_{f}|$ , and $1-|\overline{\bf S}_{f}|$ for the anti-parallel case.

Finally, the polarization of the emitted photon is determined by considering that the average polarization is in a mixed state. The basis for the emitted photon is chosen as two orthogonal pure states with the Stokes parameters $\hat{\bm{\xi}}^{\pm}\equiv\pm$ ( $\overline{\xi}_{1}$ , $\overline{\xi}_{2}$ , $\overline{\xi}_{3}$ )/ $\overline{\xi}_{0}$ , where $\overline{\xi}_{0}\equiv\sqrt{(\overline{\xi}_{1})^{2}+(\overline{\xi}_{2})^{2}+(\overline{\xi}_{3})^{2}}$ . The probabilities for the photon emission in these states $W_{fi}^{\pm}$ are given by Eq. (51). A stochastic procedure is defined using the 4th random number $r_{4}$ : if $W_{fi}^{+}/\overline{W}_{fi}\geq r_{4}$ , the polarization state $\hat{\bm{\xi}}^{+}$ will be chosen; otherwise, the polarization state will be assigned as $\hat{\bm{\xi}}^{-}$ . Here, $\overline{W}_{fi}\equiv W_{R}F_{0}$ and $W^{\pm}_{fi}\equiv W_{R}(F_{0}+\sum_{j=1,3}\xi^{\pm}_{j}F_{j})$ .

Between photon emissions, the electron dynamics in the external laser field are described by the Lorentz equations, $d{\bf p}/dt=-e(\bf{E}+\bm{\beta}\times\bf{B})$ , and are simulated using the Boris rotation method, as shown in Eqs. (5-10). Due to the small emission angle for an ultrarelativistic electron, the photon is assumed to be emitted along the parental electron velocity, i.e., ${\bf p}_{f}\approx(1-\omega{\gamma}/|{\bf p}_{i}|){\bf p}i$ . Besides, in this simulation, interference effects between emissions in adjacent coherent lengths ( $l_{f}\simeq\lambda_{L}/a_{0}$ ) are negligible when the employed laser intensity is ultrastrong, i.e., $a_{0}\gg 1$ . Therefore, the photon emissions happening in each coherent length are independent of each other.

Examples of the electron dynamics and spin can be seen in Figure.8, apparently, the average value matches the MLL equations for dynamics and MLL + M-BMT equation for spins. Besides, the beam evolution is shown in Figure. 9. The energy spectra of electrons and photons, as well as the photon polarization, can be observed in Figure. 10.

III.3.2 Definition and Transformation of Stokes Parameters

In the context of NCS and the subsequent nonlinear Breit-Wheeler (NBW) pair production, the polarization state of a photon can be characterized by the polarization unit vector $\hat{\bf P}$ , which functions as the spin component of the photon wave function. An arbitrary polarization $\hat{\bf P}$ can be represented as a superposition of two orthogonal basis vectors Berestetskii, Lifshitz, and Pitaevskii (1982):

\displaystyle\hat{\bf P}=\cos(\theta_{\alpha})\hat{\bf P}_{1}+\sin(\theta_{\alpha})\hat{\bf P}_{2}\cdot e^{i\theta_{\beta}},

(60)

here, $\theta_{\alpha}$ denotes the angle between ${\bf P}$ and $\hat{{\bf P}}_{1}$ , while $\theta_{\beta}$ represents the absolute phase. In quantum mechanics, the photon polarization state corresponding to P can be described by the density matrix:

\displaystyle\rho=\frac{1}{2}\left(1+{\bm{\xi}\cdot\bm{\sigma}}\right)=\frac{1}{2}\left(\begin{matrix}1+\xi_{3}&\xi_{1}-i\xi_{2}\\ \xi_{1}+i\xi_{2}&1-\xi_{3}\end{matrix}\right)

(61)

where $\bm{\sigma}$ represents the Pauli matrix, and $\bm{\xi}=(\xi_{1},\xi_{2},\xi_{3})$ denotes the Stokes parameters, with $\xi_{1}=\sin(2\theta_{\alpha})\cos(\theta_{\beta})$ , $\xi_{2}=\sin(2\theta_{\alpha})\sin(\theta_{\beta})$ , $\xi_{3}=\cos(2\theta_{\alpha})$ .

The calculation of the probability of pair creation requires the transformation of the Stokes parameters from the initial frame of the photon ( $\hat{{\bf P}}_{1}$ , $\hat{{\bf P}}_{2}$ , $\hat{\bf n}$ ) to the frame of pair production ( $\hat{{\bf P}}_{1}^{\prime}$ , $\hat{{\bf P}}_{2}^{\prime}$ , $\hat{\bf n}$ ). The vector $\hat{{\bf P}}_{1}^{\prime}$ is given by [ ${\bf E}-\hat{\bf n}\cdot(\hat{\bf n}\cdot{\bf E})+\hat{\bf n}\times{\bf B}$ ]/ $|{\bf E}-\hat{\bf n}\cdot(\hat{\bf n}\cdot{\bf E})+\hat{\bf n}\times{\bf B}|$ , and the vector $\hat{{\bf P}}_{2}^{\prime}$ is obtained by taking the cross product of $\hat{\bf n}$ and $\hat{{\bf P}}_{1}^{\prime}$ . Here, $\hat{\bf n}$ represents the direction of propagation of the photon, and ${\bf E}$ and ${\bf B}$ denote the electric and magnetic fields, respectively. Two groups of polarization vector are connected via a rotation of angle $\psi$ :

	$\displaystyle\hat{\bf P}_{1}^{{}^{\prime}}$	$\displaystyle=$	$\displaystyle\hat{{\bf P}}_{1}{\rm cos}(\psi)+\hat{{\bf P}}_{2}{\rm sin}(\psi),$		(62)
	$\displaystyle\hat{\bf P}_{2}^{{}^{\prime}}$	$\displaystyle=$	$\displaystyle-\hat{{\bf P}}_{1}{\rm sin}(\psi)+\hat{{\bf P}}_{2}{\rm cos}(\psi).$		(63)

Thus, the Stokes parameters with respect to the vectors $\hat{\bf P}_{1}^{{}^{\prime}}$ , $\hat{\bf P}_{2}^{{}^{\prime}}$ and $\hat{\bf n}$ are the follows:

$\displaystyle\xi_{1}^{{}^{\prime}}$	$\displaystyle=$	$\displaystyle\xi_{1}{\rm cos}(2\psi)-\xi_{3}{\rm sin}(2\psi),$
$\displaystyle\xi_{2}^{{}^{\prime}}$	$\displaystyle=$	$\displaystyle\xi_{2},$
$\displaystyle\xi_{3}^{{}^{\prime}}$	$\displaystyle=$	$\displaystyle\xi_{1}{\rm sin}(2\psi)+\xi_{3}{\rm cos}(2\psi),$	(64)

which is equivalent to a rotation Wistisen and Uggerhøj (2013); Bragin et al. (2017):

\left(\begin{array}[]{l}\xi^{\prime}_{1}\\ \xi^{\prime}_{2}\\ \xi^{\prime}_{3}\end{array}\right)=\left(\begin{array}[]{ccc}\cos 2\psi&0&-\sin 2\psi\\ 0&1&0\\ \sin 2\psi&0&\cos 2\psi\end{array}\right)\left(\begin{array}[]{l}\xi_{1}\\ \xi_{2}\\ \xi_{3}\end{array}\right)\equiv{\mathrm{ROT}}(\psi)\cdot\bm{\xi}.

(65)

III.4 Nonlinear Breit-Wheeler pair production

When the energy of a photon exceeds the rest mass of an electron-positron pair, i.e., $\omega_{\gamma}\geq 2m_{e}c^{2}$ , and is subjected to an ultraintense field of $a_{0}\gg 1$ , the related nonlinear quantum parameter $\chi_{\gamma}$ can reach unity. Here, $\chi_{\gamma}\equiv\frac{e\hbar}{m^{2}c^{3}}\sqrt{|F^{\mu\nu}k_{\nu}|^{2}}$ and is approximately equal to $2\omega_{\gamma}a_{0}/m^{2}_{e}$ 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 ( $\omega_{\gamma}+n\omega_{L}\rightarrow e^{+}+e^{-}$ ) Bell and Kirk (2008). Refs. Chen et al. (2019); Wan et al. (2020b); Xue et al. (2022); Li et al. (2022); Chen et al. (2022) proposed the spin- and polarization-resolved NBW MC method, and we followed the detailed methods in Ref. Xue et al. (2022).

III.4.1 NBW Probability

The polarization-resolved NBW probability rate with dependence on the positron energy is given by

\displaystyle\frac{{\rm d}^{2}W^{\pm}_{\rm pair}}{{\rm d}\varepsilon_{+}{\rm d}t}=\frac{1}{2}(G_{0}+\xi_{1}G_{1}+\xi_{2}G_{2}+\xi_{3}G_{3}),

(66)

where the polarization-independent term $G_{0}$ and polarization-related terms $G_{1,2,3}$ are given by

$\displaystyle G_{0}$	$\displaystyle=$	$\displaystyle\frac{W_{0}}{2}\left\{{\rm IntK}_{\frac{1}{3}}(\rho)+{\frac{\varepsilon_{-}^{2}+\varepsilon_{+}^{2}}{\varepsilon_{-}\varepsilon_{+}}}{\rm K}_{\frac{2}{3}}(\rho)+\left[{\rm IntK}_{\frac{1}{3}}(\rho)-2{\rm K}_{\frac{2}{3}}(\rho)\right]\left({\bf S}_{-}\cdot{\bf S}_{+}\right)+\right.$	(67)
		$\displaystyle\left.{\rm K}_{\frac{1}{3}}(\rho)\left[-\frac{\varepsilon_{\gamma}}{\varepsilon_{+}}\left({\bf S}_{+}\cdot\hat{{\bf b}}_{+}\right)+\frac{\varepsilon_{\gamma}}{\varepsilon_{-}}\left({\bf S}_{-}\cdot\hat{{\bf b}}_{+}\right)\right]+\left[\frac{\varepsilon_{-}^{2}+\varepsilon_{+}^{2}}{\varepsilon_{-}\varepsilon_{+}}{\rm IntK}_{\frac{1}{3}}(\rho)\right.\right.$
		$\displaystyle\left.\left.-\frac{(\varepsilon_{+}-\varepsilon_{-})^{2}}{\varepsilon_{-}\varepsilon_{+}}{\rm K}_{\frac{2}{3}}(\rho)\right]({\bf S}_{+}\cdot\hat{\bf v}_{+})({\bf S}_{-}\cdot\hat{\bf v}_{+})\right\},$
$\displaystyle G_{1}$	$\displaystyle=$	$\displaystyle\frac{W_{0}}{2}\left\{{\rm K}_{\frac{1}{3}}(\rho)\left[-\frac{\varepsilon_{\gamma}}{\varepsilon_{-}}({\bf S}_{+}\cdot\hat{{\bf a}}_{+})+\frac{\varepsilon_{\gamma}}{\varepsilon_{+}}({\bf S}_{-}\cdot\hat{{\bf a}}_{+})\right]+\frac{\varepsilon_{+}^{2}-\varepsilon_{-}^{2}}{2\varepsilon_{-}\varepsilon_{+}}{\rm K}_{\frac{2}{3}}(\rho)({\bf S}_{-}\times{\bf S}_{+})\cdot\hat{\bf v}_{+}\right.$	(68)
		$\displaystyle\left.-\frac{\varepsilon_{\gamma}^{2}}{2\varepsilon_{-}\varepsilon_{+}}{\rm IntK}_{\frac{1}{3}}(\rho)\left[({\bf S}_{+}\cdot\hat{{\bf a}})({\bf S}_{-}\cdot\hat{{\bf b}})+({\bf S}_{-}\cdot\hat{{\bf a}}_{+})({\bf S}_{+}\cdot\hat{{\bf b}}_{+})\right]\right\},$	(68)
$\displaystyle G_{2}$	$\displaystyle=$	$\displaystyle\frac{W_{0}}{2}\left\{\frac{\varepsilon_{\gamma}^{2}}{2\varepsilon_{-}\varepsilon_{+}}{\rm K}_{\frac{1}{3}}(\rho)({\bf S}_{-}\times{\bf S}_{+})\cdot\hat{{\bf a}}_{+}+\right.$	(69)
		$\displaystyle\left.\frac{\varepsilon_{+}^{2}-\varepsilon_{-}^{2}}{2\varepsilon_{-}\varepsilon_{+}}{\rm K}_{\frac{1}{3}}(\rho)\left[({\bf S}_{-}\cdot\hat{\bf v}_{+})({\bm{S}}_{+}\cdot\hat{{\bf b}}_{+})+({\bf S}_{+}\cdot\hat{\bf v}_{+})({\bf S}_{-}\cdot\hat{{\bf b}}_{+})\right]+\right.$
		$\displaystyle\left.\left[\frac{\varepsilon_{\gamma}}{\varepsilon_{-}}{\rm IntK}_{\frac{1}{3}}(\rho)-\frac{\varepsilon_{+}^{2}-\varepsilon_{-}^{2}}{\varepsilon_{-}\varepsilon_{+}}{\rm K}_{\frac{2}{3}}(\rho)\right]({\bf S}_{-}\cdot\hat{\bf v}_{+})+\right.$
		$\displaystyle\left.\left[\frac{\varepsilon_{\gamma}}{\varepsilon_{+}}{\rm IntK}_{\frac{1}{3}}(\rho)+\frac{\varepsilon_{+}^{2}-\varepsilon_{-}^{2}}{\varepsilon_{-}\varepsilon_{+}}{\rm K}_{\frac{2}{3}}(\rho)\right]({\bf S}_{+}\cdot\hat{\bf v}_{+})\right\},$
$\displaystyle G_{3}$	$\displaystyle=$	$\displaystyle\frac{W_{0}}{2}\left\{-{\rm K}_{\frac{2}{3}}(\rho)+\frac{\varepsilon_{-}^{2}+\varepsilon_{+}^{2}}{2\varepsilon_{-}\varepsilon_{+}}{\rm K}_{\frac{2}{3}}(\rho)({\bf S}_{-}\cdot{\bf S}_{+})-\right.$	(70)
		$\displaystyle\left.{\rm K}_{\frac{1}{3}}(\rho)\left[\frac{\varepsilon_{\gamma}}{\varepsilon_{+}}({\bf S}_{-}\cdot\hat{{\bf b}}_{+})-\frac{\varepsilon_{\gamma}}{\varepsilon_{-}}({\bf S}_{+}\cdot\hat{{\bf b}}_{+})\right]+\right.$
		$\displaystyle\left.\frac{\varepsilon_{\gamma}^{2}}{2\varepsilon_{-}\varepsilon_{+}}{\rm IntK}_{\frac{1}{3}}(\rho)\bigg{[}({\bf S}_{+}\cdot\hat{{\bf b}}_{+})({\bf S}_{-}\cdot\hat{{\bf b}}_{+})\right.-$
		$\displaystyle\left.({\bf S}_{+}\cdot\hat{{\bf a}}_{+})({\bf S}_{-}\cdot\hat{{\bf a}}_{+})\bigg{]}-\frac{(\varepsilon_{+}-\varepsilon_{-})^{2}}{2\varepsilon_{-}\varepsilon_{+}}{\rm K}_{\frac{2}{3}}(\rho)({\bf S}_{+}\cdot\hat{\bf v}_{+})({\bf S}_{-}\cdot\hat{\bf v}_{+})\right\},$

where $W_{0}=\alpha/\left(\sqrt{3}\pi\omega^{\prime 2}_{\gamma}\right)$ , $\omega^{\prime}_{\gamma}=\varepsilon_{\gamma}/m_{e}c^{2}$ , $\rho=2\varepsilon_{\gamma}^{2}/\left(3\chi_{\gamma}\varepsilon_{-}\varepsilon_{+}\right)=2/\left[3\delta(1-\delta)\right]$ , $\delta=\varepsilon_{+}/\varepsilon_{\gamma}$ , ${\rm IntK}_{\frac{1}{3}}(\rho)\equiv\int_{\rho}^{\infty}{\rm d}z{\rm K}_{\frac{1}{3}}(z)$ , ${\rm K}_{n}$ is the $n$ -order modified Bessel function of the second kind, $\alpha$ the fine structure constant, $\varepsilon_{\gamma}$ , $\varepsilon_{-}$ and $\varepsilon_{+}$ the energies of parent photon, created electron and positron, respectively, $\hat{\bf v}_{+}={\bf v}_{+}/|{\bf v}_{+}|$ with the positron velocity ${\bf v}_{+}$ , $\hat{{\bf a}}_{+}={\bf a}_{+}/|{\bf a}_{+}|$ with the positron acceleration ${\bf a}_{+}$ in the rest frame of positron, $\hat{{\bf b}}_{+}={\bf v}_{+}\times{\bf a}_{+}/|{\bf v}_{+}\times{\bf a}_{+}|$ , $\xi_{1}$ , $\xi_{2}$ and $\xi_{3}$ are the Stokes parameters of $\gamma$ photon, and ${\bf S}_{+}$ ( ${\bf S}_{-}$ ) denotes the positron (electron) spin vector. Note that the Stokes parameters must be transformed from the photon initial frame ( $\hat{{\bf P}}_{1}$ , $\hat{{\bf P}}_{2}$ , $\hat{\bf n}$ ) to the pair production frame ( $\hat{{\bf P}}_{1}^{\prime}$ , $\hat{{\bf P}}_{2}^{\prime}$ , $\hat{\bf n}$ ); see transformations of the Stokes parameters in Sec. III.3.2.

Summing over the electron spin, the pair production depending on the positron spin ${\bf S}_{+}$ and the photon polarization $\bm{\xi}$ is obtained:

	$\displaystyle\frac{{\rm d}^{2}W_{\rm pair}^{+}}{{\rm d}\varepsilon_{+}{\rm d}t}=W_{0}\bigg{\{}{\rm IntK}_{\frac{1}{3}}(\rho)+\frac{\varepsilon_{-}^{2}+\varepsilon_{+}^{2}}{\varepsilon_{-}\varepsilon_{+}}{\rm K}_{\frac{2}{3}}(\rho)-\frac{\varepsilon_{\gamma}}{\varepsilon_{+}}{\rm K}_{\frac{1}{3}}(\rho)({\bf S}_{+}\cdot\hat{{\bf b}}_{+})$
	$\displaystyle-\xi_{1}\bigg{[}\frac{\varepsilon_{\gamma}}{\varepsilon_{-}}{\rm K}_{\frac{1}{3}}(\rho)({\bf S}_{+}\cdot\hat{{\bf a}}_{+})\bigg{]}+\xi_{2}\left[\frac{\varepsilon_{+}^{2}-\varepsilon_{-}^{2}}{\varepsilon_{-}\varepsilon_{+}}{\rm K}_{\frac{2}{3}}(\rho)+\frac{\varepsilon_{\gamma}}{\varepsilon_{+}}{\rm IntK}_{\frac{1}{3}}(\rho)\right]\times$
	$\displaystyle\left.({\bf S}_{+}\cdot\hat{\bf v}_{+})\right]-\xi_{3}\left[{\rm K}_{\frac{2}{3}}(\rho)-\frac{\varepsilon_{\gamma}}{\varepsilon_{-}}{\rm K}_{\frac{1}{3}}(\rho)({\bf S}_{+}\cdot\hat{{\bf b}}_{+})\right]\bigg{\}}.$		(71)

It can be rewritten as:

\displaystyle\frac{{\rm d}^{2}W_{\rm pair}^{+}}{{\rm d}\varepsilon_{+}{\rm d}t}

\displaystyle=

\displaystyle W_{0}(C+{\bf S}_{+}\cdot{\bf D}),

(72)

where

$\displaystyle C$	$\displaystyle=$	$\displaystyle{\rm IntK}_{\frac{1}{3}}(\rho)+\frac{\varepsilon_{-}^{2}+\varepsilon_{+}^{2}}{\varepsilon_{-}\varepsilon_{+}}{\rm K}_{\frac{2}{3}}(\rho)-\xi_{3}{\rm K}_{\frac{2}{3}}(\rho),$	(73)
$\displaystyle{\bf D}$	$\displaystyle=$	$\displaystyle-\left(\frac{\varepsilon_{\gamma}}{\varepsilon_{+}}-\xi_{3}\frac{\varepsilon_{\gamma}}{\varepsilon_{-}}\right){\rm K}_{\frac{1}{3}}(\rho)\hat{{\bf b}}_{+}-\xi_{1}\frac{\varepsilon_{\gamma}}{\varepsilon_{-}}{\rm K}_{\frac{1}{3}}(\rho)\hat{{\bf a}}_{+}+$	(74)
		$\displaystyle\xi_{2}\bigg{[}\frac{\varepsilon_{+}^{2}-\varepsilon_{-}^{2}}{\varepsilon_{-}\varepsilon_{+}}{\rm K}_{\frac{2}{3}}(\rho)+\frac{\varepsilon_{\gamma}}{\varepsilon_{+}}{\rm IntK}_{\frac{1}{3}}(\rho)\bigg{]}\hat{\bf v}_{+}.$	(74)

When a photon decays to 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 ${\bf S}_{+}^{(\varepsilon_{+})}={\bf D}/C$ .

Similarly, summing over the positron spin, the pair production probability depending on the electron spin ${\bf S}_{-}$ and the photon polarization is obtained:

	$\displaystyle\frac{{\rm d}^{2}W_{\rm pair}^{-}}{{\rm d}\varepsilon_{+}{\rm d}t}=W_{0}(C+{\bf S}_{-}\cdot{\bf D}^{{}^{\prime}}),$		(75)
	$\displaystyle{\bf D}^{{}^{\prime}}=\left(\frac{\varepsilon_{\gamma}}{\varepsilon_{-}}-\xi_{3}\frac{\varepsilon_{\gamma}}{\varepsilon_{+}}\right){\rm K}_{\frac{1}{3}}(\rho)\hat{{\bf b}}_{+}+\xi_{1}\frac{\varepsilon_{\gamma}}{\varepsilon_{+}}{\rm K}_{\frac{1}{3}}(\rho)\hat{{\bf a}}_{+}$
	$\displaystyle-\xi_{2}\bigg{[}\frac{\varepsilon_{+}^{2}-\varepsilon_{-}^{2}}{\varepsilon_{-}\varepsilon_{+}}{\rm K}_{\frac{2}{3}}(\rho)-\frac{\varepsilon_{\gamma}}{\varepsilon_{-}}{\rm IntK}_{\frac{1}{3}}(\rho)\bigg{]}\hat{\bf v}_{+}.$		(76)

The pair production probability, relying solely on photon polarization, is determined by summing over both positron and electron spins:

\frac{{\rm d}^{2}W_{\rm pair}}{{\rm d}\varepsilon_{+}{\rm d}t}=2W_{0}\left\{{\rm IntK}_{\frac{1}{3}}(\rho)+\frac{\varepsilon_{-}^{2}+\varepsilon_{+}^{2}}{\varepsilon_{-}\varepsilon_{+}}{\rm K}_{\frac{2}{3}}(\rho)-\xi_{3}{\rm K}_{\frac{2}{3}}(\rho)\right\}.

(77)

III.4.2 MC algorithm

The algorithm for simulating pair creation with polarization is illustrated in Figure. 11. At every simulation step $\Delta t$ , a pair is generated with positron energy $\varepsilon_{+}=r_{1}\varepsilon_{\gamma}$ when the probability density $P\equiv{\rm d^{2}}W_{\rm pair}/{\rm d}\varepsilon_{+}{\rm d}t\cdot\Delta t$ of pair production is greater than or equal to a random number $r_{2}$ within the range [0,1). Here, ${\rm d^{2}}W_{\rm pair}/{\rm d}\varepsilon_{+}{\rm d}t$ is computed using Equation (77). The momentum of the created positron (electron) is parallel to that of the parent photon, and the energy of the electron $\varepsilon_{-}$ is determined as $\varepsilon_{\gamma}-\varepsilon_{+}$ . 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 SQA, where $P_{1,2,3,4}$ is computed from Equation (66). 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 Figure. 12.

III.5 High-energy bremsstrahlung

The high energy bremsstrahlung is another important emission mechanism, which can also be modeled using a MC collision model Wan et al. (2017). The MC collision model was tested using the Geant4 code Agostinelli and et al. (2003), and the results are presented in the following section. The bremsstrahlung emission is described as described by the cross-section in Ref. Tsai (1974)

$\displaystyle\frac{d\sigma_{eZ}}{d\omega}(\omega,y)=$	$\displaystyle\frac{\alpha r_{0}^{2}}{\omega}\{(\frac{4}{3}-\frac{4}{3}y+y^{2})$	(78)
	$\displaystyle\times[Z^{2}(\phi_{1}-\frac{4}{3}\mathrm{ln}Z-4f)+Z(\psi_{1}-\frac{8}{3}\mathrm{ln}Z)]$
	$\displaystyle+\frac{2}{3}(1-y)[Z^{2}(\phi_{1}-\phi_{2})+Z(\psi_{1}-\psi_{2})]\},$

where $y=\hbar\omega/E_{e}$ is the energy ratio of the emitted photon to the incident electron, $r_{0}$ is the classical electron radius, functions $\phi_{1,2}$ and $\psi_{1,2}$ depend on the screening potential by atomic electrons, while the Coulomb correction term is denoted by $f$ . When the atomic number of the target is greater than $5$ , we use Eqs. (3.38-3.41) from Ref. Tsai (1974) to calculate these functions. However, for targets with $Z<5$ , the approximated screen functions are unsuitable and require modification.

The PENELOPE code PEN utilizes another method that involves tabulated data from Ref. Seltzer and Berger (1986). This method transforms the “scaled” bremsstrahlung differential cross-section (DCS) to differential cross-section by using the following equation PEN :

\frac{d\sigma_{br}}{d\omega}=\frac{Z^{2}}{\beta^{2}}\frac{1}{\omega}\chi(Z,E_{e},y),

(79)

where $\beta=v/c$ is the normalized electron velocity. Integrating this expression over the photon frequencies yields a tabulated total cross-section $\sigma_{br}(E_{e},y)$ MC simulation, i.e, the direct sampling method can be used.

The DCS for electron and positron are related by the equation

\frac{d\sigma_{br}^{+}}{d\omega}=F_{p}(Z,E_{e})\frac{d\sigma_{br}^{-}}{d\omega},

(80)

where $F_{p}(Z,E_{e})$ is an analytical approximation factor that can be found in Ref. PEN . This reference demonstrates a high level of accuracy, with a difference of only approximately $0.5\%$ compared to Ref. Kim et al. (1986).

The bremsstrahlung implementation is based on a 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^{\Delta s/\lambda}$ , where $\Delta s=v\Delta t$ , $v=|\vec{v}|$ is the incident particle velocity, $\Delta t$ is the time interval, $\lambda=1/n\sigma(E_{e})$ , $n$ represents the target particle density and $\sigma(E_{e})$ is the total cross-section, respectively. A random number $r_{1}$ is then generated and compared to $P_{br}$ . If $r_{1}<P_{br}$ , 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 $\sigma_{br}(E_{e})$ to obtain the energy ratio $y$ through $\sigma(y,E_{e})=\sigma(E_{e})r_{2}$ . Finally, a photon with energy $\hbar\omega=E_{e}y$ and momentum direction $\vec{k}/|\vec{k}|=\mathbf{v}/|\mathbf{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 PEN . The implementation of Bethe-Heitler pair production follows a similar process.

The implementation of Bremsstrahlung emission was tested using the Geant4 software Agostinelli and et al. (2003), which is widely used for modeling high-energy particle scattering with detectors. In this study, we utilized electron bunches of 1 GeV and 100 MeV with $10^{5}$ primaries to collide with a 5 mm Au target with $Z=79$ , $\rho=19.3$ $\mathrm{g/cm^{3}}$ and a 5 mm Al target with $Z=13$ , $\rho=2.7$ $\mathrm{g/cm^{3}}$ . We disabled the field updater and weighting procedure in the PIC code, and only enabled 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 $\mathrm{MeV}$ electron bunch normally incident onto the aluminum and gold slabs, while similar distributions for a 1 $\mathrm{GeV}$ electron bunch are shown in Figure. 14.

III.6 Vacuum birefringence

Another important process for polarized photons in ultraintense laser matter interactions is vacuum birefringence (VB), in addition to the NBW processes. In this paper, we utilize Eq.(4.26) from Ref. Shore (2007) to calculate the refractive index $n$ for a photon with arbitrary energy $\omega$ (wavelength $\lambda$ ) in a constant weak EM field [ $|E|(|B|)\ll E_{cr}$ ]. We include the electric field and assume relativistic units $c=\hbar=1$ . The resulting expression is:

n\approx 1-\frac{\alpha\chi_{\gamma}^{2}m^{2}}{16\pi\omega^{2}}\int_{-1}^{1}d\upsilon(1-\upsilon^{2})\left\{\begin{matrix}\frac{1}{2}(1+\frac{1}{3}\upsilon^{2})\\ 1-\frac{1}{3}\upsilon^{2}\end{matrix}\right\}\left[\pi x^{4/3}\mathrm{Gi}^{\prime}(x^{2/3})-i\frac{x^{2}}{\sqrt{3}}\mathrm{K}_{2/3}\left(\frac{2}{3}x\right)\right],

(81)

where $\alpha$ is the fine structure constant, $m$ is the electron mass, $\chi_{\gamma}$ is the nonlinear quantum parameter as defined before, $x=\frac{4}{(1-\upsilon^{2})\chi_{\gamma}}$ , ${\rm Gi}^{\prime}(x)$ is the derivative of the Scorer’s function, and ${\rm K}_{n}(x)$ is the $n$ th-order modified Bessel function of the second kind [38]. ${\bf E}_{\rm red,\perp}={\bf E}_{\perp}+\hat{k}\times{\bf B}_{\perp}$ is the transverse reduced field (acceleration field for electrons). The first and second columns correspond to the eigenmodes parallel and perpendicular to the reduced field, respectively. By extracting a factor of

\mathcal{D}=\frac{\alpha}{90\pi}\left(\frac{e|{\bf E}_{\rm red\perp}|}{m^{2}}\right)^{2}\equiv\frac{\alpha}{90\pi}\frac{\chi_{\gamma}^{2}}{\omega^{2}/m^{2}}.

Eq. (81) arrives at

	$\displaystyle\mathrm{Re}(n)$	$\displaystyle=$	$\displaystyle 1-\frac{45}{4}\mathcal{D}\int_{0}^{1}d\upsilon(1-\upsilon^{2})\left\{\begin{matrix}\frac{1}{2}(1+\frac{1}{3}\upsilon^{2})\\ 1-\frac{1}{3}\upsilon^{2}\end{matrix}\right\}\left[\pi x^{4/3}\mathrm{Gi}^{\prime}(x^{2/3})\right],$		(82)
	$\displaystyle\mathrm{Im}(n)$	$\displaystyle=$	$\displaystyle\frac{45}{4}\mathcal{D}\int_{0}^{1}d\upsilon(1-\upsilon^{2})\left\{\begin{matrix}\frac{1}{2}(1+\frac{1}{3}\upsilon^{2})\\ 1-\frac{1}{3}\upsilon^{2}\end{matrix}\right\}\left[\frac{x^{2}}{\sqrt{3}}\mathrm{K}_{2/3}\left(\frac{2}{3}x\right)\right].$		(83)

In the weak field limit of $\chi_{\gamma}\ll 1$ , the imaginary part associated with pair production, is negligible. Now we define

\begin{array}[]{l}M(\chi_{\gamma})=-\frac{45}{4}\int_{0}^{1}d\upsilon(1-\upsilon^{2})\left\{\begin{matrix}\frac{1}{2}(1+\frac{1}{3}\upsilon^{2})\\ 1-\frac{1}{3}\upsilon^{2}\end{matrix}\right\}\left[\pi x^{4/3}\mathrm{Gi}^{\prime}(x^{2/3})\right],\leavevmode\nobreak\ \mathrm{yielding}\\ \mathrm{Re}(n)=1+M(\chi_{\gamma})\mathcal{D}\equiv 1+M(\chi_{\gamma})\frac{\alpha}{90\pi}\frac{\chi_{\gamma}^{2}}{\omega^{2}/m^{2}}.\end{array}

(84)

The numerical results of $M(\chi_{\gamma})$ and comparisons with the low-energy-limit ( $\omega_{\gamma}\ll m$ ) constants are given in Figure. 15.

While, in the limit of $\chi_{\gamma}\ll 1$ , the real part simplifies to

\displaystyle\mathrm{Re}(n)=1+\mathcal{D}\left\{\begin{matrix}4_{+}\\ 7_{-}\end{matrix}\right\},

(85)

and can be used to simulate the vacuum birefringence (VB) effect with good accuracy for $\chi_{\gamma}\ll 1$ . Note that these results are identical to those in Refs. Shore (2007); Karbstein (2013); Dinu et al. (2014). For large $\chi_{\gamma}$ , two interpolated refractive indexes are used.

The phase retardation between two orthogonal components is given by $\delta\phi=\phi_{+}-\phi_{-}=\Delta n\frac{2\pi l}{\lambda}=-3\mathcal{D}\frac{2\pi l}{\lambda}$ with $l$ denotes the propagation length, and the VB effect is equivalent to a rotation of the Stokes parameters:

\left(\begin{array}[]{l}\xi^{\prime}_{1}\\ \xi^{\prime}_{2}\\ \xi^{\prime}_{3}\end{array}\right)=\left(\begin{array}[]{ccc}\cos\delta\phi&-\sin\delta\phi&0\\ \sin\delta\phi&\cos\delta\phi&0\\ 0&0&1\end{array}\right)\left(\begin{array}[]{l}\xi_{1}\\ \xi_{2}\\ \xi_{3}\end{array}\right)\equiv{\mathrm{QED}}(\delta\phi)\cdot\bm{\xi}.

(86)

The VB effect of the probe photons in the Particle-In-Cell code is simulated with the following algorithms Wan et al. (2022):

1 Initialization part;

2 PIC initialization;

3 foreach photon in photonList do

4 photon.

\bm{\xi}

\bm{\xi}_{0}

;

5 photon.

\hat{a}_{\pm}

= (

\hat{x},\hat{y}

);

7 end foreach

8evolution part;

9 while not final step do

10 do PIC loop …;

11 foreach photon do

12 get

\mathbf{E},\mathbf{B}

;

13 get

\theta

, and

\hat{a}_{+}(\theta)\parallel\mathbf{E}_{\rm red\perp},\hat{a}_{-}(\theta)=\hat{k}\times\hat{a}(\theta)_{+}

;

14 rotate

\bm{\xi}

from

\hat{a}_{\pm}

\hat{a}_{\pm}(\theta)

via Eqs. (65);

15 calculate new

\bm{\xi}

via Eq. (86) ;

16 update the polarization basis: photon.

\hat{a}_{\pm}=\hat{a}_{\pm}(\theta)

17 end foreach

19 end while

20Post-processing;

21 select a detection plane (polarization basis), for instance $\hat{x},\hat{y}$ ;

22 foreach $x,y$ in the detector do

23 foreach photon in this area do

24 rotate

\bm{\xi}

from

\hat{a}_{\pm}

(\hat{x},\hat{y})

via Eqs. (65);

25 average all

\bm{\xi}

26 end foreach

28 end foreach

Algorithm 1 VB effect in SLIPs

See an example of VB effect in Figure. 16.

IV Framework of SLIPs

These physical processes have been incorporated into the 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 Figure. 17, SLIPs utilize 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 is employed to generate a SimBox that comprises all ParticleList and Fields, and define the FieldSolver, EOMSolver and initialize QED processes.

The internal data structure of SLIPs is constructed using the open-source numerical library, Armadillo C++ Sanderson and Curtin (2016, 2018). String expressions are parsed using the Exprtk library exp . The data is then dumped using serial-hdf5 and merged with external Python scripts to remove ghost cells.

The spin-resolved processes, i.e., tagged as Spin-QED in the diagram in Figure. 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: $\texttt{Lorentz \& T-BMT}\Rightarrow\texttt{Radiative correction}\Rightarrow\texttt{VB}\Rightarrow\texttt{NBW}\Rightarrow\texttt{NCS \& Bremss}$ .

V Polarized particles simulations

In this section, we present known results that were calculated from the single-particle mode using the 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. Li et al. (2019) and Wan et al. (2020b).

V.1 Polarized electron/positron simulation

To simulate the generation of spin-polarized electrons, we utilized an elliptically polarized laser with an intensity of $a_{0}=30$ , a wavelength of $\lambda_{0}=1\mathrm{\mu m}$ , and an ellipticity of $a_{y,0}/a_{x,0}=3\%$ . This laser was directed towards an ultrarelativistic electron bunch with an energy of 10 GeV, which was produced through laser-wakefield acceleration. The resulting polarized electrons are depicted in Figure. 19, and shows good agreement with the previously published results in Ref. Li et al. (2022).

V.2 Polarized $\gamma$ -photons via NCS

The polarization state of emitted photons can be determined in the spin/polarization-resolved NCS. Here, follow Ref. Li et al. (2022), we utilized a linearly polarized (LP) laser to collide with an unpolarized electron bunch to generate LP $\gamma$ -photons. Additionally, we used an LP laser to collide with a longitudinally polarized electron bunch to generate circularly polarized (CP) $\gamma$ -photons, which were also observed in a previous study (Ref. Li et al. (2020)). The final polarization states of LP and CP $\gamma$ -photons are presented in Figs. 20 and 21, respectively.

V.3 Laser-plasma interactions

Finally, we present a simulation result demonstrating the interaction between an ultra-intense laser with a normalized intensity of $a_{0}=1000$ and a fully ionized $2\mathrm{\mu m}$ -thick aluminum target. Note, this configuration, previously examined in Ref. Ridgers et al. (2012) with a thickness of $1\mathrm{\mu m}$ , employs a thicker target in the present study to enhance the SF-QED processes. When the laser is directed towards 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 $\gamma$ photons via NCS. The EM field distribution and number density of target electrons, NBW positrons, and NCS $\gamma$ photons are shown in Figure. 22, both of which show good consistency with Ref. Ridgers et al. (2012). 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, $\bm{e}_{1}\equiv\bm{\beta}\times\bm{\dot{\beta}}$ and $\bm{e}_{2}\equiv\hat{\bf n}\times\bm{e}_{1}$ , where $\hat{\bf n}$ denotes the momentum direction of the photon. The polarization angle-dependence observed in this study is consistent with prior literature. However, the average linear polarization degree is approximately 60% ( $\bar{\xi_{3}}\simeq 0.6$ ), as illustrated in Figs. 23(b) and (d). Notably, low-energy photons contribute primarily to the polarization, as demonstrated in Figs. 23(a) and (c). Additionally, during the subsequent NBW process, the self-generated strong magnetic field couples with the laser field dominate the positrons’ SQA. As a result, the positrons’ polarization is aligned with the $z$ direction, contingent on their momentum direction, as shown in Figure. 24. These findings constitute a novel contribution to the investigation of polarization-resolved laser-plasma interactions.

VI Outlook

Computer simulation techniques for laser and plasma interactions are constantly evolving, not only in 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 techniques not only provides opportunities for experimental verification of SF-QED processes in the high-energy density regime (which serves as a micro-astrophysics lab) but also presents challenges in 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.

VII Acknowlegement

The work is supported by the National Natural Science Foundation of China (Grants No. 12275209, 12022506, and U2267204), Open Foundation of Key Laboratory of High Power Laser and Physics, Chinese Academy of Sciences (SGKF202101), the Foundation of Science and Technology on Plasma Physics Laboratory (No. JCKYS2021212008), and the Shaanxi Fundamental Science Research Project for Mathematics and Physics (Grant No. 22JSY014).

References

Piazza et al. (2012) A. D. Piazza, C. Müller, K. Z. Hatsagortsyan, and C. H. Keitel, “Extremely high-intensity laser interactions with fundamental quantum systems,” Rev. Mod. Phys. 84, 1177–1228 (2012).
Bell and Kirk (2008) A. R. Bell and J. G. Kirk, “Possibility of prolific pair production with high-power lasers,” Phys. Rev. Lett. 101, 200403 (2008).
Gonsalves et al. (2019) A. J. Gonsalves, K. Nakamura, J. Daniels, C. Benedetti, C. Pieronek, T. C. H. de Raadt, S. Steinke, J. H. Bin, S. S. Bulanov, J. van Tilborg, C. G. R. Geddes, C. B. Schroeder, C. Tóth, E. Esarey, K. Swanson, L. Fan-Chiang, G. Bagdasarov, N. Bobrova, V. Gasilov, G. Korn, P. Sasorov, and W. P. Leemans, “Petawatt laser guiding and electron beam acceleration to 8 gev in a laser-heated capillary discharge waveguide,” Phys. Rev. Lett. 122, 084801 (2019).
Esarey, Schroeder, and Leemans (2009) E. Esarey, C. B. Schroeder, and W. P. Leemans, “Physics of laser-driven plasma-based electron accelerators,” Rev. Mod. Phys. 81, 1229–1285 (2009).
Ritus (1985) V. I. Ritus, “Quantum effects of the interaction of elementary particles with an intense electromagnetic field,” J. Russ. Laser Res. 6, 497–617 (1985).
Baier, Katkov, and Strakhovenko (1998) V. N. Baier, V. M. Katkov, and V. M. Strakhovenko, Electromagnetic Processes at High Energies in Oriented Single Crystals (World Scientific, 1998).
Jackson (2021) J. D. Jackson, Classical electrodynamics (John Wiley & Sons, 2021).
(8) I. M. T. A. A. Sokolov, Radiation from Relativistic Electrons, American Institute of Physics Translation Series (American Institute of Physics).
Del Sorbo et al. (2017) D. Del Sorbo, D. Seipt, T. G. Blackburn, A. G. R. Thomas, C. D. Murphy, J. G. Kirk, and C. P. Ridgers, “Spin polarization of electrons by ultraintense lasers,” Phys. Rev. A 96, 043407 (2017).
Li et al. (2019) Y.-F. Li, R. Shaisultanov, K. Z. Hatsagortsyan, F. Wan, C. H. Keitel, and J.-X. Li, “Ultrarelativistic electron-beam polarization in single-shot interaction with an ultraintense laser pulse,” Phys. Rev. Lett. 122, 154801 (2019).
King and Tang (2020) B. King and S. Tang, “Nonlinear compton scattering of polarized photons in plane-wave backgrounds,” Phys. Rev. A 102, 022809 (2020).
Li et al. (2020) Y.-F. Li, R. Shaisultanov, Y.-Y. Chen, F. Wan, K. Z. Hatsagortsyan, C. H. Keitel, and J.-X. Li, “Polarized Ultrashort Brilliant Multi-GeV $\gamma$ -Rays via Single-Shot Laser-Electron Interaction,” Phys. Rev. Lett. 124, 014801 (2020).
Gong, Hatsagortsyan, and Keitel (2021) Z. Gong, K. Z. Hatsagortsyan, and C. H. Keitel, “Retrieving transient magnetic fields of ultrarelativistic laser plasma via ejected electron polarization,” Phys. Rev. Lett. 127, 165002 (2021).
Wan et al. (2020a) F. Wan, Y. Wang, R.-T. Guo, Y.-Y. Chen, R. Shaisultanov, Z.-F. Xu, K. Z. Hatsagortsyan, C. H. Keitel, and J.-X. Li, “High-energy $\gamma$ -photon polarization in nonlinear breit-wheeler pair production and $\gamma$ polarimetry,” Phys. Rev. Research 2, 032049 (2020a).
Xue et al. (2020) K. Xue, Z.-K. Dou, F. Wan, T.-P. Yu, W.-M. Wang, J.-R. Ren, Q. Zhao, Y.-T. Zhao, Z.-F. Xu, and J.-X. Li, “Generation of highly-polarized high-energy brilliant $\gamma$ -rays via laser-plasma interaction,” Matter Radiat. Extremes 5, 054402 (2020).
Tang, King, and Hu (2020) S. Tang, B. King, and H. Hu, “Highly polarised gamma photons from electron-laser collisions,” Phys. Lett. B 809, 135701 (2020).
Song, Wang, and Li (2022) H.-H. Song, W.-M. Wang, and Y.-T. Li, “Dense polarized positrons from laser-irradiated foil targets in the QED regime,” Phys. Rev. Lett. 129, 035001 (2022).
Arber et al. (2015) T. D. Arber, K. Bennett, C. S. Brady, A. Lawrence-Douglas, M. G. Ramsay, N. J. Sircombe, P. Gillies, R. G. Evans, H. Schmitz, A. R. Bell, and C. P. Ridgers, “Contemporary particle-in-cell approach to laser-plasma modelling,” Plasma Phys. Control. Fusion 57, 113001 (2015).
Birdsall and Langdon (2018) C. Birdsall and A. Langdon, Plasma Physics via Computer Simulation (CRC Press, 2018).
Gonoskov et al. (2015) A. Gonoskov, S. Bastrakov, E. Efimenko, A. Ilderton, M. Marklund, I. Meyerov, A. Muraviev, A. Sergeev, I. Surmin, and E. Wallin, “Extended particle-in-cell schemes for physics in ultrastrong laser fields: Review and developments,” Phys. Rev. E 92, 023305 (2015).
Fonseca et al. (2002) R. A. Fonseca, L. O. Silva, F. S. Tsung, V. K. Decyk, W. Lu, C. Ren, W. B. Mori, S. Deng, S. Lee, T. Katsouleas, and J. C. Adam, “Osiris: A three-dimensional, fully relativistic particle in cell code for modeling plasma based accelerators,” in Computational Science — ICCS 2002, edited by P. M. A. Sloot, A. G. Hoekstra, C. J. K. Tan, and J. J. Dongarra (Springer Berlin Heidelberg, Berlin, Heidelberg, 2002) pp. 342–351.
Burau et al. (2010) H. Burau, R. Widera, W. Hönig, G. Juckeland, A. Debus, T. Kluge, U. Schramm, T. E. Cowan, R. Sauerbrey, and M. Bussmann, “Picongpu: A fully relativistic particle-in-cell code for a gpu cluster,” IEEE Transactions on Plasma Science 38, 2831–2839 (2010).
Derouillat et al. (2018a) J. Derouillat, A. Beck, F. Pérez, T. Vinci, M. Chiaramello, A. Grassi, M. Flé, G. Bouchard, I. Plotnikov, N. Aunai, J. Dargent, C. Riconda, and M. Grech, “Smilei : A collaborative, open-source, multi-purpose particle-in-cell code for plasma simulation,” Comput. Phys. Commun. 222, 351–373 (2018a).
Wu et al. (2020) D. Wu, W. Yu, S. Fritzsche, and X. T. He, “Particle-in-cell simulation method for macroscopic degenerate plasmas,” Phys. Rev. E 102, 033312 (2020).
Li et al. (2022) Y.-F. Li, Y.-Y. Chen, K. Z. Hatsagortsyan, and C. H. Keitel, “Helicity transfer in strong laser fields via the electron anomalous magnetic moment,” Phys. Rev. Lett. 128, 174801 (2022).
Danson et al. (2015) C. Danson, D. Hillier, N. Hopps, and D. Neely, “Petawatt class lasers worldwide,” High Power Laser Sci. Eng. 3, e3 (2015).
(27) “Apollon multi-PW laser Users Facility,” http://www.polytechnique.edu.
Zou et al. (2015) J. Zou, C. Le Blanc, D. Papadopoulos, G. Chériaux, P. Georges, G. Mennerat, F. Druon, L. Lecherbourg, A. Pellegrina, P. Ramirez, and et al., “Design and current progress of the apollon 10 pw project,” High Power Laser Sci. Eng. 3, e2 (2015).
(29) “The Extreme Light Infrastructure (ELI),” http://www.eli-beams.eu/en/facility/lasers/.
Gan et al. (2021) 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, and Z. Xu, “The Shanghai Superintense Ultrafast Laser Facility (SULF) Project,” in Progress in Ultrafast Intense Laser Science XVI, edited by K. Yamanouchi, K. Midorikawa, and L. Roso (Springer International Publishing, Cham, 2021) pp. 199–217.
Buneman (1967) O. Buneman, “Time-reversible difference procedures,” J. Comput. Phys. 1, 517–535 (1967).
Boris and Shanny (1971) J. P. Boris and R. A. Shanny, “Proceedings of the conference on the numerical simulation of plasmas (4th) held at the naval research laboratory, washington, d.c. on 2, 3 november 1970,” (1971), 10.21236/ada023511.
Qin et al. (2013) H. Qin, S. Zhang, J. Xiao, J. Liu, Y. Sun, and W. M. Tang, “Why is boris algorithm so good?” Phys. Plasmas 20, 084503 (2013).
Kane Yee (1966) Kane Yee, “Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propagat. 14, 302–307 (1966).
Esirkepov (2001a) T. Esirkepov, “Exact charge conservation scheme for Particle-in-Cell simulation with an arbitrary form-factor,” Computer Physics Communications 135, 144–153 (2001a).
Esirkepov (2001b) T. Esirkepov, “Exact charge conservation scheme for particle-in-cell simulation with an arbitrary form-factor,” Comput. Phys. Commun. 135, 144–153 (2001b).
Dirac (1938) P. A. M. Dirac, “Classical theory of radiating electrons,” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 167, 148–169 (1938), https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1938.0124 .
Landau and Lifshitz (1999) L. D. Landau and E. Lifshitz, The classical theory of fields (1999).
Ekman, Heinzl, and Ilderton (2021) R. Ekman, T. Heinzl, and A. Ilderton, “Reduction of order, resummation, and radiation reaction,” Phys. Rev. D 104, 036002 (2021).
Ekman (2022) R. Ekman, “Reduction of order and transseries structure of radiation reaction,” Phys. Rev. D 105, 056016 (2022).
Ilderton and Torgrimsson (2013) A. Ilderton and G. Torgrimsson, “Radiation reaction in strong field qed,” Physics Letters B 725, 481–486 (2013).
Seipt and Thomas (2023) D. Seipt and A. G. R. Thomas, “Kinetic theory for spin-polarized relativistic plasmas,” (2023), arXiv:2307.02114 [physics.plasm-ph] .
Neitz and Di Piazza (2014) N. Neitz and A. Di Piazza, “Electron-beam dynamics in a strong laser field including quantum radiation reaction,” Phys. Rev. A 90, 022102 (2014).
Tamburini et al. (2010) M. Tamburini, F. Pegoraro, A. D. Piazza, C. H. Keitel, and A. Macchi, “Radiation reaction effects on radiation pressure acceleration,” New J. Phys. 12, 123005 (2010).
Bulanov et al. (2013a) S. V. Bulanov, T. Z. Esirkepov, M. Kando, J. K. Koga, T. Nakamura, S. S. Bulanov, A. G. Zhidkov, Y. Kato, and G. Korn, “On extreme field limits in high power laser matter interactions: radiation dominant regimes in high intensity electromagnetic wave interaction with electrons,” in High-Power, High-Energy, and High-Intensity Laser Technology; and Research Using Extreme Light: Entering New Frontiers with Petawatt-Class Lasers, edited by J. Hein, G. Korn, and L. O. Silva (SPIE, 2013).
Nikishov and Ritus (1964) A. I. Nikishov and 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).
Sokolov et al. (2009) I. V. Sokolov, N. M. Naumova, J. A. Nees, G. A. Mourou, and V. P. Yanovsky, “Dynamics of emitting electrons in strong laser fields,” Phys. Plasmas 16, 093115 (2009).
Piazza, Hatsagortsyan, and Keitel (2010) A. D. Piazza, K. Z. Hatsagortsyan, and C. H. Keitel, “Quantum radiation reaction effects in multiphoton compton scattering,” Phys. Rev. Lett. 105, 220403 (2010).
Thomas et al. (2012) A. G. R. Thomas, C. P. Ridgers, S. S. Bulanov, B. J. Griffin, 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).
Sokolov et al. (2010) I. V. Sokolov, J. A. Nees, V. P. Yanovsky, N. M. Naumova, 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).
Bulanov et al. (2013b) S. V. Bulanov, T. Z. Esirkepov, M. Kando, J. K. Koga, T. Nakamura, S. S. Bulanov, A. G. Zhidkov, Y. Kato, and G. Korn, “On extreme field limits in high power laser matter interactions: radiation dominant regimes in high intensity electromagnetic wave interaction with electrons,” in High-Power, High-Energy, and High-Intensity Laser Technology and Research Using Extreme Light: Entering New Frontiers with Petawatt-Class Lasers, edited by J. Hein, G. Korn, and L. O. Silva (SPIE, 2013).
Niel et al. (2018) F. Niel, C. Riconda, F. Amiranoff, R. Duclous, and M. Grech, “From quantum to classical modeling of radiation reaction: A focus on stochasticity effects,” Phys. Rev. E 97, 043209 (2018).
Derouillat et al. (2018b) J. Derouillat, A. Beck, F. Pérez, T. Vinci, M. Chiaramello, A. Grassi, M. Flé, G. Bouchard, I. Plotnikov, N. Aunai, J. Dargent, C. Riconda, and M. Grech, “Smilei : A collaborative, open-source, multi-purpose particle-in-cell code for plasma simulation,” Computer Physics Communications 222, 351–373 (2018b).
Avetissian (2015) H. K. Avetissian, Relativistic Nonlinear Electrodynamics: The QED Vacuum and Matter in Super-Strong Radiation Fields, Vol. 88 (Springer, 2015).
Piazza (2008) A. D. Piazza, “Exact solution of the landau-lifshitz equation in a plane wave,” Letters in Mathematical Physics 83, 305–313 (2008).
Hu and Keitel (1999) S. X. Hu and C. H. Keitel, “Spin signatures in intense laser-ion interaction,” Phys. Rev. Lett. 83, 4709–4712 (1999).
Walser et al. (2002) M. W. Walser, D. J. Urbach, K. Z. Hatsagortsyan, S. X. Hu, and C. H. Keitel, “Spin and radiation in intense laser fields,” Phys. Rev. A 65, 043410 (2002).
Mocken and Keitel (2008) G. R. Mocken and C. H. Keitel, “Fft-split-operator code for solving the dirac equation in 2+1 dimensions,” Computer Physics Communications 178, 868–882 (2008).
Bauke et al. (2014) H. Bauke, S. Ahrens, C. H. Keitel, and R. Grobe, “Relativistic spin operators in various electromagnetic environments,” Phys. Rev. A 89, 052101 (2014).
Wen, Keitel, and Bauke (2017) M. Wen, C. H. Keitel, and H. Bauke, “Spin-one-half particles in strong electromagnetic fields: Spin effects and radiation reaction,” Phys. Rev. A 95, 042102 (2017).
Baĭer (1972) V. N. Baĭer, “Radiative polarization of electron in storage rings,” Soviet Physics Uspekhi 14, 695–714 (1972).
Guo et al. (2020) R.-T. Guo, Y. Wang, R. Shaisultanov, F. Wan, Z.-F. Xu, Y.-Y. Chen, K. Z. Hatsagortsyan, and J.-X. Li, “Stochasticity in radiative polarization of ultrarelativistic electrons in an ultrastrong laser pulse,” Phys. Rev. Research 2, 033483 (2020).
Kirk, Bell, and Arka (2009) J. G. Kirk, A. R. Bell, and I. Arka, “Pair production in counter-propagating laser beams,” Plasma Phys. Controlled Fusion 51, 085008 (2009).
Ridgers et al. (2014) C. Ridgers, J. Kirk, R. Duclous, T. Blackburn, C. Brady, K. Bennett, T. Arber, and A. Bell, “Modelling gamma-ray photon emission and pair production in high-intensity laser–matter interactions,” J. Comput. Phys. 260, 273–285 (2014).
Liu et al. (2022) W.-Y. Liu, K. Xue, F. Wan, M. Chen, J.-X. Li, F. Liu, S.-M. Weng, Z.-M. Sheng, and J. Zhang, “Trapping and acceleration of spin-polarized positrons from $\gamma$ photon splitting in wakefields,” Phys. Rev. Research 4, l022028 (2022).
Green and Harvey (2015) D. Green and 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).
Wan et al. (2020b) F. Wan, R. Shaisultanov, Y.-F. Li, K. Z. Hatsagortsyan, C. H. Keitel, and J.-X. Li, “Ultrarelativistic polarized positron jets via collision of electron and ultraintense laser beams,” Phys. Lett. B 800, 135120 (2020b).
Ritus (1979) V. Ritus, Tr. Fiz. Inst. Akad. Nauk SSSR 111, 6 (1979).
Yokoya et al. (2011) K. Yokoya et al., “Cain, version 2.42,” KEK, Tsukuba, Japan (2011).
Berestetskii, Lifshitz, and Pitaevskii (1982) V. B. Berestetskii, E. M. Lifshitz, and L. P. Pitaevskii, Quantum Electrodynamics: Volume 4, Vol. 4 (Butterworth-Heinemann, 1982).
Wistisen and Uggerhøj (2013) T. N. Wistisen and U. I. Uggerhøj, “Vacuum birefringence by compton backscattering through a strong field,” Phys. Rev. D 88, 053009 (2013).
Bragin et al. (2017) S. Bragin, S. Meuren, C. H. Keitel, and A. D. Piazza, “High-energy vacuum birefringence and dichroism in an ultrastrong laser field,” Phys. Rev. Lett. 119, 250403 (2017).
Chen et al. (2019) Y.-Y. Chen, P.-L. He, R. Shaisultanov, K. Z. Hatsagortsyan, and C. H. Keitel, “Polarized positron beams via intense two-color laser pulses,” Phys. Rev. Lett. 123, 174801 (2019).
Xue et al. (2022) K. Xue, R.-T. Guo, F. Wan, R. Shaisultanov, Y.-Y. Chen, Z.-F. Xu, X.-G. Ren, K. Z. Hatsagortsyan, C. H. Keitel, and J.-X. Li, “Generation of arbitrarily polarized GeV lepton beams via nonlinear breit-wheeler process,” Fundamental Research 2, 539–545 (2022).
Chen et al. (2022) Y.-Y. Chen, K. Z. Hatsagortsyan, C. H. Keitel, and R. Shaisultanov, “Electron spin- and photon polarization-resolved probabilities of strong-field qed processes,” Phys. Rev. D 105, 116013 (2022).
Wan et al. (2017) F. Wan, C. Lv, M. Jia, H. Sang, and B. Xie, “Photon emission by bremsstrahlung and nonlinear compton scattering in the interaction of ultraintense laser with plasmas,” EUROPEAN PHYSICAL JOURNAL D 71 (2017), 10.1140/epjd/e2017-70805-7.
Agostinelli and et al. (2003) S. Agostinelli and et al., “Geant4—a simulation toolkit,” Nucl. Instrum. Meth. A. 506, 250–303 (2003).
Tsai (1974) Y.-S. Tsai, “Pair production and bremsstrahlung of charged leptons,” Rev. Mod. Phys. 46, 815–851 (1974).
(79) “PENELOPE 2018: A code system for Monte Carlo simulation of electron and photon transport,” in Workshop Proceedings, Barcelona, Spain, 28 January – 1 February 2019 (OECD).
Seltzer and Berger (1986) S. M. Seltzer and 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,” Atom. Data Nucl. Data 35, 345–418 (1986).
Kim et al. (1986) L. Kim, R. H. Pratt, S. M. Seltzer, and M. J. Berger, “Ratio of positron to electron bremsstrahlung energy loss: An approximate scaling law,” Phys. Rev. A 33, 3002–3009 (1986).
Shore (2007) G. M. Shore, “Superluminality and UV completion,” Nucl. Phys. B 778, 219–258 (2007).
Karbstein (2013) F. Karbstein, “Photon polarization tensor in a homogeneous magnetic or electric field,” Phys. Rev. D 88, 085033 (2013).
Dinu et al. (2014) V. Dinu, T. Heinzl, A. Ilderton, M. Marklund, and G. Torgrimsson, “Vacuum refractive indices and helicity flip in strong-field QED,” Phys. Rev. D 89, 125003 (2014).
Wan et al. (2022) F. Wan, T. Sun, B.-F. Shen, C. Lv, Q. Zhao, M. Ababekri, Y.-T. Zhao, K. Z. Hatsagortsyan, C. H. Keitel, and J.-X. Li, “Enhanced signature of vacuum birefringence in a plasma wakefield,” (2022), arXiv:2206.10792 [physics.plasm-ph] .
Sanderson and Curtin (2016) C. Sanderson and R. Curtin, “Armadillo: a template-based c++ library for linear algebra,” The Journal of Open Source Software 1, 26 (2016).
Sanderson and Curtin (2018) C. Sanderson and R. Curtin, “A user friendly hybrid sparse matrix class in c++,” in Mathematical Software – ICMS 2018 (Springer International Publishing, 2018) pp. 422–430.
(88) “C++ mathematical expression library,” .
Ridgers et al. (2012) C. P. Ridgers, C. S. Brady, R. Duclous, J. G. Kirk, K. Bennett, T. D. Arber, A. P. L. Robinson, and A. R. Bell, “Dense electron-positron plasmas and ultraintense $\gamma$ rays from laser-irradiated solids,” Phys. Rev. Lett. 108, 165006 (2012).

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

Abstract

I Introduction

II PIC Algorithm

II.1 Particle pushing

II.2 Field solving

II.3 Current deposition

II.4 Force deposition

II.5 Particle shape function

II.6 Time ordering

III QED Algorithm

III.1 Radiative particle pusher

III.1.1 Landau-Lifshitz equation

III.1.2 Modified Landau-Lifshitz equation

III.1.3 Algorithms of the Radiative Pusher

III.2 Spin dynamics

III.2.1 T-BMT equation

III.2.2 Radiative T-BMT equation

III.2.3 Algorithms of simulating the spin precession

III.3 Nonlinear Compton scattering

III.3.1 Spin-resolved/summed nonlinear Compton scattering

III.3.2 Definition and Transformation of Stokes Parameters

III.4 Nonlinear Breit-Wheeler pair production

III.4.1 NBW Probability

III.4.2 MC algorithm

III.5 High-energy bremsstrahlung

III.6 Vacuum birefringence

IV Framework of SLIPs

V Polarized particles simulations

V.1 Polarized electron/positron simulation

V.2 Polarized γ\gamma-photons via NCS

V.3 Laser-plasma interactions

VI Outlook

VII Acknowlegement

References

V.2 Polarized $\gamma$ -photons via NCS