Particle integrator for particle-in-cell simulations of ultra-high intensity laser-plasma interactions
PParticle integrator for particle-in-cell simulationsof ultra-high intensity laser-plasma interactions
Kavin Tangtartharakul, Guangye Chen, and Alexey Arefiev Department of Mechanical and Aerospace Engineering, University of California San Diego, La Jolla, CA 92093,USA Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA (Dated: 23 February 2021)
Particle-in-cell codes are the most widely used simulation tools for kinetic studies of ultra-intense laser-plasmainteractions. Using the motion of a single electron in a plane electromagnetic wave as a benchmark problem,we show surprising deterioration of the numerical accuracy of the PIC algorithm with increasing normalizedwave amplitude for typical time-step and grid sizes. Two significant sources of errors are identified: strongacceleration near stopping points and the temporal field interpolation. We propose adaptive electron sub-cycling coupled with a third order temporal interpolation of the magnetic field and electric field as an efficientremedy that dramatically improves the accuracy of the particle integrator.
I. INTRODUCTION
The ELI-NP laser facility [1] has recently demonstratedthat their laser system is able to achieve the projectedlaser power of 10 PW [2]. Even based on conservative es-timates, such a powerful laser pulse will be able to achieveintensities exceeding 5 × watt/cm . The correspond-ing normalized laser amplitude a ≡ | e | E /mcω (1)would exceed a ≈ E is the peak amplitudeof the laser electric field, ω is the laser frequency, c is thevacuum speed of light, and m and e are the electron massand charge. We define ω = 2 πc/λ , where λ = 820 nmis the vacuum wavelength of the laser. The significanceof reaching such high values of a is that the motion oflaser-irradiated electrons inside a target would becomeultra-relativistic.Simulations of laser-matter interactions at a > a (cid:29) λ . In this regime,the electron experiences alternating acceleration and de-celeration while periodically coming to a complete stopduring its predominantly forward motion. Even for givenfields with an analytical expression, the Boris algorithmfails to recover the particle dynamics near the stoppingpoints if a is increased for a fixed time-step ∆ t [5]. Thisleads to significant errors along the rest of the electrontrajectory and impacts the maximum energy gain. One needs to dramatically reduce the time-step ∆ t ina PIC algorithm with a (cid:29) x/c ∆ t close tounity, where ∆ x is the cell size. Otherwise, the numericaldispersion errors for electromagnetic waves representingthe laser become significant and this can adversely affectthe electron acceleration as well. For example, 80 cellsper wavelength were used to correctly reproduce electronacceleration at a ≈ . a (cid:29) a [5]. Two key components are used in com-bination to improve the accuracy: 1) adaptive temporalsub-cycling of the conventional particle pusher employ-ing the Boris algorithm [3]; 2) third-order (or higher),rather than the conventional first-order, temporal in-terpolation of electric and magnetic fields used by thepusher. Our key finding is that the Boris algorithm usingsub-cycling [5] may perform poorly when implementedinto a standard PIC algorithm. The electromagneticfields must be temporally interpolated with high-orderpolynomials for the sub-cycling to yield significant bene-fits. When using appropriately high-order temporal fieldinterpolations, the sub-cycling improves simulation per-formance (see Fig. 2) at a (cid:29) a (cid:29) a r X i v : . [ phy s i c s . p l a s m - ph ] F e b FIG. 1: Trajectories of an electron irradiated by a wave with a = 100 (left column) and the corresponding momentum (rightcolumn). The gray scale curves are the analytical solution. The color-coded curves are results from perfectly propagating PICalgorithm with particle sub-cycling ( λ/ ∆ tc = T / ∆ t = 75, λ/ ∆ x ≈ . c ∆ t/ ∆ x = 1 / .
01 = 0 .
99, and Ψ max = 0 . st order temporal interpolation of the fields; panels (c) and (d) are obtained using 3 rd ordertemporal interpolation of the fields. The black point represents the stopping point. The color gradient represents themagnitude of the E and B fields. to clearly identify the impact of temporal interpolation.In Section IV, we show that implementing particle sub-cycling into a PIC algorithm that utilizes linear temporalinterpolation is not sufficient to dramatically improve thesimulation. We then explain, in Section V, why the lin-ear interpolation is responsible for the failure of the sub-cycling algorithm. In Section VI, we describe a higherorder temporal interpolation procedure for the fields andshow that the third order interpolation dramatically im-proves the results for the model problem when used to-gether with the sub-cycling. Finally, in Section VII, wequantify the performance of the improved PIC algorithmby comparing it with that of a conventional PIC algo-rithm for the model problem. II. DESCRIPTION OF MODEL PROBLEM
In this section we review the solution of a well-knownproblem where a single electron is irradiated in vacuumby a linearly polarized electromagnetic wave [10] . Thisproblem provides an important test case for evaluatingthe performance of particle integrators at a (cid:29) x -axis with the wave electric fielddirected along the y -axis. The field evolution in this caseis fully described by a normalized vector potential thatonly has a y component and that is only a function of thephase variable ξ ∝ t − x/c . We denote this componentas a , such that E y = B z = − mωc | e | dadξ , (2) a ( ξ ) ≡ a exp (cid:2) − ( ξ − ξ ) / (2 σ ) (cid:3) sin( ξ ) , (3) ξ ≡ ω ( t − x/c ) , (4)where t is the time and a is the peak amplitude. Inthe specific examples that follow we set ξ = − π and σ = 8 π . Note that the phase velocity in this case is equalto the speed of light, v ph = c . The group velocity is alsoequal to c , such that the pulse preserves its shape duringpropagation.An analytical solution exists for an electron irradiatedby the described wave. Appendix A provides details onhow to obtain a general solution, so here we only sum-marize the key results. We take an electron that is atrest prior to the arrival of the wave, i.e. at a = 0. Thenthe electron momentum in the wave is given by p x /mc = a ( ξ ) / ≥ , (5) p y /mc = a ( ξ ) . (6)The relativistic γ -factor defined as γ = (cid:112) p /m c is γ = 1 + a ( ξ ) / . (7)As the electron moves in the wave, the following quantityremains conserved R ≡ γ − p x /mc = 1 . (8)It can be shown that R = γω dξdt , (9)so R represents the electron dephasing rate in the elec-tron’s instantaneous rest frame [11].An example of the described analytical solution inEqs. (5) and (6) is shown by the gray curve in Fig. 1bfor a = 100. The corresponding electron trajectory isshown by the gray curve in Fig. 1a, where we set the elec-tron initial location as x = y = 0. These figures illustratethat p x dominates the trajectory at a (cid:29)
1, which leadsto a large displacement in the x direction. In Fig. 1,numerical solutions, that will be discussed in the nextsection, are plotted on top of the analytic solution withthe color-coding indicating the strength of the electric, | E y | , and magnetic, | B z | , fields at the electron location.The electron reaches its most energetic state when the E and B fields vanish, which corresponds to | a ( ξ ) | = a .And when the field strength is at a local maximum, theelectron reaches a motionless position. We call these po-sitions the stopping points. We single out these pointsbecause the electron dynamics in their vicinity during nu-merical integration has strong impact on the subsequentelectron acceleration. III. PIC ALGORITHM AND TEMPORALDISCRETIZATION
The focus of the paper is on simulations of electron ac-celeration by an ultra-high intensity laser pulse and theaccuracy improvements that can be achieved within theframework of a PIC algorithm. The details of a stan-dard PIC algorithm are well-known and can be found inRef. [3]. Here we review the key aspects of the temporaldiscretization with a fixed time-step ∆ t to provide thenecessary context for the discussion that follows.The model problem detailed in Sec. II can be simulatedusing a one-dimensional (1D) version of a standard PICcode. The difficulties mentioned in Sec. I can already beobserved in 1D, so we limit our analysis to a 1D version ofa PIC code in order to make the discussion more compact.In a standard PIC algorithm, there are two types ofquantities: those known at t = n ∆ t and denoted by n ; and those known at t = ( n + 1 / t and denoted by n + 1 /
2, where n is an integer (the simulation starts at n = 0). The field equations are discretized in space andtime using the Yee scheme. In our problem, there areonly E y and B z field components, so we have E y | i +1 / n = E y | i +1 / n − − c ∆ t ∆ x (cid:16) B z | i +1 n − / − B z | in − / (cid:17) , (10) B z | in +1 / = B z | in − / − c ∆ t ∆ x (cid:16) E y | i +1 / n − E y | i − / n (cid:17) . (11)Here the superscript represents the spatial discretizationwith a cell size ∆ x , such that E y | i +1 / is the field at x = x i +1 / ≡ ( i + 1 / x and B z | i is the field at x = x i ≡ i ∆ x . The Boris particle pusher that we use toadvance the considered electron requires the knowledgeof B z at t = n ∆ t . The standard approach, which in thispaper is called “standard PIC algorithm”, uses a lineartemporal interpolation: B z | in = 12 (cid:16) B z | in +1 / + B z | in − / (cid:17) . (12)A shape function S is used to obtain the fields E y and B z at the electron position x n (we denote these fields witha superscript p ): E py | n = (cid:88) i E y | i +1 / n S ( x i +1 / − x n ) , (13) B pz | n = (cid:88) i B z | in S ( x i − x n ) , (14)In this paper, we employ a triangular or tent shape func-tion [3]). The electron is advanced using the standardBoris particle pusher detailed below: p n +1 / = p + + q E p | n ∆ t , (15) p + = p − + p (cid:48) × s , (16) p − = p n − / + q E p | n ∆ t , (17) p (cid:48) = p − + p − × Ψ , (18) s = 2 Ψ , (19) Ψ = q B p | n ∆ t γ n m e c , (20) γ n = (cid:115) (cid:18) p − m e c (cid:19) , (21) γ n +1 / = (cid:115) (cid:18) p n +1 / m e c (cid:19) , (22) v n +1 / = p n +1 / m e γ n +1 / , (23) r n +1 = r n + ∆ t ( v n +1 / ) . (24)In the standard PIC algorithm, the temporal dis-cretization of the fields can impact the electron dynamicsin two distinct ways: through numerical dispersion and FIG. 2: Relative dephasing, (cid:15) dephasing , and energy, (cid:15) energy , errors for different versions of the PIC algorithm over a range oflaser amplitudes a . Panel (a) is obtained using the standard PIC with and without the sub-cycling. Panel (b) is obtainedusing the perfectly propagating PIC with and without the sub-cycling. The data sets shown with red and green markers inboth panels are obtained by replacing the 1 st order interpolation with the 3 rd order temporal interpolation for the fieldsdetailed in Sec. VI. All simulations use ∆ t = λ/ c and ∆ x ≈ λ/ .
26 to satisfy c ∆ t/ ∆ x = 0 . through errors introduced by the temporal interpolation.It is well-understood that the numerical dispersion of theYee scheme alters the electron dephasing rate and grad-ually distorts the envelope of the considered pulse. Itsnegative effects can be greatly reduced by setting the ra-tio ∆ x/c ∆ t close to unity in 1D simulations without dra-matically increasing the resolution. In contrast to that,the errors introduced by linearly interpolating the fieldscan only be reduced by reducing the time-step ∆ t .In order to clearly identify the impact of the tempo-ral interpolation on electron dynamics, we introduce areduced version of the PIC algorithm where the fields ofthe laser propagate according to the analytical expres-sions given by Eqs. (2) - (4), using the same ∆ t and ∆ x as the standard PIC algorithm. As a result, it does notsolve the Maxwell’s equations, but retains the need forthe temporal and spatial interpolations of the field fromgrid points to particles. In what follows, we refer to thisalgorithm as the “perfectly propagating PIC”. The algo-rithm is particularly useful for simulating the describedmodel problem because only the fields of the laser needto be accounted for in this case and this is exactly whatis done in the reduced algorithm. IV. PERFORMANCE OF PIC ALGORITHMS WITHFIRST-ORDER TEMPORAL INTERPOLATION
In this section, we examine how the standard and theperfectly propagating 1D PIC algorithms perform withthe increase of a when employing the first-order tem-poral field interpolation. We also introduce the conceptof temporal sub-cycling and evaluate its impact on theperformance of both algorithms.A robust metric for evaluating the performance of aPIC algorithm is the maximum relative error in the de-sphasing rate, (cid:15) dephasing ≡ (cid:12)(cid:12)(cid:12)(cid:12) R − R th R th (cid:12)(cid:12)(cid:12)(cid:12) max = | R − | max . (25)Here R th = 1 is the analytical result for the model prob-lem and R = γ − p x /mc (26)is the dephasing calculated using the numerical valuesof the electron momentum. The advantage of using (cid:15) dephasing is that it is straightforward to compute regard-less of the numerical algorithm.Fig. 2a shows that (cid:15) dephasing increases almost mono-tonically with the increase of a in a simulation that usesthe standard 1D PIC algorithm (purple points). Thetime-step is set at ∆ t = T /
75, where T = 2 π/ω isthe period of the considered wave. We keep the ratio c ∆ t/ ∆ x = 0 .
99 fixed as we vary a to ensure that thenumerical dispersion remains unaffected by the scan. Ac-cordingly, the grid size is set set to ∆ x ≈ λ/ .
26. Thisresolution is similar to that used in Ref. [8] to correctlyreproduce electron acceleration at a ≈ .
5. At a > (cid:15) dephasing are very similar to those shown in Fig. 2a. Wehave additionally examined how the energy gain deviatesfrom the analytical solution along the electron trajectory.The corresponding metric is the maximum relative errorin the energy gain defined as (cid:15) energy ≡ (cid:12)(cid:12)(cid:12)(cid:12) γ − γ th γ th (cid:12)(cid:12)(cid:12)(cid:12) max , (27)where γ th is the analytical result given by Eq. (7). Thiscomparison is possible due to the fact that the perfectlypropagating PIC uses analytically prescribed fields andthe theoretical value of γ is known for these fields. Fig. 2bshows that (cid:15) energy also increases almost monotonicallywith a (purple points). At a >
40, the errors in the en-ergy gain become significant, with the maximum energyunderestimated by a factor of three at a = 100.The results for the perfectly propagating PIC rule outthe numerical wave propagation as the primary cause forthe poor performance at a (cid:29)
1. The observed errorsare clearly associated with the particle integrator. It waspreviously shown that the Boris algorithm itself performspoorly even when the fields acting on the particle are pre-scribed using the analytical form without any temporalinterpolation [5]. In this case, the analytical solution pro-vides E n and B n directly for the Boris algorithm. It wasfound that the errors originate at the stopping points ofthe trajectory where large rotations of the momentumare performed in a single time-step [5]. It was shownthat the errors in the energy gain can be dramaticallyreduced by sub-cycling the Boris algorithm’s time-stepsnear the stopping points [5].In an attempt to improve the performance of the con-sidered PIC algorithms, we first implement a similar sub-cycling algorithm as described in Ref. [5]. The details ofimplementation are shown as Algorithm 1 below. We de-fine a maximum rotation angle, Ψ max , that the particle’smomentum is allowed to experience during that singletime-step. Conveniently, the Boris algorithm splits theevolution of the particle’s momentum into three steps:half of an electric field push, a magnetic field rotation,and another half of an electric field push. We can thus use that rotation generated by the magnetic field to de-termine the size of the sub-cycling time-step ∆ t ∗ for theBoris algorithm. The momentum rotation over the time-step ∆ t ∗ is Ψ = − | e | B ∆ t ∗ γmc . (28)We set ∆ t ∗ = ∆ t/ k where k is the smallest whole num-ber that satisfies the condition | Ψ | < Ψ max , with ∆ t being the original time-step and also the time-step forthe Yee scheme.The switch to a new time-step must be consistent withthe leapfrog scheme used to advance the particle’s po-sition in the standard PIC algorithm. In order to keepthe momentum and position of the sub-cycled particlestaggered by half a time-step, the first momentum pushwhen changing ∆ t ∗ from ∆ t ∗ old to ∆ t ∗ new must be per-formed using ∆ t ∗ = (∆ t ∗ old + ∆ t ∗ new ) / t ∗ = ∆ t ∗ new . This is a gen-eral approach, so it equally applies to those cases where∆ t ∗ old = ∆ t or ∆ t ∗ new = ∆ t . Algorithm 1:
Sub-Cycling ∆ t ∗ old = ∆ t ∗ new (∆ t ∗ new from previous ∆ t time-step)∆ t ∗ new = ∆ t while | Ψ | > Ψ max do ∆ t ∗ new = ∆ t ∗ new /
4Ψ = −| e | B ∆ t ∗ new / γmc ; counter = (∆ t ∗ old / ∆ t ∗ new ) for j = 1 : counter do ∆ t ∗ = (∆ t ∗ old + ∆ t ∗ new ) / p j +1 / from Eq. (15)using ∆ t ∗ instead of ∆ t ∆ t ∗ = ∆ t ∗ new find x j +1 from Eqs. (22) - (24)using ∆ t ∗ instead of ∆ t ∆ t ∗ old = ∆ t ∗ new In order to synchronize the motion of sub-cycled parti-cles with the field discretization, we set, as stated earlier,∆ t ∗ = ∆ t/ k , where k can be different for different parti-cles. The temporal and spatial step sizes for the fields arefixed to ∆ t and ∆ x for the duration of the entire simula-tion. The sub-cycling algorithm is used only to advanceparticles, with the | Ψ | < Ψ max condition applied to eachelectron individually to determine the corresponding ∆ t ∗ .The condition is only applied once at the beginning of ev-ery ∆ t time-step to determine the value of ∆ t ∗ for thewhole ∆ t time-step. Since the fields are not known atthe sub-cycled time-steps ∆ t ∗ , we use a linear interpola-tion procedure for both the electric and magnetic fields.Specifically, we interpolate the electric and magnetic fieldfields to t + ∆ t ∗ / t to t + ∆ t ∗ . This allows the Boris algorithm toupdate the momentum using field values halfway of each∆ t ∗ . Spatially, the fields are calculated by the particle’s FIG. 3: Dephasing error along the electron trajectory in themodel problem at a = 50 for different temporal fieldinterpolation orders [1 st (a), 2 nd (b), and 3 rd (c)] inperfectly propagating PIC with and without the sub-cycling.All simulations use ∆ t = λ/ c and ∆ x ≈ λ/ .
52 tosatisfy c ∆ t/ ∆ x = 0 .
99. The maximum rotation angle Ψ max for the sub-cycling is shown in parentheses. triangular shape function. Once the entire time inter-val ∆ t has been sub-cycled for each electron, the fieldsmust be updated and individual time-steps ∆ t ∗ can bechanged again to satisfy the new | Ψ | < Ψ max rotationcondition.Fig. 2 shows the results from the standard and per-fectly propagating PIC algorithm that employ the de-scribed sub-cycling with Ψ max = 0 .
01 (blue markers).The time-step ∆ t and the grid-size ∆ x are the same asthe ones used for the simulations without the sub-cycling(purple markers). It is unclear if the performance is bet- ter when using the sub-cycling. Figures 1a and 1b il-lustrate the electron trajectory and momentum obtainedusing the standard PIC with sub-cycling at a = 100.The trajectory is symmetric, but it is greatly distorted,whereas the maximum energy gain is roughly three timeslower than what it should be according to the analyticalsolution. The fact that the sub-cycling delivers a dra-matic improvement only in the case when the fields areprescribed analytically at the electron location and at theexact time of the momentum update [5] strongly suggeststhat there are additional sources of error associated withthe temporal interpolation. V. IMPACT OF TEMPORAL INTERPOLATION ONTHE PARTICLE INTEGRATOR
To start this section, we look to qualitatively under-stand why the use of the sub-cycling algorithm deliversonly marginal improvement. We examine the time evo-lution of the dephasing rate during a single simulation ofthe model problem with the perfectly propagating PIC.In Fig. 3a, the relative dephasing error, ( R − R th ) /R th ,is shown for every time-step. The purple curve repre-sents a simulation without the sub-cycling, whereas theblue curve represents a simulation with the sub-cyclingusing Ψ max = 0 .
01. The sharp downward spikes alongthese curves coincide with the stopping points along theelectron trajectory. Even though the sub-cycling curveshows noticeable decrease in the length of the spikes, itstill exhibits a significant departure from R th = 1, whichis similar to that of the curve without the sub-cycling.In this section, we will show that the just mentionedaccumulation of the dephasing error is associated withtemporal interpolation of temporally discretized electricand magnetic fields. As explained in Sec. III, a commonapproach is to use linear or first order temporal interpo-lation for the magnetic field. In the case of sub-cycling,we use the same interpolation procedure for the electricfield. For the rest of the paper, this approach will bereferred to as the 1 st order interpolation.An example of the error produced by the 1 st orderinterpolation in the perfectly propagating PIC algorithmis shown in the lower panel of Fig. 4. B a /B is shown inthe upper panel, where B a is the analytical solution forthe magnetic field of a propagating laser while B > | ( B a − B t ) /B | , where B t is the temporally interpolatedfield using the analytical field B a known at a half time-step before and after on a fixed spatial grid. The biggesterrors occur close to the maxima and minima of B a wherethe field has a quadratic rather than linear dependenceon the phase variable ξ and thus time t .Though the errors shown in Fig. 4 are seemingly small,they can have a profound impact on the energy gain byaltering the electron dynamics. We illustrate this aspectby considering the temporal field discretization used bythe perfectly propagating PIC without the sub-cycling. FIG. 4: Errors induced by 1 st to 5 th order temporalinterpolation in the perfectly propagating PIC. Upper panel:magnetic field B a of the pulse. Lower panel: interpolationerror, where B t is interpolated to t = n ∆ t using B a at t = ( n + 1 / t and at preceding time-steps. The temporaldiscretization is performed using ∆ t = λ/ c . This algorithm uses an analytically propagated field thatrequires a 1 st order temporal interpolation of the mag-netic field to advance the electron momentum. Analyt-ically, as given by Eq. (2), the magnetic field acting onthe electron is equal to the electric field. This is becausethe laser propagation is in a vacuum, which results in thephase velocity, v ph , being equal to the speed of light, i.e. v ph = c . Problematically, the 1 st order temporal inter-polation effectively reduces the magnetic field amplitude.As a result, this interpolation makes the field configura-tion similar to that of a laser pulse with a superluminalphase velocity: B = Ec/v ph . (29)It can be shown that the integral of motion in such apulse [12] is γ − v ph c p x mc = const (30)rather than γ − p x /mc = const. Eq. (30) is a manifesta-tion of the changes in the electron dynamics due to theinterpolation.We now use the analogy between the interpolated fieldand that of a superluminal laser pulse to estimate theinduced error in the dephasing R . For simplicity, we ne-glect spatial discretization in the derivation that follows.The extrapolated field is B t = B a − ∆ B, (31)so that, after taking into account that E = B a , we obtain B t /E = 1 − ∆ B/B a . (32) It then follows from Eq. (29) that the expression on theright-hand side can be interpreted as c/v ph , such that theeffective phase velocity is given by v ph /c ≈ B/B a . (33)Note that the second term on the right-hand side isnever negative during the 1 st order interpolation, sothat v ph /c ≥
1. We substitute the expression given byEq. (33) into Eq. (30) to find that γ − p x mc = const + ∆ BB a p x mc . (34)Recall that we are using the analytical form for thefields, so the dephasing defined by Eq. (9) is equal to R = γ − p x /mc . It then follows from Eq. (34) that thereduction in the magnetic field strength due to the inter-polation, ∆ B/B a >
0, increases the dephasing betweenthe electron and the wave, R = γ − p x mc = const + ∆ BB a p x mc . (35)An important consequence of this result is that evensmall interpolation errors can lead to significant errorsin the dephasing for an electron with an ultra-relativisticlongitudinal momentum, i.e. p x (cid:29) mc .In order to validate the mechanism linking the tem-poral interpolation with the dephasing errors, we re-examine the plot from Fig. 3a for the algorithm withoutthe sub-cycling. Fig. 5 zooms in around a stopping pointat tω ≈ R − R th ) /R th , whereas Fig. 5e shows the error thatwe expect due to the temporal interpolation according toEq. (35). In order to be consistent with the general PICapproach, we apply the triangular particle shape to calcu-late both B a and ∆ B . At tω < R − R th ) /R th remarkablywell, which confirms our estimates and the assertion thatmost of the errors at p x /mc (cid:29) R − R th ) /R th at the stopping point. To aid thecomparison, we only use the sub-cycling for the consid-ered stopping point, i.e. no sub-cycling is applied at tω < st ordertemporal interpolation.It must be pointed out that the sub-cycling algorithminherits the field interpolation error. Temporal interpo-lation during the sub-cycling causes ∆ B/B a to rapidlyoscillate (see Fig. 5b). To understand the nature of theoscillations, recall that the magnetic field is known at n − / n + 1 /
2, which means that there is no mag-netic field interpolation error at those time indices. Theresult is that ∆ B during the sub-cycling oscillates be-tween 0 and some interpolation errors. FIG. 5: Performance of perfectly propagating PIC using the1 st order interpolation with and without the sub-cyclingnear a stopping point. (a): longitudinal electron momentum.(b): the error between the interpolated magnetic field andthe analytical magnetic field. (c) and (d): relative error inthe dephasing rate. (e): prediction for the relative error inthe dephasing rate given by Eq. (35). The sub-cycling(Ψ max = 0 .
01) is only performed near the stopping pointshown in these plots. Both simulations use a = 50,∆ t = λ/ c and ∆ x ≈ λ/ . We further investigate the impact of the sub-cycling byapplying the sub-cycling algorithm with Ψ max = 0 . max = 0 .
01. However, the number oftime integration steps increases by a factor of 4 to roughly2 . × steps. In terms of the number of sub-cyclingsteps, the increase is by a factor of 21 from 72,000 stepsto 1 . × steps. Even though the sub-cycling doessolve some local problems at the stopping points, it isunable to prevent accumulation of ( R − R th ) /R th evenwhen using a computationally expensive Ψ max in the sub-cycling algorithm. This accumulation leads to significantdeviations from the analytical solution (for example, seeFigs. 1a and 1b for Ψ max = 0 . st order temporal field interpolation has theeffect that is similar to that of a superluminal laser pulse,resulting siginificant errors for ultra-relativistic particles.The effective superluminosity is particularly evident fromthe momentum plot in Fig. 1 that resembles that for anelectron in a plane wave with v ph > c [13]. VI. HIGHER-ORDER FIELD INTERPOLATION
In order to reduce the errors caused by the 1 st order(linear) temporal interpolation of the fields, we changethe interpolation procedure by increasing the order of thepolynomial using Lagrange interpolation. As discussedlater in this section, we find that a 2 nd order interpo-lation produces only marginal improvements whereas a3 rd order interpolation dramatically reduces numericalerrors.The 1 st order interpolation for the magnetic field with-out the sub-cycling has the form B z | n = 12 (cid:0) B z | n +1 / + B z | n − / (cid:1) . (36)However, the information about the magnetic field isknown at t < ( n − / t , so we can leverage it to moreaccurately describe B z | n : B z | n = p (cid:88) j =0 B z | n +1 / j − p L p,j , (37)where p is the order of the interpolation and L p,j are theLagrange coefficients defined as L p,j = p (cid:89) k =0 ,k (cid:54) = j p − k − / j − k . (38)For example, the 2 nd order interpolation is given by B z | n = 18 (cid:0) − B z | n − / + 6 B z | n − / + 3 B z | n +1 / (cid:1) . (39)In order to apply the interpolation during the sub-cycling, we generalize Eqs. (37) and (38) to the case whenthe field is interpolated to t = ∆ t ( n − /
2) + ∆ t (cid:48) . Thenthe corresponding field that we denote as B (cid:48) z is given by B (cid:48) z = p (cid:88) j =0 B z | n +1 / j − p L Bp,j , (40)where L Bp,j = p (cid:89) k =0 ,k (cid:54) = j p − k − t (cid:48) / ∆ tj − k . (41)These expressions are used for ∆ t ∗ ≤ ∆ t (cid:48) ≤ ∆ t , where∆ t ∗ is the sub-cycling step. As previously stated, thesub-cycling algorithm also requires the electric field tobe temporally interpolated. The interpolated expressionat t = ∆ t ( n − /
2) + ∆ t (cid:48) is given by E (cid:48) y = p (cid:88) j =0 E y | n + j − p L Ep,j , (42)where L Ep,j = p (cid:89) k =0 ,k (cid:54) = j p − k − / t (cid:48) / ∆ tj − k . (43)We first implement the 2 nd order temporal field inter-polation to the perfectly propagating PIC algorithm asdescribed by Eq. (39). The dephasing errors for simu-lations with and without the sub-cycling are shown inFig. 3b. The simulation setup and parameters are thesame as those used in Fig. 3a, which uses the 1 st orderinterpolation instead. The use of the 2 nd order interpo-lation together with the sub-cycling significantly reducesthe length of the downward spikes at the stopping points.However, there is still a large error accumulation over theduration of the entire simulation. This behavior is quali-tatively different from what we observe for the 1 st orderinterpolation where R/R th − nd order interpo-lation that uses two data points from the “past” ( n − / n − /
2) and only one data point from the “future”( n + 1 / n + 1 / n + 3 /
2) and only one data point from the“past” ( n − / R/R th − R/R th − rd order temporal interpolationof the fields is sufficient to achieve a dramatic improve-ment. As shown in Fig. 3c, the dephasing error is reduceddramatically when the 3 rd order interpolation is used to-gether with the sub-cycling. Increasing the order of theinterpolation does not resolve the asymmetry issue, but itdoes reduce the magnitude of the errors (see lower panelof Fig. 4).The role of the higher-order temporal interpolation isto reduce the temporal interpolation errors. The curvefor the simulation without the sub-cycling in Fig. 3cshows that the 3 rd order interpolation makes the erroraccumulations mostly happen near the stopping points.This means that the interpolation errors are no longerthe main source of the dephasing error. This conclusionis reaffirmed as we see that the errors near the stoppingpoints are mitigated when adding the sub-cycling algo-rithm to the 3 rd order interpolation. The resulting R inFig. 3c is very close to R th along the entire trajectory,with very slight downward ticks at the stopping points.Figures 1a and 1b confirm that the PIC algorithm thatuses both the sub-cycling and the 3 rd order interpolationis able to recover the electron trajectory and the momen-tum evolution.We conclude this section by performing a broad laseramplitude scan by increasing a from 5 to 100. In Fig. 2b,the green markers show (cid:15) energy for a perfectly propagat-ing PIC that uses both the sub-cycling and the 3 rd orderinterpolation. The errors are dramatically reduced com-pared to the algorithm that uses the 3 rd order interpo-lation without the sub-cycling. Fig. 2a confirms that asimilar trend is in place for the standard PIC algorithmwhere the fields are numerically propagated. These scans FIG. 6: Relative dephasing error, (cid:15) dephasing , for differentversions of the PIC algorithm over a range of laseramplitudes a . All simulation results are obtained using thestandard PIC algorithm. Each data set representssimulations using a different combination of temporal fieldinterpolation orders, with/without sub-cycling, anddiscretization resolution. All simulations use c ∆ t/ ∆ x = 0 . max = 0 .
01 is with sub-cycling. show that a higher order interpolation is necessary for thesub-cycling algorithm to yield significant benefits.
VII. IMPROVEMENTS DUE TO rd AND th ORDERTEMPORAL INTERPOLATION
In this section, we quantify improvements in precisiondelivered by implementing a higher-order temporal fieldinterpolation and the sub-cycling into a standard PICalgorithm. We compare the results for 3 rd and 5 th orderinterpolations to those obtained using the standard 1 st order interpolation.Fig. 6 shows (cid:15) dephasing as a function of a for simula-tions using the 5 th order temporal interpolation with andwithout the sub-cycling. The 3 rd order points in thisfigure are exactly the same points as found in Fig. 2a.Without the sub-cycling, there is essentially no changein (cid:15) dephasing when increasing the interpolation order from3 rd to 5 th . When the algorithm employs the sub-cycling,0 FIG. 7: Improvements due to the 3 rd order interpolation, (cid:15) dephasing /(cid:15) dephasing , and required increase in the number oftime steps, δ steps , for PIC algorithm with sub-cycling(Ψ max = 0 . c ∆ t/ ∆ x = 0 . the 5 th order interpolation does reduce the errors com-pared to the 3 rd order interpolation.Continuing in Fig. 6, it is instructive to determine theresolution needed in the standard PIC algorithm (1 st or-der interpolation and no sub-cycling) to achieve low val-ues of (cid:15) dephasing . We use the results of the new algo-rithm (3 rd order interpolation with sub-cycling, shownby the green points) for comparison. Scanning over5 ≤ a ≤
100 and using a resolution of ∆ t = λ/ c ,the new algorithm maintained (cid:15) dephasing < . t = λ/ c is required to achieve similar results (shownby the yellow points). We kept c ∆ t/ ∆ x = 0 .
99 duringboth of these scans. Therefore, not only the number oftime-steps has increased by a factor of 8, but the numberof cells has also increased by the same factor. This com-parison shows that improving the precision of the stan-dard PIC algorithm through a straightforward reductionof ∆ t and ∆ x can become prohibitively expensive.To have a more direct comparison for how importantthe high order temporal field interpolation is in ensuringeffective sub-cycling, we introduce χ = (cid:15) dephasing (cid:14) (cid:15) dephasing (44) χ = (cid:15) dephasing (cid:14) (cid:15) dephasing , (45)where the superscripts indicate the order of the temporalinterpolation in a simulation with the sub-cycling. Fig-ures 7 and 8 show χ and χ as functions of a for threedifferent values of λ/ (∆ tc ) which denotes the resolutionof both (cid:15) dephasing and (cid:15) dephasing in the χ calculation. FIG. 8: Improvements due to the 5 th order interpolation, (cid:15) dephasing /(cid:15) dephasing , and required increase in the number oftime steps, δ steps , for PIC algorithm with sub-cycling(Ψ max = 0 . c ∆ t/ ∆ x = 0 . A significant feature of Fig. 7a is the maximum of χ at a ≈
25 for ∆ t = λ/ c . The location of the maximumshifts to higher values of a as we increase the resolution.The implication of this trend is that the resolution mustbe adjusted according to the value of a in order to main-tain the effectiveness of the higher order interpolation.We find that χ is reduced with the increase of a pri-marily during the ultra-relativistic electron motion. Thiswas confirmed by first observing that an almost identi-cal trend exists for the perfectly propagating PIC. Wethen replaced the fields in the perfectly propagating PICby their analytical values without the interpolation alongthe ultra-relativistic parts of the trajectory, which caused χ to become monotonically increasing. We note that c − v x decreases with the increase of a , where v x is thelongitudinal electron velocity along the ultra-relativisticpart of the trajectory. The electron motion becomes moresensitive to errors in the effective phase velocity causedinterpolation as we increase a (see Eq. 35).An alternative to increasing the resolution is to in-crease the order of the interpolation. Fig. 8a shows thatthe interpolation errors are reduced sufficiently to pre-vent a rollover of χ at a ≤
100 for ∆ t = λ/ c . Thepenalty for using the 5 th order interpolation is the needto store more information about the electric and mag-netic fields. However, one can use relatively large cellsand time-steps, the saving of which might be significantwhen using this approach for 3D simulations.When the resolution λ/ (∆ tc ) becomes more fine in theimproved PIC algorithm, (cid:15) dephasing unsurprisingly de-1creases. Though we see in Fig. 8a that χ is reduced withlarger λ/ (∆ tc ). This shows that simulations with coarserresolutions have the most potential to be improved byincreasing the order of the temporal field interpolationgiven that the interpolation is sufficiently high-order.Figures 7 and 8 also show the relative increase in thenumber of time-steps when using the sub-cycling. Wedefine δ steps = N ∗ N ∗ + N , (46)where N is the number of ∆ t sized time-steps taken dur-ing the simulation. ∆ t is the fixed time-step size of theYee scheme and the without sub-cycling Boris algorithm. N ∗ is the number of sub-cycling time-steps (∆ t ∗ < ∆ t )taken during the simulation. As expected, by using thesame sub-cycling condition, δ steps is mostly unaffectedbecause the interpolation error χ is small in the con-sidered examples and the electron trajectories are verysimilar for the 3 rd and 5 th order interpolations.As a increases, δ steps decreases. This happens be-cause the travel time between the stopping points in-creases with a . The sub-cycling is applied only in thevicinity of the stopping points, which means that the sub-cycling has to be applied less frequently at higher a . Atfiner resolutions, δ steps also decreases as Ψ is scales pro-portionally with the time-step size. VIII. SUMMARY AND DISCUSSION
In the 1D PIC tests where the errors induced by thefield solver are negligible, the standard linear interpo-lation of the fields temporally can introduce significanterrors to the motion of charged particles. These errorscan be interpreted as an effective increase in the phasevelocity of the laser fields, and are strongly pronouncedin the ultra-relativistic regime. We have shown that ahigh order temporal interpolation for electric and mag-netic fields is required in addition to particle sub-cyclingnear stopping points to accurately reproduce electron dy-namics in ultra-high intensity laser pulses. The need forthe high-order interpolation is a general requirement forsimulations at a (cid:29)
1, so our findings equally apply toimplementations that employ a different particle pusherwith an improved performance at a (cid:29)
1, such as theone detailed in Ref. [9].We have found that a third rather than second ordertemporal interpolation is needed in order to fully leveragethe benefits of the sub-cycling algorithm and thus dra-matically improve the simulation accuracy. While thehigher order temporal field interpolations would requireincreases in memory usage, the sub-cycling allows forfairly coarse discretizations to achieve results comparableto standard simulations using very fine discretizations.In the considered problem, the accuracy of the solutionis primarily impacted by the temporal field interpolationand the errors introduced by the particle pusher, while the errors introduced by the field solver are relativelyinconsequential. If, in a more general case, the errorsfrom the field solver become significant, then it may bejustifiable to consider using a higher order solver. Sincethe errors introduced by the field solver are separate fromthe time-interpolation errors, the two issues should betreated separately.
ACKNOWLEDGEMENTS
The work of K.T. and A.A. was supported by the Na-tional Science Foundation (PHY 1821944). G. C.’s workwas supported by the Exascale Computing Project (grantno. 17-SC-20-SC), a collaborative effort of the U.S. De-partment of Energy Office of Science and the NationalNuclear Security Administration. C. N. Danson, C. Haefner, J. Bromage, T. Butcher, J.-C. F.Chanteloup, E. A. Chowdhury, A. Galvanauskas, L. A. Gizzi,J. Hein, D. I. Hillier, and et al., High Power Laser Science andEngineering , e54 (2019). F. Lureau, G. Matras, O. Chalus, C. Derycke, T. Morbieu,C. Radier, O. Casagrande, S. Laux, S. Ricaud, G. Rey, andet al., High Power Laser Science and Engineering , e43 (2020). C. Birdsall and A. Langdon,
Plasma physics via computer simu-lation , The Adam Hilger series on plasma physics (McGraw-Hill,1985). K. Yee, IEEE Transactions on antennas and propagation , 302(1966). A. Arefiev, G. Cochran, D. Schumacher, A. Robinson,and G. Chen, Physics of Plasmas , 013103 (2015),https://doi.org/10.1063/1.4905523. D. F. Gordon, B. Hafizi, and J. Palastro,AIP Conference Proceedings , 050002 (2017),https://aip.scitation.org/doi/pdf/10.1063/1.4975863. A. P. L. Robinson, K. Tangtartharakul, K. Weichman,and A. V. Arefiev, Physics of Plasmas , 093110 (2019),https://doi.org/10.1063/1.5115993. A. Sorokovikova, A. V. Arefiev, C. McGuffey, B. Qiao, A. P. L.Robinson, M. S. Wei, H. S. McLean, and F. N. Beg, Phys. Rev.Lett. , 155001 (2016). D. Gordon and B. Hafizi, Computer Physics Communications , 107628 (2021). J. V. Shebalin, IEEE transactions on plasma science , 390(1988). A. V. Arefiev, V. N. Khudik, and M. Schollmeier, Physics ofPlasmas , 033104 (2014), https://doi.org/10.1063/1.4867491. V. Khudik, A. Arefiev, X. Zhang, and G. Shvets, Physics ofPlasmas , 103108 (2016), https://doi.org/10.1063/1.4964901. A. Robinson, A. Arefiev, and V. Khudik, Physics of Plasmas ,083114 (2015), https://doi.org/10.1063/1.4928893. APPENDIXA. Electron Dynamics in a Plane Wave
The equations of motion for an electron have the form: d p dt = −| e | E − | e | γmc [ p × B ] , (47) d r dt = cγ p mc . (48)2We are considering a plane electromagnetic wave prop-agating along the x -axis in the positive direction. Thefields can then be described using a normalized vectorpotential a that is only a function of a phase variable ξ ,with E y = B z = − mωc | e | dadξ , (49) ξ ≡ ω ( t − x/c ) . (50)The three components of Eq. (47) now read dp x dt = − | e | γmc p y B z , (51) dp y dt = −| e | E y + | e | γmc p x B z , (52) dp z dt = 0 . (53)It follows from Eq. (53) that p z is conserved. In this workwe consider an electron that is initially at rest, so that p z = 0 (54)during the electron motion in the laser pulse. Next, wesubstitute the expressions for the fields given by Eq. (49)into Eq. (52) and take into account that v = p /γm tofind that dp y dt = mωc (cid:16) − v x c (cid:17) dadξ . (55)This equation can be simplified even further, because dξdt = ∂ξ∂t + dxdt ∂ξ∂x = ω (cid:16) − v x c (cid:17) . (56)We use this relation in Eq. (55) to obtain an equationthat determines the evolution of p y for a given vectorpotential a : ddt (cid:16) p y mc − a (cid:17) = 0 . (57)We then have p y /mc = a ( ξ ) (58)for an electron that is immobile prior to the arrival of thelaser pulse.Equations (51) - (53) have another integral of motion, ddt (cid:16) γ − p x mc (cid:17) = 0 , (59)where γ = (cid:112) p /m c . In order to show this, we addEq. (51) multiplied by p x to Eq. (52) multiplied by p y ,which yields 12 dp dt = −| e | E y p y . (60) Taking into account the definition for γ , we find that dγdt = − | e | E y p y γm c . (61)On the other hand, it follows from Eq. (51) that dp x dt = − | e | γmc p y E y , (62)where we used Eq. (49) to replace B z with E y . It cannow be directly verified that Eq. (59) holds.Using the same initial conditions as before (immobileelectron), we find from Eq. (59) that γ − p x mc = 1 . (63)This relation reduces to an equation for p x after we takeinto account that p y = amc and p z = 0: γ = (cid:18) p x m c + a (cid:19) / . (64)The solution of the resulting equation is p x /mc = a / . (65)Equations (65), (58), and (54) fully describe the evolutionof the momentum for a given normalized vector potential a = a ( ξ ). It also follows from Eq. (63) that γ = 1 + a / . (66) B. Energy Conservation Analysis
It directly follows from the equations of motion for anelectron [see Eq. (47)] that dε k dt = −| e | ( E · v ) , (67)where ε k ≡ ( γ − mc (68)is the kinetic energy, v = p /γm is the electron velocity,and γ = (1+ p /m e c ) / is the relativistic factor. In ourtest problem, the electron is initially at rest, with ε k = 0at t = 0. Then the kinetic energy at time t is given by ε k ( t ) = −| e | (cid:90) t ( E · v ) dt (cid:48) . (69)The particle pusher that we use in this manuscript doesnot automatically guarantee that Eq. (69) is satisfied. Itis thus insightful to examine the discrepancy between thework that is done by the laser electric field (right-handside) and the change in the kinetic energy (left-hand side)3in our simulations. In order to aid the comparison, wedefine the kinetic energy in our algorithm as ε n +1 / k = mc ( γ n +1 / − . (70)The work, w , done on the electron by the electric field asthe velocity changes from v n − / to v n +1 / is w n +1 / = q ∆ t E n · ( v n − / + v n +1 / ) . (71)The net work done from the beginning of the simulationsis then ε n +1 / w ≡ n (cid:88) i =0 w i +1 / . (72) FIG. 9: Energy conservation in a PIC algorithm with 3 rd order temporal interpolation without the sub-cyclingalgorithm. (top panel) The kinetic energy, ε k , and the network done, ε w , normalized by ε max . (bottom panel) Thedifference between the work done and the change in kineticenergy over one ∆ t time-step, ∆ ε k − ∆ ε w , normalized by∆ ε t . The simulation parameters are set to a = 100, λ/ (∆ tc ) = 75, and c ∆ t/ ∆ x = 0 . In Fig. 9, the top panel shows ε k and ε w in a PIC algo-rithm that employs the 3 rd order temporal field interpola-tion without the sub-cycling. The simulation parametersare provided in the caption. The curves are normalizedby ε max = mc a /
2, which is the theoretical maximumkinetic energy of an electron with the considered initialconditions (we neglected the p y contribution as p x (cid:29) p y for a (cid:29) ε k /ε max and ε w /ε max remain signif-icantly below unity along the electron trajectory, whichindicates that the algorithm fails to correctly reproduce the energy gain. On the other hand, there is no visible FIG. 10: Energy conservation in a PIC algorithm with 3 rd order temporal interpolation with the sub-cycling algorithm.(top panel) The kinetic energy, ε k , and the net work done, ε w , normalized by ε max . (bottom panel) The differencebetween the work done and the change in kinetic energyover one ∆ t time-step, ∆ ε k − ∆ ε w , normalized by ∆ ε t . Thesimulation parameters are set to a = 100, λ/ (∆ tc ) = 75, c ∆ t/ ∆ x = 0 .
99, and Ψ max = 0 . difference between the two curves. We thus conclude thatthe lack of energy conservation of the Boris pusher is notthe primary factor causing the shortfall in the maximumenergy gain.The bottom panel in Fig. 9 shows the difference be-tween the change in kinetic energy over one time-step,∆ ε k , and and the work done over one time-step, ∆ ε w ,in the same simulation. The difference is normalized by∆ ε t = | e | cE ∆ t , which is an upper limit for the workdone over time interval ∆ tt