Accurately simulating nine-dimensional phase space of relativistic particles in strong fields
Fei Li, Viktor K. Decyk, Kyle G. Miller, Adam Tableman, Frank S. Tsung, Marija Vranic, Ricardo A. Fonseca, Warren B. Mori
AAccurately simulating nine-dimensional phase space of relativistic particlesin strong fields
Fei Li a, ∗ , Viktor K. Decyk b , Kyle G. Miller b , Adam Tableman b , Frank S. Tsung b , Marija Vranic c ,Ricardo A. Fonseca c,d , Warren B. Mori a,b, ∗ a Department of Electrical Engineering, University of California Los Angeles, Los Angeles, CA 90095, USA b Department of Physics and Astronomy, University of California Los Angeles, Los Angeles, CA 90095, USA c GOLP/Instituto de Plasma e Fus˜ao Nuclear, Instituto Superior T´ecnico, Universidade de Lisboa, Lisbon, Portugal d ISCTE - Instituto Universit´ario de Lisboa, 1649–026, Lisbon, Portugal
Abstract
Next-generation high-power laser systems that can be focused to ultra-high intensities exceeding10 W/cm are enabling new physics regimes and applications. The physics of how these lasersinteract with matter is highly nonlinear, relativistic, and can involve lowest-order quantum effects.The current tool of choice for modeling these interactions is the particle-in-cell (PIC) method.In the presence of strong electromagnetic fields, the motion of charged particles and their spin isaffected by radiation reaction (either the semi-classical or the quantum limit). Standard (PIC)codes usually use Boris or similar operator-splitting methods to advance the particles in standardphase space. These methods have been shown to require very small time steps in the strong-fieldregime in order to obtain accurate results. In addition, some problems require tracking the spinof particles, which creates a nine-dimensional (9D) particle phase space, i.e., ( x , u , s ). Therefore,numerical algorithms that enable high-fidelity modeling of the 9D phase space in the strong-fieldregime (where both the spin and momentum evolution are affected by radiation reaction) aredesired. We present a new particle pusher that works in 9D and 6D phase space (i.e., with andwithout spin) based on analytical rather than leapfrog solutions to the momentum and spin advancefrom the Lorentz force, together with the semi-classical form of radiation reaction in the Landau-Lifshitz equation and spin evolution given by the Bargmann-Michel-Telegdi equation. Analyticalsolutions for the position advance are also obtained, but these are not amenable to the staggeringof space and time in standard PIC codes. These analytical solutions are obtained by assuming alocally uniform and constant electromagnetic field during a time step. The solutions provide the9D phase space advance in terms of a particle’s proper time, and a mapping is used to determinethe proper time step duration for each particle as a function of the lab frame time step. Due to theanalytical integration of particle trajectory and spin orbit, the constraint on the time step neededto resolve trajectories in ultra-high fields can be greatly reduced. The time step required in a PICcode for accurately advancing the fields may provide additional constraints. We present single-particle simulations to show that the proposed particle pusher can greatly improve the accuracy ofparticle trajectories in 6D or 9D phase space for given laser fields. We have implemented the newpusher into the PIC code Osiris . Example simulations show that the proposed pusher providesimprovement for a given time step. A discussion on the numerical efficiency of the proposed pusheris also provided.
Keywords: particle pusher, laser-plasma interaction, radiation reaction, Landau-Lifshitzequation, Bargmann-Michel-Telegdi equation, spin precession, particle-in-cell algorithm1 a r X i v : . [ phy s i c s . c o m p - ph ] A ug . Introduction With the recent advent of petawatt-class lasers and a roadmap for multi-petawatt-class lasersystems [1, 2, 3], laser intensities exceeding 10 W/cm will soon become available. These laserswill open a new door for research avenues in plasma physics, including plasma-based accelera-tion [4, 5, 6, 7] in the strong-field regime, the coupling of laser-plasma interactions and quantumelectrodynamics (QED) [8], and the ability to mimic some astrophysical phenomena (e.g., gamma-ray bursts and supernova explosions) in the laboratory. The physics of how ultra-high-intensitylasers interact with matter is highly nonlinear, relativistic, and involves non-classical processes suchas radiation reaction and quantum effects. Simulations will be a critical partner with experimentsto unravel this physics. The electromagnetic particle-in-cell (PIC) algorithm [9, 10, 11] has beensuccessfully applied to the research of plasma or charged-particle beams interacting with radiationfor nearly half a century. With moderate radiation (laser) parameters, e.g., eA/ ( m e c ) (cid:38)
1, where A is the vector potential of the laser, PIC simulations have proven to be a reliable tool. However,in the strong-field regime where eA/ ( m e c ) (cid:29)
1, accurate modeling becomes much more chal-lenging. Developing high-fidelity PIC simulation algorithms requires a comprehensive and deepunderstanding of each aspect of the numerical algorithm and the physical problem itself. To im-prove the simulation accuracy and reliability, much effort has already been undertaken to mitigatevarious numerical errors; these include improper numerical dispersion, errors to the Lorentz forcefor a relativistic particle interacting with a laser, numerical Cerenkov radiation and an associatedinstability [12, 13, 14, 15, 16], finite-grid instability [17, 18, 19, 20] and spurious fields surroundingrelativistic particles [21].In this article, we address inaccuracies and challenges for the particle pusher used as part ofa PIC code. The pusher has been found to be one of the major factors that prevent high-fidelityPIC simulations in the strong-field regime. Most PIC codes use the standard Boris scheme [22] orone of its variants [23, 24] for the particle push. These later variants correct a shortcoming of theBoris push for a particle moving relativistically where E + v × B ≈
0. In the standard Boris splitalgorithm, the velocity can change even when the Lorentz force vanishes. However, when the fields(forces) are large, these algorithms require small time steps to provide sufficient accuracy.Gordon et al. [25] showed that it is possible to construct an analytic or exact covariant non-splitting pusher. This method assumes the fields (forces) are constant during an interval of theproper time and then advances the particle momentum using analytic solutions. Since this methodpushes particles in the proper time rather than in the observer’s time, it cannot be directly appliedto PIC simulations and can only be used for single-particle tracking. Gordon et al. also discussedhow to include radiation reaction (RR), but used a form of RR that is challenging to incorporate.In some very recent work [26], Gordon and Hafizi propose a more compact in form which they calla special unitary particle pusher. This method provides a method to obtain solutions to all ordersof the time step that maintains Lorentz invariance. They show that with a second-order-accuratemapping from the simulation time step to the proper time step, this pusher can be comparable tothe standard Boris pusher in the push rate. P´etri [27] (who does not seem to be aware of the earlier ∗ Corresponding author
Email addresses: [email protected] (Fei Li), [email protected] (Warren B. Mori)
Preprint submitted to arxiv.org August 6, 2020 ork of Gordon et al.) recently proposed a different implementation of the exact pusher that relieson Lorentz transforming into the particles rest frame and that includes a mapping between theproper and observer time step, allowing the pusher to be applicable for PIC simulations. However,P´etri did not consider RR in his implementation.In strong fields, the motion of charged particles will be significantly impacted by the RR forceand its accompanied energy loss. Therefore, determining how to accurately model the RR effectis also of crucial importance when in the strong-field regime. The Lorentz-Abraham-Dirac (LAD)equation describes the radiation reaction in the semi-classical perspective [28]. However, this equa-tion has unphysical runaway solutions that can be avoided by instead using the Landau-Lifshitz(LL) equation [29], which was shown to contain all physical solutions of the LAD equation [30].There are other models appropriate for numerical implementation (their comparison can be foundin [31]), and most of them give similar results when applied to the semi-classical interaction config-urations accessible with near-term laser technology. However, only the LAD and LL models wereshown to be consistent classical limits of the QED description of an electron interacting with astrong plane wave [32].The numerical integration of electron motion in strong fields requires a very fine temporalresolution, especially when the electron is not ultra-relativistic [31]. This is independent of thechoice of the radiation reaction model and is true even when the radiation reaction is turnedoff. This calls for solutions like sub-cycling [33] or the exact pusher proposed in this manuscript.The usual way to implement the additional RR force in PIC codes is (1) to integrate the particletrajectory using a pusher (splitting or exact) solely for the Lorentz force and then (2) to add animpulse from the RR force separately [34, 25]. This splitting process is simple to implement but canlead to the accumulation of errors in simulations with a large number of time steps, even thoughthe RR effect is perturbative. The particle pusher proposed in this article extends the analyticalsolution to include the LL equation. The only assumption required is that the fields are constantwithin each simulation time step. Therefore, the proposed pusher is free of numerical errors causedby splitting the operator for the Lorentz force. Although the motivation for developing an analyticpusher was to handle ultra-high fields, such a pusher will also accurately model the motion of arelativistic particle when E + v × B ≈
0. The analytic pusher (or any sub-cycling approach [33])will exhibit some errors in particle trajectories from assuming constant fields during a time step.Therefore, the time step must properly resolve the evolution of the fields as well.The PIC method is also beginning to be used to study the production of spin-polarized particlebeams. Furthermore, there is also a growing interest [35, 36, 37] in how RR affects particle-spindynamics in strong fields. Particle spin precession follows the Bargmann-Michel-Telegdi (BMT)equation [28], in which the phase space ( x , u ) is used to evaluate the spin s , where x and u areposition and momentum, respectively. Therefore, the effect of radiation reaction on the phasespace trajectories ( x , u ) will be also manifested in the behavior of spin dynamics. However, thereis far less literature directly related to the numerical schemes of the spin “push” than those ofthe momentum “push”. A typical numerical method [38] is similar to the Boris scheme: the spinorbit motion is approximated to be a pure rotation with a frequency that is evaluated with thetime-centered values of the electromagnetic fields and particle momentum. This Boris-like schemeis subject to large numerical errors in the strong-field regime as will be shown later. In this article,we also derive semi-analytic solutions to the BMT equation by utilizing the analytic expressionsof particle momentum without radiation reaction to advance the spin within an interval of time.The RR is then included as an impulse. During the next interval of time, the initial conditions3or the analytical update of the momentum are thus different, impacting the spin evolution duringsubsequent time intervals. Obtaining a fully analytic solution to the BMT equation in the presenceof radiation reaction is extremely difficult and likely impossible; however, the RR force can stillbe accurately included via the aforementioned splitting correction method. We note that thissemi-analytic approach can also be applied when quantum effects for RR are included. We leavecomparisons of examples with QED for a later publication.Although we are focusing on finding analytical solutions for both the momentum and positionadvance during intervals of time where the fields are constant, it is still important to relate this tothe leapfrog time indices in a standard PIC code. In most PIC codes, the position and momentum(proper velocity) are staggered in time such that x are known at half-integer values of time and u are known at integer values of time. For a given time step n , the fields are assumed constant duringthe particle push for the interval of time between n ∆ t and ( n + 1)∆ t ; the field values are assumedto be given at time ( n + 1 / t , requiring that particle positions are also assumed to be knownat time ( n + 1 / t . This implicitly assumes that the particle’s position does not change duringa time step. Under these conditions, we look to analytically advance the momentum forward fromtime n ∆ t to ( n + 1)∆ t . Although we may wish to then advance the particle position analyticallyto time ( n + 3 / t (assuming d x /dt = u /γ ), this can only be done during time intervals for which u is known—only until ( n + 1)∆ t for this example. Therefore, the proposed analytical pusher for astandard PIC code is really only doing an analytical advance of the particle momentum. However,the pusher may still lead to significant improvements in accuracy since the momentum advancecan lead to much larger errors than the position advance. This is easy to see by noting that theparticle’s speed is limited by the speed of light, from which it follows that during a time step aparticle can only move a fraction of a cell for any field strength. On the other hand, for ultra-strongfields the change in the proper velocity during a time step can be many orders of magnitude, i.e.,∆ u/u (cid:29)
1, during a time step. Nevertheless, the leapfrog advance in the position still leads tonoticeable errors compared to an analytic advance in position, as will be shown in a later section.The use of the analytic solutions together with the PSATD—or new concepts where the positionand momentum are defined at the same time—may lead to new PIC time-indexing algorithms.The remainder of the paper is organized as follows: In Section 2 and in Appendix A, wederive the equations for an analytic push of the Lorentz force and introduce the mathematicalformalism that can be extended to include the LL and BMT equations. In Section 3, we use themathematical formalism from Appendix A to obtain an analytic particle pusher for the 6D phasespace including the LL equation. These solutions are exact if the fields are constant during aninterval of proper time. In both sections, we also show how to obtain a mapping between the timestep in the lab frame and the proper time step. In Section 4, we derive the analytic solutionsto the BMT equation by employing the analytical solutions of momentum obtained in Section 3.The workflow and implementation of the proposed pusher for the 9D phase space are describedin Section 5. In Section 6, we first show simulation results using the proposed pusher for a singleparticle in an ultra-intense laser field propagating in vacuum, along with a comparison of resultsusing the standard Boris and Higuera-Cary pushers along with a Boris-like scheme for the spinpush. It is shown that the conventional numerical methods lead to large errors in the advance of9D phase space, while the proposed method provides accurate results. We then conduct full PICsimulations using
Osiris [39, 40] to investigate the difference in collective particle behavior usingthe proposed and conventional pushers. The performance of the proposed and regular pushers iscompared in Section 7. A summary and directions for future work are given in Section 8.4 . Particle motion in constant and uniform fields without radiation reaction
In this section, we will present a derivation of exact solutions to both the momentum andposition updates for constant fields. Analytic expressions can be obtained in various ways. P´etri [27]introduced a Lorentz-boosted frame where the E and B fields are parallel, for which analyticsolutions are possible. The analytic solutions then need to be transformed back to the lab frame.He also provided a mapping between the boosted (proper) and lab frame time steps. Gordonet al. [25] showed that the momentum update can be solved analytically in a covariant formand described a matrix representation of the analytic solution. However, Gordon et al. neitherprovided a mapping between the proper and lab frame time steps nor addressed special cases thatneed to be considered. As noted above, Gordon and Hafizi, [26] very recently proposed a specialunitary pusher which provides a second order accurate mapping in the absence of RR. Althoughthe underlying mathematics for obtaining solutions is different, each of the above approaches yieldsthe same net result for the cases considered. However, the forms for the solutions can have differentdegrees of algorithmic complexity. Here, we will present another method for finding an analyticexpression that is more compact and easier to implement into a PIC code. We use the covariantform for the equations of motion. More importantly, the mathematical formalism we use will beextended to include the LL and the BMT equations in later sections.The covariant form of the equation of motion without radiation reaction isd u µ d τ = qmc F µν u ν , (2.1)where u µ is the four-velocity, τ is the proper time, q is the particle charge and m is the particlemass. The field tensor is written as F µν = E E E E B − B E − B B E B − B . (2.2)To avoid rewriting constant factors, we use normalized physical quantities, i.e., τ → ω τ , q → qe , m → mm e and F µν → eF µν m e ω c , where e is the elementary charge, m e is the rest mass of electron and ω is a characteristic reference frequency which, for instance, can be chosen to be the electronplasma frequency or the laser frequency. In addition to the above normalization, we also absorbthe charge-to-mass ratio into F µν , i.e., F µν → qm F µν , to further simplify the expressions. Unlessotherwise specified, for the remainder of the paper we will use F to denote the tensor F µν . Thenormalized equation of motion is then given byd u d τ = F u. (2.3)If the elements of F are all constant in τ , then it is clear that this equation is easily solved if weknow the eigenvalues ( λ ) and eigenvectors of F . In Appendix A, it is shown that the field tensorhas four eigenvalues that come in pairs. One pair is real, given by λ = ± κ , and the other pair ispurely imaginary, given by λ = ± i ω , where κ = 1 √ (cid:114) I + (cid:113) I + 4 I , ω = 1 √ (cid:114) −I + (cid:113) I + 4 I , (2.4)5nd I = | E | − | B | , I = E · B (2.5)are Lorentz invariants.In order to obtain general solutions, it is important to project the initial values of the position, x ν , and proper velocity, u ν , four-vectors onto the eigenvectors. To facilitate this, the vector spaceof F can be split into two subspaces that are each expanded by two eigenvectors, i.e., S κ =span { e κ , e − κ } and S ω = span { e i ω , e − i ω } , where e λ is the eigenvector associated with the eigenvalue λ . It can be shown that (see Appendix A) S κ and S ω are mutually orthogonal in the sense of thefour-vector inner product. In this article, the four-vector inner product denoted by ( ·|· ) is definedas the contraction of two four-vectors, i.e., ( U | V ) = U µ V µ or ( U | V ) = U T GV in the matrix form,where G ≡ diag { , − , − , − } is the metric tensor. The modulus or length of a four-vector V isthus defined as | V | ≡ (cid:112) ( V | V ). In this article, we will decompose some physical quantities into S κ and S ω to simplify the mathematical derivation. The decomposition or projection can be achievedby applying the projection operator to a physical four-vector of interest, U , i.e., U κ = P κ U and U ω = P ω U where P κ and P ω are the projection operators defined as (see Appendix A), P κ = ω I + F κ + ω , P ω = κ I − F κ + ω , (2.6)where I is the 4 × u κ and u ω separately. Taking the propertime derivative of both sides of Eq. (2.3) and using the properties F u κ = κ u κ and F u ω = − ω u ω (see Appendix A), the equation of motion can be decomposed intod u κ d τ = κ u κ (2.7)and d u ω d τ = − ω u ω . (2.8)The solutions are given by u κ ( τ ) = u κ cosh( κτ ) + κ − F u κ sinh( κτ ) (2.9)and u ω ( τ ) = u ω cos( ωτ ) + ω − F u ω sin( ωτ ) , (2.10)where u κ = u κ ( τ = 0) and u ω = u ω ( τ = 0), with each obtained via u κ = P κ u and u ω = P ω u .The four-position can then be obtained by directly integrating the expressions for the propervelocity over the proper time to give x κ ( τ ) − x κ = κ − (cid:2) u κ sinh( κτ ) + κ − F u κ (cosh( κτ ) − (cid:3) , (2.11) x ω ( τ ) − x ω = ω − (cid:2) u ω sin( ωτ ) − ω − F u ω (cos( ωτ ) − (cid:3) , (2.12)where the initial position components x κ and x ω are likewise obtained by x κ = P κ x and x ω = P ω x .These equations represent analytic expressions for how to advance the particle four-velocityand four-position from initial to final values during an interval of the proper time in absence of6adiation reaction. The 6D phase space evolution could therefore be advanced during an interval ofthe proper time, i.e., a proper time step ∆ τ . However, in a PIC simulation, the fields are advancedusing the lab-frame time step ∆ t . We therefore need to advance forward the phase space for fixedlab-frame steps rather than a fixed proper time step for each particle. Although ∆ τ changes foreach simulation (observer) time step, a mapping between the lab and proper time intervals forfixed fields can be found using the time-like component of Eqs. (2.11) and (2.12):∆ t = x (∆ τ ) − x (∆ τ ) = x κ (∆ τ ) + x ω (∆ τ ) − x κ − x ω , (2.13)or the explicit form,∆ t = κ − (cid:20) u κ sinh( κτ ) + u · E κ (cosh( κτ ) − (cid:21) + ω − (cid:20) u ω sin( ωτ ) − u · E ω (cos( ωτ ) − (cid:21) . (2.14)This is a transcendental equation consisting of trigonometric and hyperbolic functions. We usuallyneed to resort to some root-finding algorithms such as the Newton-Raphson method with second-order precision or the Householder method with higher precision to seek the solution. We emphasizethat Eq. (2.14) is only valid for κ (cid:54) = 0 and ω (cid:54) = 0, while Eq. (2.13) is generally correct. Therefore,in the situation where only κ or only ω vanishes, the x κ − x κ or x ω − x ω terms in Eq. (2.13) shouldbe replaced with the corresponding expressions that will be presented in Sections 2.1 and 2.2. Sincethe subspace decomposition fails when κ and ω simultaneously vanish, the complete solution fordisplacement, i.e., x − x (given in Section 2.3), should be used for the time step mapping.The above analytic solutions are mathematically well behaved for non-vanishing eigenvalues.However, computational precision might be lost somewhat due to round-off errors when κ and/or ω become numerically small. Therefore, we use Taylor expansions of the general solutions whenever κ or ω become smaller than some threshold value, (cid:15) th . The threshold should be properly selectedso that the Taylor expansion is sufficiently accurate, e.g., reaching the machine precision.One, but not both, of the eigenvalues κ or ω vanishes when E and B are orthogonal to oneanother ( I = 0) but not equal in amplitude ( I (cid:54) = 0). When I >
0, i.e., | E | > | B | , we have ω → κ (cid:54) = 0, while for the opposite case when I <
0, i.e., | E | < | B | , we have κ → ω (cid:54) = 0.When both I → I →
0, i.e., E and B are equal in amplitude and orthogonal, then both κ → ω →
0. Next we describe how we handle these special cases separately. I → and I > , or equivalently ω → and κ (cid:54) = 0In this case, it should be noted that the subspace decomposition is still valid because theprojection operators in Eq. (2.6) are well-defined unless κ and ω vanish simultaneously. If E and B are orthogonal with | E | > | B | , we have ω = 0. The exact solution of u ω can be obtained byexpanding Eq. (2.10) and then taking the limit ω →
0. The Taylor expansion of Eq. (2.10) is u ω ( τ ) = u ω + F u ω τ − (cid:18) u ω + 16 F u ω τ (cid:19) ω τ + O [( ωτ ) ] , (2.15)where we note that F u ω is O ( ω ). The exact solution of the four-position can be obtained eitherby taking the limit of the Taylor expansion of Eq. (2.12) or by integrating the Taylor expansionfor u ω , both giving x ω ( τ ) − x ω = u ω τ + 12 F u ω τ − (cid:18) u ω + 124 F u ω τ (cid:19) ω τ + O ( ω τ ) . (2.16)The total value for u and x can be obtained by summing with full expressions for u κ and x κ ,respectively. 7 .2. Special case: I → and I < , or equivalently κ → and ω (cid:54) = 0If on the other hand E and B are orthogonal with | B | > | E | , we have κ →
0, and the expressionsfor u κ and x κ are analogously given by u κ ( τ ) = u κ + τ F u κ + (cid:18) u κ + 16 F u κ τ (cid:19) κ τ + O [( κτ ) ] (2.17)and x κ ( τ ) − x κ = u κ τ + 12 F u κ τ + (cid:18) u κ + 124 F u κ τ (cid:19) κ τ + O ( κ τ ) . (2.18) I → and I → , or equivalently ω → and κ → E and B are mutually orthogonal and equal in magnitude, I → I → κ and ω approach zero and the subspace decomposition is thus no longer valid.Nevertheless, the full solutions are still valid, and thus we can obtain the Taylor expansions for u by adding Eqs. (2.9) and (2.10) together and then taking the limit κ, ω →
0. This results in u ( τ ) = u + F u τ + 12 F u τ + 16 F u τ + O [( κτ ) , ( ωτ ) ] . (2.19)The four-position can be similarly obtained by taking the Taylor expansion of the sum of Eqs. (2.11)and (2.12), x ( τ ) − x = u τ + 12 F u τ + 16 F u τ + 124 F u τ + O ( κ τ , ω τ ) . (2.20)It should be pointed out that F → κ, ω → u for F → κτ ) n ( ωτ ) − n , ( n = 0 , , ,
3. Particle motion in constant and uniform fields with radiation reaction
In this section, we will derive the exact solutions to the LL equation by utilizing the orthogonal-ity of S κ and S ω introduced in the previous section. We will see that by splitting the four-velocity u into components belonging to the two subspaces, i.e., u κ and u ω , the integration of the LL equationis greatly simplified. The covariant form for the LL equation can be written asd u µ d τ = qmc F µν u ν + 2 q m c (cid:18) ∂F µν ∂x i u ν u i − qmc F µν F νi u i + qmc ( F ij u j )( F ki u k ) u µ (cid:19) , (3.1)where x i is the four-position. We are investigating cases where the fields are assumed constant inthe proper time during a time step. Furthermore, it has been shown by others [34] that the firstterm in the parentheses with the partial derivatives of x i can be neglected. This is referred to asthe reduced Landau-Lifshitz model. After normalizing all quantities as described in Section 2, thereduced LL equation can be written asd u d τ = F u + σ q m (cid:2) F u − ( u | F u ) u (cid:3) , (3.2)where σ is a dimensionless parameter defined as σ = e ω m e c .8y utilizing the subspace decomposition for u , i.e., u = u κ + u ω , and recalling the relations F u κ = κ u κ and F u ω = − ω u ω , it can be shown that the contraction ( u | F u ) becomes ( u | F u ) = κ ( u | u κ ) − ω ( u | u ω ) = κ | u κ | − ω | u ω | . Substituting this result into the reduced LL equation (3.2)and using the fact that the four-velocity has unit length, i.e., | u | = | u κ | + | u ω | = 1, we obtaintwo decoupled nonlinear differential equations,d u κ d τ = F u κ + α (1 − | u κ | ) u κ , (3.3)d u ω d τ = F u ω − α (1 − | u ω | ) u ω , (3.4)where α ≡ σ q m ( κ + ω ). To solve the nonlinear ordinary differential equation (3.3) [Eq. (3.4)can be solved in an analogous manner], we first construct a trial solution as the product of theamplitude of u κ and a four-vector w κ , i.e., u κ = | u κ ( τ ) | w κ . This implies that w κ is also enforcedto have unit length. With this assumption Eq. (3.3) can be separated into two ODEs asd w κ d τ = F w κ , (3.5)d | u κ | d τ = α (1 − | u κ | ) | u κ | . (3.6)The first ODE is exactly the unperturbed Lorentz equation. It implies that the modulus of w κ does not change, which can be justified by left multiplying w κ on both sides of the equation andusing the property ( w κ | F w κ ) = 0 as described in Appendix A. As discussed in Section 2, theunperturbed Lorentz equation (3.5) has the solution w κ ( τ ) = w κ cosh( κτ ) + κ − F w κ sinh( κτ ) , (3.7)where w κ = w κ ( τ = 0).Eq. (3.6) can be directly integrated to obtain a solution to the amplitude equation, | u κ ( τ ) | = | u κ | (cid:112) | u κ | + | u ω | e − α τ . (3.8)Combining the solutions for w κ and | u κ | yields u κ ( τ ) = 1 (cid:112) | u κ | + | u ω | e − α τ (cid:2) u κ cosh( κτ ) + κ − F u κ sinh( κτ ) (cid:3) . (3.9)The solution to u ω can be obtained in an analogous way and is given by u ω ( τ ) = 1 (cid:112) | u ω | + | u κ | e α τ (cid:2) u ω cos( ωτ ) + ω − F u ω sin( ωτ ) (cid:3) . (3.10)Although there is no simple and closed-form expression for the four-position if radiation reactionis included, it is still possible to obtain approximate expressions with sufficiently high accuracyas long as the “friction” coefficient α is much less than ∆ τ . Simple estimates can show thatthis premise is often true for problems of interest. According to its definition, we know that α = ( π r e λ ) q m (cid:112) I + 4 I ≤ ( π r e λ ) q m ( | E | + | B | ), where r e is the classical electron radius. The9quality is true if and only if E and B are parallel, i.e., E · B = | E || B | . For example, assumingthe characteristic length λ ∼ µ m and the normalized field strengths E and B are on the orderof 10 , we get α ∼ − . In simulations, the time step must be properly selected to sufficientlyresolve the characteristic time scales, say ∆ t ∼ .
1, and thus ∆ τ ∼ ∆ t/γ ≤
1. Therefore the upperlimit of α τ is on the order of 10 − when 0 < τ < ∆ τ , and keeping only the first term in the Taylorexpansions of the denominator in Eqs. (3.9) and (3.10) is consequently valid. The four-position x κ can be approximately given by integrating the lowest-order expansion of Eq. (3.9), x κ ( τ ) − x κ = κ − (cid:2) u κ sinh( κτ ) + κ − F u κ cosh( κτ ) (cid:3) (1 + α | u ω | τ ) − κ − F u κ − α | u ω | κ − (cid:2) u κ (cosh( κτ ) −
1) + κ − F u κ sinh( κτ ) (cid:3) + O [( α τ ) ] . (3.11)Similarly, we have x ω ( τ ) − x ω = ω − (cid:2) u ω sin( ωτ ) − ω − F u ω cos( ωτ ) (cid:3) (1 − α | u κ | τ ) + ω − F u ω − α | u κ | ω − (cid:2) u ω (cos( ωτ ) −
1) + ω − F u ω sin( ωτ ) (cid:3) + O ( α τ ) . (3.12)These expressions can then be used to approximately obtain the time step mapping using Eq. (2.13),and the fast root-finding algorithms mentioned previously in Section 2 are still applicable.Next we will discuss the special cases where κ and/or ω vanish as in Sections 2.1–2.3. I → and I > , or equivalently ω → and κ (cid:54) = 0When E and B are orthogonal and | E | > | B | , u ω can be obtained by simply expanding Eq. (3.10)and then taking the limit ω →
0. However, instead of conducting a full expansion, we only expandthe terms inside the bracket of Eq. (3.10) because the leading factor outside the bracket does notcause singularity when ω →
0. This also helps to keep the form simple. The expansion can bewritten as u ω ( τ ) = 1 (cid:112) | u ω | + | u κ | e α τ (cid:20) u ω + F u ω τ − (cid:18) u ω + 16 F u ω τ (cid:19) ω τ (cid:21) + O [( ωτ ) ] . (3.13)The solution of u κ still follows Eq. (3.9). Conducting a Taylor expansion in ω for Eq. (3.12), thefour-position can be written as x ω ( τ ) − x ω = u ω τ + 12 ( F u ω − α | u κ | u ω ) τ − α | u κ | F u ω τ − (cid:18) u ω −
124 (3 α | u κ | u ω − F u ω ) τ − α | u κ | F u ω τ (cid:19) ω τ + O ( α τ , ω τ ] . (3.14)The solution of x κ still follows Eq. (3.11). I → and I < , or equivalently κ → and ω (cid:54) = 0Similarly, if E and B are orthogonal and | B | > | E | , we have κ = 0 and the Taylor expansionsof u κ and x κ are given by u κ ( τ ) = 1 (cid:112) | u κ | + | u ω | e − α τ (cid:20) u κ + F u κ τ + (cid:18) u κ + 16 F u κ τ (cid:19) κ τ (cid:21) + O [( κτ ) ] (3.15)10nd x κ ( τ ) − x κ = u κ τ + 12 ( F u κ + α | u ω | u κ ) τ + 13 α | u ω | F u κ τ + (cid:18) u κ + 124 (3 α | u ω | u κ + F u κ ) τ + 130 α | u ω | F u κ τ (cid:19) κ τ + O ( α τ , κ τ ) . (3.16) I → and I → , or equivalently ω → and κ → κ → ω → E and B are orthogonal and equal in magnitude). The analytic solution can be sought byfirst replacing u κ and u ω with P κ u and P ω u , respectively, and then taking the limit ω, κ → | u κ | and | u ω | also needs to be expressed in terms of u . It can be shown(see Appendix B) that | u κ | = ω − | F u | κ + ω , | u ω | = κ + | F u | κ + ω . (3.17)Adding Eqs. (3.9) and (3.10) together yields u ( τ ) = e α τ (cid:2) u κ cosh( κτ ) + κ − F u κ sinh( κτ ) (cid:3) + e − α τ (cid:2) u ω cos( ωτ ) + ω − F u ω sin( ωτ ) (cid:3)(cid:112) | u κ | e α τ + | u ω | e − α τ . (3.18)In order to keep a simple form, we take the Taylor expansion of the numerator and denomina-tor separately, rather than seeking a full expansion. Inserting Eqs. (3.17) and (A.5), the Taylorexpansion is given by u ( τ ) = u + ( F u + ˜ σ F u ) τ + F u τ (cid:112) − σ y τ + ( ω − κ )˜ σ τ + O ( κ τ , ω τ , κ ω τ )+ ( ω − κ )( u + F u τ )˜ σ τ + F u (cid:0) τ + ˜ σ τ (cid:1) + O ( κ τ , ω τ , κ ω τ ) (cid:112) − σ y τ + ( ω − κ )˜ σ τ + O ( κ τ , ω τ , κ ω τ ) , (3.19)where ˜ σ ≡ σ q m and y ≡ | F u | .Similarly, we can obtain the expression for x ( τ ) by adding Eqs. (3.11) and (3.12) and inserting(3.17) and (A.5), giving x ( τ ) − x = u τ + 12 (cid:0) F u + ˜ σ y u + ˜ σ F u (cid:1) τ + 13 (cid:18) ˜ σ y F u + 12 F u (cid:19) τ + 18 ˜ σ y F u τ + 13 (cid:18) ˜ σ + 18 τ + 110 ˜ σ y τ (cid:19) F u τ + O (˜ σ τ , κ τ , ω τ ) . (3.20)
4. Spin precession in uniform and constant fields
In this section, we will derive the semi-analytic solutions to the particle four-spin vector inuniform and constant fields by utilizing the analytic expression of the four-velocity in absence ofRR. After obtaining analytic solutions for the spin evolution based on the analytic evolution of u s µ d τ = qmc (cid:20) g F µν s ν − c (cid:16) g − (cid:17) ( u i F ij s j ) u µ (cid:21) , (4.1)where g is the Land´e g-factor and is dimensionless.The four-spin s µ here is described in the observer frame, and hence its time-like component isnonzero. However, as an intrinsic property, it is more conventional to investigate the spin precessiondynamics in the particle rest frame. Therefore, we need to transform s µ to the particle rest frameafter solving the BMT equation. Using normalized units and absorbing the qm factor into F asdone in the previous two sections, the BMT equation can be written asd s d τ = (1 + a ) F s − a ( u | F s ) u, (4.2)where a ≡ g − a (cid:39) . f ≡ ( u | F s ) and show that it can be analyticallysolved even without knowing how s evolves. We define f ( τ ) = ( u κ | F s κ ) + ( u ω | F s ω ) ≡ f κ ( τ ) + f ω ( τ )and then split Eq. (4.2) into S κ and S ω based on the eigenvalues of F as was done for the propervelocity: d s κ d τ = (1 + a ) F s κ − af ( τ ) u κ , (4.3)d s ω d τ = (1 + a ) F s ω − af ( τ ) u ω . (4.4)Combining Eqs. (3.3) and (4.3) and using the fact that F s κ = κ s κ and ( u κ | F u κ ) = 0, it followsthat the time derivative of f κ is d f κ d τ = aκ ( u κ | s κ ) . (4.5)Similarly, the time derivative of f ω is d f ω d τ = − aω ( u ω | s ω ) . (4.6)The quantity I ≡ ω f κ − κ f ω is an invariant. This can be readily verified by taking the appropriatelinear combination Eqs. (4.5) and (4.6),d I d τ = aω κ [( u κ | s κ ) + ( u ω | s ω )] = aω κ ( u | s ) = 0 . (4.7)Here we have also used the fact ( u | s ) = 0, which follows from the fact that the time-like componentof four-spin in the particle rest frame is zero, i.e., according to the Lorentz transformation s (cid:48) =12 s − u · s ≡ ( u | s ) = 0. Taking the time derivative of Eq. (4.5) and substituting in Eqs. (3.3) and(4.3) gives d f κ d τ = a κ (cid:0) | u ω | f κ − | u κ | f ω (cid:1) . (4.8)We can similarly get the second-order ODE for f ω ,d f ω d τ = − a ω (cid:0) | u κ | f ω − | u ω | f κ (cid:1) . (4.9)Adding these two ODEs together and using the relations | u κ | + | u ω | = 1 and f = f κ + f ω , wefinally arrive at d f d τ = − a Ω f + a I , (4.10)where Ω = ω | u κ | − κ | u ω | = ω | u κ | − κ | u ω | (note that | u κ | and | u ω | are constant withoutRR). It should be noted that Ω is always positive due to | u κ | ≥ | u ω | ≤ f ( τ ) = − Ω − ( I − f Ω ) cos( a Ω τ ) + ( a Ω) − ˙ f sin( a Ω τ ) + Ω − I , (4.11)where we have used the initial conditions f = ( u κ | F s κ ) + ( u ω | F s ω ) and ˙ f = a [ κ ( u κ | s κ ) − ω ( u ω | s ω )]. After obtaining the solution to f ( τ ), we insert it back into Eqs. (4.3) and (4.4) tosolve for s κ and s ω . Eqs. (4.3) and (4.4) can now be treated as inhomogeneous ODEs, and thecomplete solutions are the sum of the homogeneous and inhomogeneous solutions, i.e., s κ = ¯ s κ + ˜ s κ and s ω = ¯ s ω + ˜ s ω . The homogeneous solutions satisfy¨¯ s κ = (1 + a ) κ ¯ s κ , ¨¯ s ω = − (1 + a ) ω ¯ s ω . (4.12)We impose the initial conditions ¯ s λ (0) = s λ and ˙¯ s λ (0) = (1 + a ) F s λ ( λ = κ, ω ) to the aboveODEs, which implies that the inhomogeneous solutions must satisfy the initial conditions ˜ s λ (0) = 0and ˙˜ s λ (0) = − af u λ ( λ = κ, ω ). Solving the above homogeneous ODEs yields¯ s κ ( τ ) = s κ cosh[(1 + a ) κτ ] + κ − F s κ sinh[(1 + a ) κτ ] , (4.13)¯ s ω ( τ ) = s ω cos[(1 + a ) ωτ ] + ω − F s ω sin[(1 + a ) ωτ ] . (4.14)Since the inhomogeneous terms in Eqs. (4.3) and (4.4) include u κ and u ω , we can thus constructthe trial solution of ˜ s κ as the linear combination of u κ and ˙ u κ , and that of ˜ s ω as the linearcombination of u ω and ˙ u ω , i.e.,˜ s κ = C κ ( τ ) u κ + D κ ( τ ) ˙ u κ , ˜ s ω = C ω ( τ ) u ω + D ω ( τ ) ˙ u ω . (4.15)According to the initial conditions to which ˜ s κ , ˜ s ω and their time derivatives must be consist with,we have that the coefficients in Eq. (4.15) must meet the initial conditions C λ (0) = 0, ˙ C λ (0) = − af , D λ (0) = 0 and ˙ D λ (0) = 0 ( λ = κ, ω ). A set of first-order ODEs for these coefficients can be foundby inserting Eq. (4.15) into (4.3) and (4.4) and comparing the coefficients of u κ , ˙ u κ , u ω and˙ u ω . Therein we have used ¨ u κ = κ u κ and ¨ u ω = − ω u ω . The detailed process for solving these13oefficients is tedious and can be found in Appendix C. Here, we directly list the final results. Forthe general case, we have C κ = ˙ f cos( a Ω τ ) − cosh( aκτ ) a ( κ + Ω ) + h Ω κ sin( a Ω τ ) − h κ Ω sinh( aκτ ) κ Ω( κ + Ω ) ,D κ = ˙ f k sin( a Ω τ ) − Ω sinh( aκτ ) aκ Ω( κ + Ω ) − h Ω κ cos( a Ω τ ) + h κ Ω cosh( aκτ ) κ Ω ( κ + Ω ) + I κ Ω , (4.16)and C ω = ˙ f cos( aωτ ) − cos( a Ω τ ) a ( ω − Ω ) + h ω Ω sin( aωτ ) − h Ω ω sin( a Ω τ ) ω Ω( ω − Ω ) ,D ω = ˙ f Ω sin( aωτ ) − ω sin( a Ω τ ) aω Ω( ω − Ω ) − h ω Ω cos( aωτ ) − h Ω ω cos( a Ω τ ) ω Ω ( ω − Ω ) − I ω Ω , (4.17)where we have defined h κ ≡ I + f κ , h ω ≡ I − f ω and h Ω ≡ I − f Ω to simplify the notation.The singularities appearing in these coefficients and in Eqs. (4.13) and (4.14) must be treatedwith care. Apart from the singularities caused by either κ → ω →
0, the characteristicfrequency Ω in the denominators of Eqs. (4.16) and (4.17) will also bring about singularities de-pending on the values of κ and ω . In the following discussion about these special cases, we willseek the exact solutions or the leading terms of the Taylor expansions for C κ , D κ , C ω and D ω bytaking the limits of Eqs. (4.16) and (4.17). ω (cid:54) = 0 and κ (cid:54) = 0In the general case where ω (cid:54) = 0 and κ (cid:54) = 0, the singularity can only be caused by the factor ω − Ω in the expressions of C ω and D ω when Ω → ω . The u κ and u ω included in the inhomogeneoussolutions still follow Eqs. (2.9) and (2.10), and ¯ s κ and ¯ s ω still follow Eqs. (4.13) and (4.14). Theexact solution of the coefficients C ω and D ω can be found by taking the limit Ω → ω of Eq. (4.17): C ω = ah ω τ ω cos( aωτ ) + 12 (cid:32) h ω − I ω − ˙ f τω (cid:33) sin( aωτ ) ,D ω = − I ω + (cid:32) I ω + ˙ f τ ω (cid:33) cos( aωτ ) − (cid:32) ˙ f aω − ah ω τ ω (cid:33) sin( aωτ ) . (4.18) ω (cid:54) = 0 and κ → u κ . The homogeneous solution ¯ s κ can befound by taking the limit κ → s κ = s κ + (1 + a ) F s κ τ + (cid:18) s κ + 16 (1 + a ) F s κ τ (cid:19) (1 + a ) κ τ + O ( κ τ ) . (4.19)Since the factor ω − Ω in the expressions of C ω and D ω will always cause singularity whenΩ → ω , independent of the values of κ and ω , the discussion of these coefficients here (and the onethat follows) can be divided into two parts according to whether or not Ω → ω .14f Ω (cid:57) ω , there is no singularity in C ω and D ω , so Eq. (4.17) is still valid. Taking the limit κ → C κ and D κ are given by C κ = − aτ I Ω + ˙ f a Ω (cos( a Ω τ ) −
1) + h Ω Ω sin( a Ω τ ) ,D κ = − ˙ f τ Ω − a τ I − h Ω Ω (cos( a Ω τ ) −
1) + ˙ f a Ω sin( a Ω τ ) . (4.20)When Ω → ω , C ω and D ω follow Eq. (4.18), and Eq. (4.20) can be used to calculate C κ and D κ simply by replacing Ω with ω . ω → and κ (cid:54) = 0In this case, u ω in the inhomogeneous solution ˜ s ω should be evaluated using Eq. (2.15). Thehomogeneous solution ¯ s ω can be obtained by taking the limit ω → s ω = s ω + (1 + a ) F s ω τ − (cid:18) s ω + 16 (1 + a ) F s ω τ (cid:19) (1 + a ) ω τ + O ( ω τ ) . (4.21)If Ω (cid:57) ω →
0, only C ω and D ω present singularity, and the solutions are given by C ω = − aτ I Ω + ˙ f a Ω (cos( a Ω τ ) −
1) + h Ω Ω sin( a Ω τ ) ,D ω = − ˙ f τ Ω − a τ I − h Ω Ω (cos( a Ω τ ) −
1) + ˙ f a Ω sin( a Ω τ ) , . (4.22)If Ω → ω →
0, all the coefficients present singularities. The solutions can be written as C κ = aτ I κ − ˙ f aκ (cosh( aκτ ) − − (cid:18) I κ + f κ (cid:19) sinh( aκτ ) ,D κ = ˙ f τκ + a τ I κ − (cid:18) I κ + f κ (cid:19) (cosh( aκτ ) − − ˙ f aκ sinh( aκτ ) (4.23)and C ω = − aτ (cid:18) f + 12 ˙ f τ + 16 a τ I (cid:19) ,D ω = − a τ (cid:18) f + 16 ˙ f τ + 124 a τ I (cid:19) . (4.24) ω → and κ → ω → κ →
0, it can be shown that C κ = C ω ≡ C and D κ = D ω ≡ D , and thecoefficients are given by C = − aτ (cid:18) f + 12 ˙ f τ (cid:19) , D = − a τ (cid:18) f + 16 ˙ f τ (cid:19) . (4.25)As previously stated, the subspace decomposition fails in this situation. The homogeneous solution¯ s can be found by adding Eqs. (4.19) and (4.21) together, applying the relations s κ = P κ s and s ω = P ω s and then taking the limit κ, ω →
0, i.e.,¯ s = s + (1 + a ) τ F s + 12 (1 + a ) τ F s . (4.26)15ote that since F = 0 when both κ and ω vanish, Eq. (4.26) actually gives the exact solution. Theinhomogeneous solution ˜ s can be also obtained by simply adding ˜ s κ and ˜ s ω and then evaluating at κ, ω →
0, i.e., ˜ s = C ( τ ) u ( τ ) + D ( τ ) ˙ u ( τ ), where u ( τ ) and ˙ u ( τ ) can be found from Eq. (2.19).
5. Algorithm workflow
In this section, we will introduce the algorithm workflow using the analytical expressions of 9Dphase space obtained in previous sections. Depending on the purposes and program implementa-tions of PIC codes, we provide four distinct algorithms that are characterized by different choicesfrom the subset of the analytical solutions. Figure 1 shows the numerical workflow of the fouralgorithm implementations. The numbers of the requisite equations for each algorithm have beensummarized in each block.For existing PIC codes, the momentum and position of the particles are staggered in time andthe fields are known at the same time as the position. As a result, only the red and blue pathsin Fig. 1 are possible without significantly reworking the PIC algorithm. Maintaining the leapfrogadvance of the position thus permits modification of only the momentum update, and the fieldsolve and current deposit do not have to be modified. If updating the spin is unimportant, thenthe red path is desirable, which uses the analytic solution for the momentum with RR included andthe leapfrog advance for the position. The blue path should be used when including spin dynamics,and the details for this method are similar to those of the yellow path described below.For PIC codes that define position and momentum at the same points in time, the yellow pathshould be used when the evolution of the full 9D phase space is important; otherwise the greenpath should be selected. The yellow path utilizes the analytical solutions to ( x , u , s ) without RR,but the effects of RR are incorporated by splitting the change in momentum due to RR into twohalf-impulses that are applied before and after the analytic solution without RR is used. In thetime interval n ∆ t < t < ( n + 1)∆ t , the first half-impulse can be applied to u via u − = u n + ∆ t f RR ( u n ) , (5.1)after which u − is pushed to u + using the analytical solutions for a full time step. The quantities x and s are also analytically advanced a full time step (for the blue path only s is analyticallyadvanced), where u − is used for the u values in the pertinent equations. The other RR half-impulseis then applied via u n +1 = u + + ∆ t f RR ( u + ) . (5.2)Here the RR force is evaluated as follows: f RR ( u ) = σ q γm [ F u − ( u | F u ) u ] spatial , (5.3)where the subscript “spatial” refers to the space-like component of a four-vector.If the position and momentum are staggered in time, the analytical expressions of x can nolonger be used; in the time interval n ∆ t < t < ( n + 1)∆ t where the solution to u is known, x isknown only within n ∆ t < t < ( n + 1 / t . Therefore, the positions need to be advanced in theconventional leapfrog manner, i.e. x n + = x n + + u n +1 ∆ t/γ n +1 . (5.4)16 dvance (2.9) (2.10)(2.9) (2.15)(2.10) (2.17)(2.19) advance (3.9) (3.10)(3.9) (3.13)(3.10) (3.15)(3.19) advance (2.11) (2.12)(2.11) (2.16)(2.12) (2.18)(2.20) advance (3.11) (3.12)(3.11) (3.14)(3.12) (3.16)(3.20) half RR impulse (5.1) half RR impulse (5.2) advance (4.13)-(4.17)+Sec. 4.1(4.13)(4.15)(4.16)+Sec. 4.3(4.14)(4.15)(4.17)+Sec. 4.2Sec. 4.4 calculate eigenvaluessubspace projection (2.4)(2.6) leapfrog advance (5.4) mapping time step (2.13) P3P4 P5P6
Figure 1: Numerical workflow of four algorithm implementations. The relevant equation numbers are summarizedin each block, and the blue, red, yellow and green paths correspond to the pusher combinations P3–P6, respectively.The red and green paths analytically advance particle momentum with radiation reaction (RR) included, but withoutconsidering spin. The blue and yellow paths analytically advance particle momentum and spin without consideringRR, but incorporate RR by applying two half-impulses at the beginning and end of the advance. The blue and redpaths advance the position through the standard leapfrog scheme used in most PIC codes, whereas the yellow andgreen paths advance the position analytically (requires a new time-indexing PIC algorithm). s can be defined on the same time grid pointsas u so that the analytical solutions still applies.
6. Example simulations
In this section, we will compare different particle pushers through a series of particle-trackingsimulations where (1) a single particle interacts with an ultra-intense laser pulse in prescribed fieldsand (2) many particles collectively interact with self-consistent fields in an
Osiris
PIC simulation.As we have multiple options to advance the particle position, momentum and spin, the followingschemes (P1–P6) will be investigated to see how accurately they advance the ( x , u , s ) phase space:1. P1 – The Boris pusher is used to advance the particle momentum, and the position isadvanced in a leapfrog manner with second-order accuracy in ∆ t . The RR force is addedaccording to the splitting method addressed in Section 5. Vieira’s scheme [38] is used toadvance the spin.2. P2 – The setup is identical to P1 except the Higuera-Cary pusher [24] is used to advance theparticle momentum.3. P3 to P6 refer to the blue, red, yellow and green paths in Fig. 1, respectively. In this section, we compare the various pushers using a particle-tracking code in which thefields are prescribed. This permits using the analytic position update as well. We first consider aone-dimensional case in which a laser pulse propagates in vacuum. Test particles are initialized infront of the laser pulse. The plane-wave laser is linearly polarized in the ˆ2-direction and moves inthe ˆ1-direction. The normalized vector potential is given by A = a cos (cid:18) πφ ω τ FWHM (cid:19) cos φ ˆ e (6.1)when the phase φ ≡ ω t − k x is within [ − ω τ FWHM , ω τ FWHM ], and vanishes otherwise. Here, τ FWHM is defined as the full-width-at-half-maximum of the field envelope, ω is the laser fre-quency and a is the strength parameter which is connected with the peak intensity via a =0 . (cid:113) I [10 W/cm ] λ [ µ m]. In all of the following comparisons, a pulse duration of τ FWHM =50 ω − is chosen, and the field is expressed analytically according to Eq. (6.1). We assume thelaser wavelength to be 0.8 µ m and set the reference frequency to be the laser frequency ω so thatthe dimensionless radiative damping parameter σ ≈ . × − .In the first set of simulations, the test particle has an initial momentum of p = − m e c (the negative sign means it counter-propagates relative to the laser), and the initial spin is alongthe positive ˆ1-direction. We tracked the transverse momentum p , phase φ and transverse spin s during the particle-wave interaction for various values of ∆ t . For relatively weak laser intensitieswhere a is on the order of unity, it is found that all the aforementioned numerical schemes providenearly identical and correct phase space trajectories. This is not the case for higher intensities.Figure 2 shows the results for a = 300 ( I ∼ . × W/cm ) for two values of ∆ t . Theblack dashed line is obtained using a fourth-order Runge-Kutta integrator with sufficiently smalltime step that it can be viewed as the “correct” result. It can be seen that the schemes which18se the split operator, i.e., standard particle pushers (P1 and P2), lead to incorrect results forboth ∆ t = 0 . ω − and ∆ t = 0 . ω − . The phase shift of particles pushed by P1 and P2 areseverely miscalculated [see Figs. 2(b) and (e)], which leads to a large deviation in the phase spacetrajectories. According to our tests, P1 and P2 do not converge until reducing ∆ t to ∼ . ω − . ForP3 and P4, which advance the position in a leapfrog manner and the momentum with the analyticalpusher (P3 also analytically advances spin), the momentum and spin oscillations and phase shiftare qualitatively correct, but quantitatively inaccurate for ∆ t = 0 . ω − [see Figs. 2(a)–(c)]. Whenthe time step is reduced to ∆ t = 0 . ω − , both P3 and P4 converge to the “correct” results asshown in Figs. 2(d)–(f). Since P5 and P6 advance both position and momentum analytically (P5also advances spin analytically), they give good agreement with the “correct” results for the twotime steps, as expected.We next tested how well these numerical schemes work with zero initial momentum as shownin Fig. 3. According to Vranic et al. [31], this situation is more sensitive to numerical noise. Due tothe energy loss during the laser-particle interaction, the particle will stay in phase for much longer,increasing the duration of interaction. We tested two time steps, ∆ t = 0 . ω − and ∆ t = 0 . ω − ,and again P1 and P2 lead to a large deviation from the correct results. For ∆ t = 0 . ω − , P3 andP4 lead to the correct phase space trajectory results for the first few cycles, but clear deviationsappear at later times due to the accumulation of numerical errors over a long duration. Goodagreement can be reached when the time step is reduced to ∆ t = 0 . ω − . As before, P5 and P6lead to excellent agreement with the correct results even for a time step typically used to accuratelysolve for the fields in laser-plasma-interaction simulations (∆ t = 0 . ω − ).We also examined the effect of using the proposed schemes for the situation where the testparticles are initialized inside the laser field. A stationary ( p = 0) test particle was initializedinside laser fields with a = 100 at a location where the laser electric field (vector potential) reachesa maximum (zero), i.e., φ = 0. The evolution of the transverse momentum p is shown in Fig. 4.We gradually reduced the time step of each scheme to examine the maximum ∆ t for which thesimulation result converges to that of the high-precision Runge-Kutta method (black dashed linesin Fig. 4). As we can see from Figs. 4(a) and (b), the schemes using the regular pushers do notconverge at a conventionally selected ∆ t (∆ t = 0 . . ω − ) that resolves the laser frequency,but require an extremely small time step ∆ t = 0 . ω − to converge. The maximum time stepfor P3 and P4 to converge is around ∆ t = 0 . ω − , as shown in Figs. 4(c) and (d). Therefore,the benefit of solely using the analytic momentum advance is twenty-fold compared to P1 and P2.However, P5 and P6 converge at an even larger time step ∆ t = 4 ω − for this specific problem, asshown in Figs. 4(e) and (f), giving a hundred-fold improvement over P3 and P4. It should be notedthat the comparison here is to show that the proposed pushers can greatly reduce the requirementfor time steps, rather than to give a rule of thumb for choosing a time step. The choice of timestep and the benefits of using the proposed pushers are problem-specific. As shown in the last section, single-particle motion in strong laser fields varies significantlywhen different particle pushers are used, even when prescribed (analytical) fields are used for thelaser. In this section, we will show that the collective behavior of a particle bunch can also varysignificantly depending on the choice of the pusher unless very small time steps are used. Wehave implemented the proposed particle pusher into
Osiris . The aforementioned P3 is adoptedbecause it uses a time-staggering layout for the particle positions and momenta, along with theresulting need for a leapfrog advance of the particle position. In the full 2D PIC simulations, a19
100 200 300 400 500 600 700 800-20002000 100 200 300 400 500 600 700 800-20-1001020300 100 200 300 400 500 600 700 800-1-0.500.51 P1P2P3P4P5P6RK4 0 100 200 300 400 500 600 700 800-20002000 100 200 300 400 500 600 700 800-20-1001020300 100 200 300 400 500 600 700 800-1-0.500.51 P1P2P3P4P5P6RK4 (a)(b)(c) (d)(e)(f)
Figure 2: Single-particle motion in a head-on collision with an ultra-intense laser ( a = 300) using various numericalschemes. The test particle has an initial longitudinal momentum p = − m e c . Evolution of (a)(d) transversemomentum p , (b)(e) phase in laser field and (c)(f) transverse spin s are compared for two values of time step. Thefield felt by the test particle is determined analytically. All proposed pushers give good agreement for ∆ t = 0 . ω − ,whereas P5 and P6 give the best agreement for ∆ t = 0 . ω − . -20002000 1 2 3 4 5 10 -35-30-25-20 P1P2 P3P4P5P6RK40 1 2 3 4 5 10 -1-0.500.51 0 1 2 3 4 5 10 -20002000 1 2 3 4 5 10 -35-30-25-20 P1P2 P3P4P5P6RK40 1 2 3 4 5 10 -1-0.500.51 (a)(b)(c) (d)(e)(f) Figure 3: Single-particle motion in an ultra-intense laser ( a = 300) using various numerical schemes. The test particleis initialized at rest. Evolution of (a)(d) transverse momentum p , (b)(e) phase in laser field and (c)(f) transverse spin s are compared for two values of time step. The field felt by the test particle is determined analytically. All proposedpushers give good agreement for ∆ t = 0 . ω − , whereas P5 and P6 give the best agreement for ∆ t = 0 . ω − . -100-50050100 P1 -100-50050100 P2 -100-50050100 P3 -100-50050100 P4 (a) (b)(c) (d) -100-50050100 P5 -100-50050100 P6 (e) (f) Figure 4: Evolution of transverse momentum p in a plane-wave ultra-intense laser ( a = 100) using various numericalschemes and values of time step. The particle is initialized at rest where the vector potential of the laser is zero. Theblack dashed line represents the result of using fourth-order Runge-Kutta method with very high precision, and thetime step is reduced for each scheme until convergence is reached. µ m wavelength plane-wave laser pulse with a = 500 and 50 ω − FWHM duration for the fieldenvelope collides head-on with an electron beam. The electron beam has a bi-Gaussian densitydistribution with rms radius σ = 10 k − , rms length σ = 15 k − (ˆ1 is the propagation direction)and initial momentum p = − m e c . The beam has zero emittance and energy spread. The cellsizes are ∆ x = 0 . k − and ∆ x = 2 k − , and the time step is ∆ t = 0 . ω − . To accurately simulatethe particle motion in the laser field, we have used a Maxwell solver with an extended stencil [42]to reduce the numerical errors arising from numerical dispersion and the interlacing of E and B fields in time.Figure 5 shows the laser field and beam density distribution. As shown in Fig. 5(a), theelectron beam initially moves toward the laser pulse from right to left. The bunch length isthen compressed by the extremely strong radiation pressure of the leading edge of the laser. Thepropagation direction of the beam is eventually reversed so that it co-moves with the laser pulseas shown in Fig. 5(b). There are significant differences in phase space between the “standard” andthe proposed numerical schemes. Figure 6 shows the x - p - p space phase for P1 [Figs. 6(a) and(c)] and P3 [Figs. 6(b) and (d)] at t = 60 ω − and t = 200 ω − . At t = 60 ω − there are only slightdifferences between the two schemes. At t = 200 ω − the phase space distribution begins to broadenfor P1 while it remains narrow for P3. Since the macro-particles have negligible charge (they aretest particles) and thus weak particle-particle interactions, the physical quantities including p and p should only be a function of x for a given time t in the 1D limit (plane-wave laser). Therefore,the narrow distribution in Fig. 6(d) is to be expected. We also conducted convergence tests usingP1 with ten-fold higher resolution in space and time. The results converged with those shown inFig. 6(d) for the larger time step using P3. -100 1000-40400 -100 1000-40400 (a) (b) laser beam laser beam Figure 5: 2D PIC simulation results using
Osiris . Snapshots of the beam density (green) and laser field (red andblue) are shown at (a) t = 0 and (b) t = 38 ω − . The bunch length is compressed and reversed by the extremelystrong radiation pressure of the laser. The scheme P3 is used to generate the plot. We also compared how the evolution of the spin precession is modified for the different schemes.The spin of the electron beam is initially polarized along the positive ˆ1-direction with a smalldivergence, as shown in Fig. 7(a). In Fig. 7 we plot the spin in the rest frame so that all theparticles move on the surface of a sphere of radius (cid:126) / s - s - s space. When the beam startsto interact with the laser field, the particles move down toward the negative ˆ1-direction along thelongitudes. Significant differences between the schemes can be seen at t = 200 ω − : the particlesadvanced analytically in momentum space and with the exact spin pusher in spin space using P3[see Fig. 7(c)] are fully congregated at the pole in the negative ˆ1-direction, while the particles23 (a) (b)(c) (d)P1 P3 Figure 6: Particle distributions in x - p - p phase space for schemes (a)(c) P1 and (b)(d) P3 at two different times.The phase space distribution remains narrow late in time when using P3, as is expected. -110 101 0-1 -1 -110 101 0-1 -1-110 101 0-1 -1 (a) (b) (c) Figure 7: Particle distributions in s - s - s space at (a) t = 0 and (b)(c) t = 200 ω − . P1 and P3 were used to obtainthe results in (b) and (c), respectively.
7. Performance
In this section, we compare the performance of the P3 and P4 implementations into
Osiris against each other and against the standard Boris push. We carried out two-dimensional simu-lations of a thermal plasma with 4 macro-particles per cell and fixed ions. The particle shapescorresponded to linear weighting/interpolation. The box was 512 ×
512 cells large, and the cellsize was 0.3 Debye lengths in each direction.
Osiris was compiled using GNU Fortran 8.3 with-O3 optimization on an Intel Core i7-8650U processor. For each time step, the time cost of var-ious procedures including the momentum advance, RR correction and spin advance, along withposition advances, field interpolation and current update (“other”) are summarized in Fig. 8. Thecomputational cost of the momentum advance in P3 is 2.4 times that of P1 (Boris scheme), andthe spin advance in P3 is 5.4 times slower than that of P1 (Vieira scheme). With the RR correctionincluded, the computational cost of the momentum advance of P3 is 1.7 times that of P1. P4 isslightly faster than P3, being only 1.5 times slower than P1. The additional cost of updating thepositions, interpolating fields from the grid onto particle positions and depositing the current ontothe grid is shown by the yellow blocks in Fig. 8.In summary, for simulations where particle spin is not considered, the schemes employingexact momentum pushers can provide competitive performance on a per-time-step basis comparedto schemes using regular pushers; this includes weak/moderate field scenarios where the regularpushers remain accurate at conventionally selected time steps. However, in moderate- to strong-field regimes, time steps 10–100 times smaller are required for regular pushers to obtain the accuracyof the analytical pushers. Therefore, these new pushers can significantly reduce the computationaltime needed for high-fidelity simulations. We note that the performance differences will becomeeven smaller as the order of the particle shape increases, since the steps encapsulated in yellow arethe same across each scheme and will take longer per particle.25 time / part. / step (ns)
P1P3P4 u adv.RR corr.others adv.
Figure 8: Performance of various pushers implemented in
Osiris . The “other” category includes field interpolation,update of particle positions and current deposition. The proposed pushers can provide competitive performance withstandard pushers when spin is not considered.
8. Conclusion
In this article, we derived the analytic solutions to the change in the four-vector of momentumand position while including a reduced form of the radiation reaction (RR). When the equationsof motion are written in covariant form, analytic solutions can be found straightforwardly if theelectric and magnetic fields are considered constant (in both space and time) over a single timestep. We obtained forms of the solutions to both the momentum (proper velocity) and position[( x , u ) phase space] using projection operators amenable to PIC codes. The trajectory of ( x , u )can be accurately computed in the strong-field regime with these explicit, closed-form expressionsusing much larger time steps than would be required for standard pushers. These expressions areanalytic, so any errors arise only from the assumption of constant and uniform fields at each timestep. When the RR is involved, these expressions are still highly accurate, except for in cases whereclassical theory fails.With an analytical solution to u and keeping the fields constant and uniform, the Bargmann-Michel-Telegdi equation can also be analytically solved, and the closed-form solutions can be usedto simulate spin precession in strong fields. Although these expressions are only perfectly accuratewithout RR, the effect can still be properly taken into account by separately including radiativeimpulse corrections to u . This semi-analytical approach can also be used when RR is modeled asa QED process.The advantage in computational efficiency (defined as the computational time to accuratesolution) of the proposed 9D phase space pusher over existing schemes was demonstrated througha series of single-particle simulations where the fields are associated with a laser. It is shown thatthe proposed pusher can yield correct or sufficiently accurate phase space trajectories with timesteps an order of magnitude smaller than for the standard split operator pushers for normalizedlaser amplitudes a on the order of at least 10 . We note that for problems where the fields varyslowly in time (including high-amplitude imposed magnetic fields), the proposed pusher will beeven more efficient than standard schemes. For example, when the laser fields are known (given)such that the position can also be analytically updated, the full analytic pusher can in some casesobtain accurate results for ω ∆ t = 4 while the standard split operator pushers require ω ∆ t = . Osiris , main-taining the leapfrog position advance. Therefore, only the momentum update needed to be modi-fied while the field solver, position update and current deposit remained unchanged. Using
Osiris ,PIC simulations were also conducted to compare the proposed numerical scheme against standardschemes for the head-on collision of a spin polarized electron beam with an ultra-intense laserpulse. The results showed significant differences in the phase space (including the spin precession)between the proposed and the standard schemes. As the time step was reduced, the standardpusher simulations converged to that of the analytical pusher case with the larger time steps. Al-though these sample simulations were all conducted in the context of laser-plasma interactions, theproposed algorithm itself is general and can be applied to many other research fields.Future work may involve the development of PIC algorithms that define the position and mo-mentum at the same time or that use predictor-corrector algorithms. We found that the proposedscheme without (with) the spin advance is only 20 (80) percent slower per particle than standardpushers (including field interpolation, momentum update, and current deposit) for linear parti-cle shapes. However, the proposed scheme can provide accurate solutions with time steps muchlarger than those required for standard pushers (depending on the field strength and configuration),generating significant speedups. For example, for some of the laser-plasma interaction examplespresented here where the laser fields are updated using the field solver, time steps as much as 10times larger can be used with the proposed scheme.
Acknowledgments
This work was supported in parts by the US Department of Energy contract number de-sc0010064, de-sc0019010 and SciDAC FNAL subcontract 644405, Lawrence Livermore NationalLaboratory subcontract B634451, and US National Science Foundation grant numbers 1806046.The work of MV was supported by the European Research Council (ERC-2015-AdG Grant No.695088) and Portuguese Science Foundation (FCT) Grant No. SFRH/BPD/119642/2016. Simu-lations were carried out on the Cori Cluster of the National Energy Research Scientific ComputingCenter (NERSC).
Appendix A. Eigensystem of the field tensor F and relevant properties The eigenvalues of the field tensor F (under the assumption that the elements are constantin τ ) are determined by the characteristic equation det( F − λI ) = 0. This leads directly to thefollowing equations for the eigenvalues, λ − I λ − I = 0 , (A.1)where I and I are the well-known Lorentz invariants [28], I = | E | − | B | , I = E · B . (A.2)From this it follows that there are two pairs of eigenvalues, λ = ± κ and λ = ± i ω , where κ = 1 √ (cid:114) I + (cid:113) I + 4 I , ω = 1 √ (cid:114) −I + (cid:113) I + 4 I . (A.3)27o facilitate the derivations of the analytic pushers, we introduce two subspaces that are definedby the eigenvectors, i.e., S κ = span { e κ , e − κ } and S ω = span { e i ω , e − i ω } , where e λ denotes theeigenvector associated with the eigenvalue λ . For general four-vectors V κ ∈ S κ and V ω ∈ S ω , thefollowing relations, F V κ = κ V κ , F V ω = − ω V ω (A.4)are satisfied. These relations can be easliy verified by expressing V κ and V ω as a linear combinationof the appropriate eigenvectors and then using the fact that F e λ = λe λ . To decompose an arbitraryfour-vector V into S κ and S ω subspaces, V is rewritten as V = V κ + V ω , and then the operator F is applied to both sides. These two equations can then be solved for V κ and V ω as V κ = ( ω V + F V ) / ( κ + ω ) , V ω = ( κ V − F V ) / ( κ + ω ) , (A.5)which indicates that P κ ≡ ( κ + ω ) − ( ω I + F ) and P ω ≡ ( κ + ω ) − ( κ I − F ) are the projectionoperators of a four-vector into S κ and S ω .According to the Cayley-Hamilton theorem [43], the field tensor F also satisfies the character-istic equation (A.1), i.e., F − I F − I I = 0, which leads to( κ I − F )( ω I + F ) = 0 . (A.6)With this property, we can prove that S κ and S ω are mutually orthogonal by explicitly taking theinner product of V κ and V ω and then substituting in Eq. (A.5) to give( κ + ω ) ( V κ | V ω ) = V T ( ω I + ( F ) T ) G ( κ I − F ) V = V T G ( ω I + F )( κ I − F ) V = 0 . (A.7)We have also used Eq. (A.6) and an obvious relation between the field tensor and its transpose, F = − GF T G .Another important relation that will be frequently used in this article is that F = 0 when I = I = 0 (or κ = ω = 0). This can be shown by explicit calculation using Eq. (2.2) to give F = I F + I F ∗ , (A.8)where F ∗ is the dual tensor defined as F ∗ = B B B B − E E B E − E B − E E . (A.9)Therefore, F vanishes when I = 0 and I = 0. Appendix B. Modulus of four-velocity
In this appendix, we will discuss the nature of the modulus of the four-velocity components in S κ and S ω . As addressed in Appendix A, the subspace components u κ and u ω can be obtained byprojecting u to the subspaces using the projection operators P κ and P ω , i.e., u κ = ω u + F uκ + ω , u ω = κ u − F uκ + ω . (B.1)28ombining this with the characteristic equation for F [Eq. (A.6)] written as F = ( κ − ω ) F + κ ω I , the modulus of u κ and u ω can be calculated as | u κ | = ω − | F u | κ + ω , | u ω | = κ + | F u | κ + ω . (B.2)Now we will prove that the modulus of the four-force has a maximum of − κ . The problemcan be more accurately defined for a given F by finding the extrema of | F u | under the restrictedcondition | u | = 1. We use the Lagrange multiplier method to handle this problem and constructthe Lagrangian function L ( u, χ ) = | F u | + χ ( | u | − χ is the Lagrange multiplier.The extremum point ( u ∗ , χ ∗ ) is determined by ∂ u L = 0 and ∂ χ L = 0. The latter equation directlygives the restricted condition | u | = 1, and the former can be written in the matrix form as12 ∂ L ∂u = F T GF u + χGu = 0 . (B.3)The existence of a non-trivial solution for u requires det( F T GF + χG ) = 0, from which χ ∗ can bedetermined. Using the fact that F T G = − GF and G has a non-zero determinant [det( G ) = − F − χI ) = 0, which is exactly the characteristic equation of F ; the solution χ ∗ is theassociated eigenvalue. Recalling that F has two pairs of eigenvalues, ± κ and ± i ω , the characteristicequation therefore has two roots, χ ∗ = κ and χ ∗ = − ω . Noticing that ( u ∗ , χ ∗ ) satisfies Eq. (B.3),the meaning of χ ∗ can be revealed by left multiplying Eq. (B.3) by u ∗ T , which gives | F u ∗ | = − χ ∗ .This indicates that − κ and ω are two extrema of | F u | . However, the extremum ω should bediscarded because | F u | < | F u | = ( u · E ) − | γ E + u × B | = ˙ γ − | ˙ u | = ( u · ˙ u ) γ − | ˙ u | ≤ | ˙ u | (cid:18) | u | γ − (cid:19) < . Therefore, | F u | has the unique extremum − κ , and we can verify | F u | ≤ − κ by substituting inan arbitrary u . With this property we can know from Eq. (B.2) that | u κ | ≥ , | u ω | ≤ . (B.4) Appendix C. Inhomogeneous solutions to Eqs. (4.3) and (4.4)
In this appendix, we will seek the inhomogeneous solutions to Eqs. (4.3) and (4.4). Noticethat the inhomogeneous terms in Eqs. (4.3) and (4.4) contain u κ and u ω , respectively, so the trialsolutions can be constructed as ˜ s κ = C κ ( τ ) u κ + D κ ( τ ) ˙ u κ , (C.1)˜ s ω = C ω ( τ ) u ω + D ω ( τ ) ˙ u ω . (C.2)Inserting the trial solution of Eq. (C.1) back into Eq. (4.3) and comparing the coefficients of termsproportional to u κ and ˙ u κ , we get two ODEs for C κ ( τ ) and D κ ( τ ),˙ C κ ( τ ) = aκ D κ ( τ ) − af ( τ ) , (C.3)˙ D κ ( τ ) = aC κ ( τ ) . (C.4)29ubstituting Eq. (4.11) into above equations, we can find out the solutions that satisfy the zeroinitial conditions, i.e., C κ (0) = 0 and D κ (0) = 0. We note that the initial conditions of ˙ C κ and ˙ D κ required by the inhomogeneous solutions, i.e., ˙ C κ (0) = − af and ˙ D κ (0) = 0, are naturally satisfiedaccording to Eqs. (C.3) and (C.4). The solutions are given by C κ = ˙ f cos( a Ω τ ) − cosh( aκτ ) a ( κ + Ω ) + h Ω κ sin( a Ω τ ) − h κ Ω sinh( aκτ ) κ Ω( κ + Ω ) ,D κ = ˙ f k sin( a Ω τ ) − Ω sinh( aκτ ) aκ Ω( κ + Ω ) − h Ω κ cos( a Ω τ ) + h κ Ω cosh( aκτ ) κ Ω ( κ + Ω ) + I κ Ω , (C.5)where h κ ≡ I + f κ and h Ω ≡ I − f Ω .The ODEs of C ω and D ω can be similarly established by inserting Eq. (C.2) into Eq. (4.4) andcomparing the coefficients: ˙ C ω ( τ ) = − aω D ω ( τ ) − af ( τ ) , (C.6)˙ D ω ( τ ) = aC ω ( τ ) . (C.7)The solutions that satisfy C ω (0) = 0 and D ω (0) = 0 are C ω = ˙ f cos( aωτ ) − cos( a Ω τ ) a ( ω − Ω ) + h ω Ω sin( aωτ ) − h Ω ω sin( a Ω τ ) ω Ω( ω − Ω ) ,D ω = ˙ f Ω sin( aωτ ) − ω sin( a Ω τ ) aω Ω( ω − Ω ) − h ω Ω cos( aωτ ) − h Ω ω cos( a Ω τ ) ω Ω ( ω − Ω ) − I ω Ω , (C.8)where h ω ≡ I − f ω . References [1] The extreme light infrastructure (ELI), .[2] Exawatt center for extreme light studies (XCELS), https://xcels.iapras.ru/ .[3] Shanghai superintense ultrafast laser facility (SULF), http://english.siom.cas.cn/Newsroom/hotnews/201907/t20190710_212831.html .[4] T. Tajima, J. M. Dawson, Laser electron accelerator, Physical Review Letters 43 (4) (1979) 267.[5] P. Chen, J. M. Dawson, R. W. Huff, T. Katsouleas, Acceleration of electrons by the interaction of a bunchedelectron beam with a plasma, Physical Review Letters 54 (7) (1985) 693–696.[6] C. Joshi, T. Katsouleas, Plasma accelerators at the energy frontier and on tabletops, Physics Today 56 (6)(2003) 47–53.[7] W. Lu, M. Tzoufras, C. Joshi, F. S. Tsung, W. B. Mori, J. Vieira, R. A. Fonseca, L. O. Silva, Generatingmulti-gev electron bunches using single stage laser wakefield acceleration in a 3d nonlinear regime, PhysicalReview Special Topics - Accelerators and Beams 10 (6) (2007) 061301.[8] M. Vranic, T. Grismayer, R. A. Fonseca, L. O. Silva, Quantum radiation reaction in head-on laser-electron beaminteraction, New Journal of Physics 18 (7) (2016) 073035.[9] J. M. Dawson, Particle simulation of plasmas, Reviews of Modern Physics 55 (2) (1983) 403–447.[10] R. W. Hockney, J. W. Eastwood, Computer simulation using particles, crc Press, 1988.[11] C. K. Birdsall, A. B. Langdon, Plasma Physics via Computer Simulation, Institute of Physics Pub., 2005.[12] B. B. Godfrey, J.-L. Vay, Numerical stability of relativistic beam multidimensional pic simulations employingthe Esirkepov algorithm, Journal of Computational Physics 248 (2013) 33–46.[13] X. Xu, P. Yu, S. F. Martins, F. S. Tsung, V. K. Decyk, J. Vieira, R. A. Fonseca, W. Lu, L. O. Silva, W. B. Mori,Numerical instability due to relativistic plasma drift in EM-PIC simulations, Computer Physics Communications184 (11) (2013) 2503–2514.
14] P. Yu, X. Xu, V. K. Decyk, F. Fiuza, J. Vieira, F. S. Tsung, R. A. Fonseca, W. Lu, L. O. Silva, W. B. Mori,Elimination of the numerical cerenkov instability for spectral EM-PIC codes, Computer Physics Communications192 (2015) 32–47.[15] P. Yu, X. Xu, A. Tableman, V. K. Decyk, F. S. Tsung, F. Fiuza, A. Davidson, J. Vieira, R. A. Fonseca,W. Lu, L. O. Silva, W. B. Mori, Mitigation of numerical Cerenkov radiation and instability using a hybrid finitedifference-FFT Maxwell solver and a local charge conserving current deposit, Computer Physics Communications197 (2015) 144–152.[16] F. Li, P. Yu, X. Xu, F. Fiuza, V. K. Decyk, T. Dalichaouch, A. Davidson, A. Tableman, W. An, F. S. Tsung,R. A. Fonseca, W. Lu, W. B. Mori, Controlling the numerical Cerenkov instability in PIC simulations usinga customized finite difference Maxwell solver and a local FFT based current correction, Computer PhysicsCommunications 214 (2017) 6–17.[17] A. B. Langdon, Effects of the spatial grid in simulation plasmas, Journal of Computational Physics 6 (2) (1970)247–267.[18] H. Okuda, Nonphysical noises and instabilities in plasma simulation due to a spatial grid, Journal of Computa-tional Physics 10 (3) (1972) 475–486.[19] M. D. Meyers, C. K. Huang, Y. Zeng, S. A. Yi, B. J. Albright, On the numerical dispersion of electromagneticparticle-in-cell code: Finite grid instability, Journal of Computational Physics 297 (2015) 565–583.[20] C. K. Huang, Y. Zeng, Y. Wang, M. D. Meyers, S. Yi, B. J. Albright, Finite grid instability and spectral fidelityof the electrostatic particle-in-cell algorithm, Computer Physics Communications 207 (2016) 123–135.[21] X. Xu, F. Li, F. S. Tsung, T. N. Dalichaouch, W. An, H. Wen, V. K. Decyk, R. A. Fonseca, M. J. Hogan, W. B.Mori, On numerical errors to the fields surrounding a relativistically moving particle in PIC codes, Journal ofComputational Physics 413 (2020) 109451.[22] J. P. Boris, R. A. Shanny, Proceedings: Fourth Conference on Numerical Simulation of Plasmas, November 2,3, 1970, Naval Research Laboratory, 1972.[23] J.-L. Vay, Simulation of beams or plasmas crossing at relativistic velocity, Physics of Plasmas 15 (5) (2008)056701.[24] A. V. Higuera, J. R. Cary, Structure-preserving second-order integration of relativistic charged particle trajec-tories in electromagnetic fields, Physics of Plasmas 24 (5) (2017) 052104.[25] D. F. Gordon, B. Hafizi, J. Palastro, Pushing particles in extreme fields, AIP Conference Proceedings 1812 (1)(2017) 050002.[26] D. Gordon, B. Hafizi, Special unitary particle pusher for extreme fields (jun 2020). arXiv:2006.10146 .[27] J. Ptri, A relativistic particle pusher for ultra-strong electromagnetic fields (2019). arXiv:1910.04591 .[28] J. D. Jackson, Classical electrodynamics, John Wiley & Sons, 2007.[29] L. D. Landau, The classical theory of fields, Vol. 2, Elsevier, 2013.[30] H. Spohn, The critical manifold of the lorentz-dirac equation, EPL (Europhysics Letters) 50 (3) (2000) 287.[31] M. Vranic, J. L. Martins, R. A. Fonseca, L. O. Silva, Classical radiation reaction in particle-in-cell simulations,Computer Physics Communications 204 (2016) 141–151.[32] A. Ilderton, G. Torgrimsson, Radiation reaction in strong field qed, Physics Letters B 725 (4-5) (2013) 481–486.[33] A. V. Arefiev, G. E. Cochran, D. W. Schumacher, A. P. L. Robinson, G. Chen, Temporal resolution criterion forcorrectly simulating relativistic electron motion in a high-intensity laser field, Physics of Plasmas 22 (1) (2015)013103.[34] M. Tamburini, F. Pegoraro, A. D. Piazza, C. H. Keitel, A. Macchi, Radiation reaction effects on radiationpressure acceleration, New Journal of Physics 12 (12) (2010) 123005.[35] Y. F. Li, R. Shaisultanov, K. Z. Hatsagortsyan, F. Wan, C. H. Keitel, J. X. Li, Ultrarelativistic electron-beampolarization in single-shot interaction with an ultraintense laser pulse, Physical Review Letters 122 (15) (2019)154801.[36] H.-H. Song, W.-M. Wang, J.-X. Li, Y.-F. Li, Y.-T. Li, Spin-polarization effects of an ultrarelativistic electronbeam in an ultraintense two-color laser pulse, Physical Review A 100 (3) (2019).[37] X. S. Geng, L. L. Ji, B. F. Shen, B. Feng, Z. Guo, Q. Q. Han, C. Y. Qin, N. W. Wang, W. Q. Wang, Y. T. Wu,X. Yan, Q. Yu, L. G. Zhang, Z. Z. Xu, Spin-dependent radiative deflection in the quantum radiation-reactionregime, New Journal of Physics 22 (1) (2020) 013007.[38] J. Vieira, C. K. Huang, W. B. Mori, L. O. Silva, Polarized beam conditioning in plasma based acceleration,Physical Review Special Topics - Accelerators and Beams 14 (7) (2011).[39] 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,et al., OSIRIS: A three-dimensional, fully relativistic particle in cell code for modeling plasma based accelerators,in: International Conference on Computational Science, Springer, 2002, pp. 342–351.
40] R. G. Hemker, Particle-in-cell modeling of plasma-based accelerators in two and three dimensions, arXiv preprintarXiv:1503.00276 (2015).[41] V. Bargmann, L. Michel, V. L. Telegdi, Precession of the polarization of particles moving in a homogeneouselectromagnetic field, Physical Review Letters 2 (1959) 435–436.[42] F. Li, K. G. Miller, X. Xu, F. S. Tsung, V. K. Decyk, W. An, R. A. Fonseca, W. B. Mori, A new fieldsolver for modeling of relativistic particle-laser interactions using the particle-in-cell algorithm, arXiv preprintarXiv:2004.03754 (2020).[43] A. S. Householder, The theory of matrices in numerical analysis, Courier Corporation, 2013.40] R. G. Hemker, Particle-in-cell modeling of plasma-based accelerators in two and three dimensions, arXiv preprintarXiv:1503.00276 (2015).[41] V. Bargmann, L. Michel, V. L. Telegdi, Precession of the polarization of particles moving in a homogeneouselectromagnetic field, Physical Review Letters 2 (1959) 435–436.[42] F. Li, K. G. Miller, X. Xu, F. S. Tsung, V. K. Decyk, W. An, R. A. Fonseca, W. B. Mori, A new fieldsolver for modeling of relativistic particle-laser interactions using the particle-in-cell algorithm, arXiv preprintarXiv:2004.03754 (2020).[43] A. S. Householder, The theory of matrices in numerical analysis, Courier Corporation, 2013.