Particle-in-cell codes are now standard tools for studying ultra-intense laser-plasma interactions. Motivated by direct laser acceleration of electrons in sub-critical plasmas, we examine temporal resolution requirements that must be satisfied to accurately calculate electron dynamics in strong laser fields. Using the motion of a single electron in a perfect plane electromagnetic wave as a test problem, we show surprising deterioration of the numerical accuracy with increasing wave amplitude *a*_{0} for a given time-step. We go on to show analytically that the time-step must be significantly less than *λ*/*ca*_{0} to achieve good accuracy. We thus propose adaptive electron sub-cycling as an efficient remedy.

## I. INTRODUCTION

Ongoing progress in laser engineering has significantly increased the maximum laser intensity available for ultra-intense laser-plasma experiments. A number of emerging applications now rely on the ability of high intensity laser beams to accelerate electrons to relativistic energies that considerably exceed the electron rest mass.^{1–3} It has therefore become critical to accurately simulate electron dynamics in an ultra-intense electromagnetic field. However, it is not always possible to determine whether the simulation has been done with required accuracy and whether the simulation results are physically correct without performing a convergence study. Such a study can be a time-consuming effort without guarantee of a conclusive outcome if the parameter space has not been narrowed down sufficiently using relevant test cases. Convergence studies are best performed based on a fundamental understanding of the numerical and physical constraints on simulation parameters.

A commonly used tool for simulating laser-plasma interactions is a particle-in-cell (PIC) code that consists of two key blocks: a wave solver on a given spatial grid and a particle pusher that uses the calculated fields to advance particles. Two key parameters, besides the number of particles, which determine the accuracy and speed of the simulation, are the grid-size (cell-size) and the time-step. The choice of the two is interrelated through the Courant criterion,^{4} which limits the maximum time-step for a given grid-size, particularly for explicit codes. Depending on the problem, different criteria are used to determine the grid-size and the time-step.

We are particularly interested in the regime where the plasma density is significantly below the critical density^{6–8} or, equivalently, where the laser frequency is significantly above the plasma frequency. In this regime, the group and phase velocities of an electromagnetic laser pulse inside the plasma are close to the speed of light, which enables electron acceleration to high energies as the electron moves forward with the pulse via direct laser acceleration. For this type of problem, one typically determines the spatial resolution first based on the wavelength of the laser pulse and its transverse dimensions. It is however not necessary to resolve the Debye length when simulating electron acceleration by an ultra-intense laser in a significantly under-dense plasma, since the energy that results from numerical heating is inconsequential (see Sec. VII for a detailed explanation). The time-step is then determined using the Courant criterion, so that the numerical scheme remains stable. It is usually chosen close to its maximum value allowed by the Courant criterion in order to reduce the numerical dispersion of the wave caused by the grid. However, we find that at sufficiently high intensities in a low density environment, the wavelength does not set the scale for accurate treatment of electron motion. Instead, the electron motion near its stopping points becomes the critical factor resulting in surprisingly stringent requirements for convergence. Given the current widespread use of a large number of differing PIC codes, often incorporating multiple algorithms, simple test cases and criteria for evaluating them are highly beneficial. In particular, since experiments employing intensities of up to 10^{21 }W/cm^{2} (Ref. 9) are now commonplace and experiments at significantly higher intensities are underway or anticipated,^{10} evaluation of PIC codes in this regime is crucial.

We revisit the dynamics of a single free electron irradiated by a high-intensity plane electromagnetic wave as a test problem for evaluating the performance of a particle-in-cell code. The electron motion becomes relativistic at large normalized wave amplitudes, $a0\u2261|e|E/mec\omega \u226b1$, where *E* is the amplitude of the wave electric field, *ω* is the wave frequency, *c* is the speed of light, and *e* and *m _{e}* are the electron charge and mass. In this regime, most of the electron energy is associated with the longitudinal motion and the maximum relativistic

*γ*-factor increases as $\gamma \u2248a02/2$. As a result of the longitudinal motion with relativistic velocity, the frequency of the transverse electron oscillations decreases by a factor of

*γ*, which can become substantial for large wave amplitudes. If, for a given time-step, the simulation reproduces the electron motion with

*a*

_{0}∼ 1 correctly, one might expect that the relativistic motion of an electron with

*a*

_{0}≫ 1 would also be correctly reproduced since the period of the oscillations increases with

*a*

_{0}. The numerical results presented below (1D-3V) show an exactly opposite trend, with the accuracy quickly deteriorating with the increase of

*a*

_{0}for a fixed time-step.

In what follows, we outline a criterion that must be considered when simulating electron acceleration by a high amplitude electromagnetic wave in an under-dense plasma. We show that the electron dynamics can be correctly reproduced only if the time-step is sufficiently small to resolve the electron motion near stopping points along the trajectory. This condition requires that the time-step in the simulation is less than 1/*a*_{0}*ω*, where *ω* is the wave frequency. This criterion is independent of constraints on spatial resolution. It becomes more stringent at higher wave amplitudes due to the fact that the acceleration is more rapid near the stopping points for larger *a*_{0}. This means that the error accumulates primarily along relatively small segments of the electron trajectory in the vicinity of the stopping points. We therefore propose adaptive electron sub-cycling as an efficient remedy. The idea is to reduce the time-step for a given electron when the acceleration can no longer be correctly reproduced using the original time-step. Our results show that sub-cycling permits a dramatic increase in accuracy with only a modest increase in the total number of time steps. Given current interest in direct laser acceleration,^{11–13} the rapidly increasing focus on using ultra-intense lasers to study radiation reaction^{14,15} and the ambition to explore QED effects with lasers in the near future,^{16} it is important that the new constraint on time-step described here be taken into account.

The rest of the paper is organized as follows. In Sec. II, we review the dynamics of a free electron irradiated by an incoming electromagnetic wave to establish the context for the analysis that follows. In Sec. III, we demonstrate using a significantly under-dense plasma that relative numerical errors grow with wave amplitude in the case of an electron accelerated by a laser pulse. In Secs. IV and V, we analyze the errors originating from the particle pusher and derive a corresponding criterion for the time-step. In Sec. VI, we show that the errors can be greatly reduced using adaptive sub-cycling which helps to better resolve the electron dynamics near stopping points. Finally, in Sec. VII, we summarize our results and discuss possible implementation of the sub-cycling in a particle-in-cell code.

## II. SINGLE ELECTRON DYNAMICS IN A PLANE WAVE

In this section, we summarize the key features of single electron dynamics in a plane wave in order to establish the context for the subsequent analysis of the numerical results. We consider a free electron irradiated by a plane wave that propagates along the *x*-axis. The wave electric field is directed along the *y*-axis and the wave magnetic field is directed along the *z*-axis. The wave propagation can be described using a normalized vector potential

where *a* is only a function of a dimensionless phase variable

Here, *λ* is the wave-length, *c* is the speed of light, *t* is the time in the laboratory frame of reference, and **e**_{y} is a unit vector. The electric and magnetic fields of the wave are given by

where *m _{e}* and

*e*are the electron mass and charge.

An initially stationary electron irradiated by this wave moves only in the (*x*, *y*)-plane according to the following equations:

where *p _{x}* and

*p*are components of the electron momentum and

_{y}is the relativistic factor. This system of equations has two integrals of motion

The second integral of motion implies that the electron dephases from the wave at a constant rate. In order to show that, we first take the derivative of the phase variable *ξ* defined by Eq. (2) with respect to time *t*, which yields

On the other hand, the proper time that we denote as *τ* and the time *t* are related by the expression *dτ*/*dt* = 1/*γ*. The proper time is the elapsed time that would be measured by the electron itself. Using this relation in Eq. (12), we find that

According to the integral of motion (11), the expression on the right-hand side is a constant and, therefore, $d\xi /d\tau $ is also a constant. This means that the phase of the field sampled by the electron increases linearly with proper time. It is then appropriate to interpret the integral of motion (11) as the corresponding dephasing rate. It should be emphasized that the value of the dephasing rate has a direct and significant impact on the maximum energy that the electron gains during acceleration by the wave.^{6,8} This aspect will play a key role in the subsequent analysis in Secs. III–VI.

If the electron is at rest (*p _{x}* =

*p*= 0) before the wave arrives (

_{y}*a*= 0), then it follows from Eq. (11) that

In momentum space, the electron always moves along a parabola $px/mec=(py/mec)2/2$, with only the maximum displacement along *p _{x}* and

*p*changing with the wave amplitude.

_{y}Equations (15) and (16) give *p _{x}* and

*p*only in terms of

_{y}*ξ*and one still needs to integrate Eq. (7) in order to find

*p*(

_{x}*x*,

*t*) and

*p*(

_{y}*x*,

*t*). One can find from Eq. (7) using the definition of

*ξ*given by Eq. (2) that

where $\gamma =1+a2(\xi )/2$ according to Eqs. (9), (15), and (16). Equations (15)–(18) allow one to implicitly determine components of the electron momentum as functions of time *t* and axial distance *x*.

To summarize, a free electron irradiated by a plane wave moves along a parabola in momentum space due to the fact that its dephasing rate $\gamma \u2212px/mec$ remains constant. The maximum electron *γ*-factor is $\gamma *=1+a02/2$ and the corresponding maximum electron energy gain is $\gamma *mec2$, where *a*_{0} is the maximum value of *a*.

## III. SIMULATION RESULTS FOR A SIGNIFICANTLY UNDER-DENSE PLASMA

In this section, we present several results from 1D-3V particle-in-cell simulations (one spatial and three velocity dimensions) in order to determine how well the single electron dynamics described in Sec. II is reproduced numerically. We initialize a low-density hydrogenic plasma slab ($ne=ni=10\u22123\u2009ncrit$, where *n*_{crit} is the critical density) with cold electrons that are irradiated by a plane electromagnetic wave. The plasma density is deliberately set very low so that space-charge effects are negligible during our runs and each electron effectively behaves as a free electron in a vacuum. The plasma is essentially acting as a convenient cold electron source. At this density, the effect of the electron currents on the field of the wave is also negligible in our runs. This setup mimics the initial conditions considered in Sec. II. The normalized pulse amplitude *a* ramps up to *a*_{0} and then remains constant. The exact profile of the electron density and the wave amplitude are not important. All the plasma electrons are equivalent in this setup and they only differ by their initial location. In what follows, we select one electron and track it throughout the simulation.

Figure 1 shows numerical results for *a*_{0} = 10 and three different grid-sizes, Δ*x*/*λ* = 1/20, 1/40, and 1/80. In all three runs, the ratio of the time-step Δ*t* to the grid-size Δ*x* was the same, with *cΔt*/Δ*x* = 0.95. The upper panel in Fig. 1 shows the trajectory of a single electron in momentum space and the lower panel shows the corresponding dephasing rate as a function of the distance traveled by the electron from its initial location. The electron data is shown with dots, because it was recorded at discrete time intervals $dt=0.075\lambda /c$. According to the analytical solution of Sec. II, the electron should be moving along a parabola $px/mec=(py/mec)2/2$, with $|py/mec|\u2264a0$ and $0\u2264px/mec\u2264a02/2$. The dephasing rate must remain constant and equal to unity, $\gamma \u2212px/mec=1$. Not surprisingly, the convergence to the analytical solution improves as we decrease Δ*x*, which in this case is equivalent to decreasing Δ*t* because their ratio is maintained.

Figure 1 shows that there is a correlation between the deviation of the dephasing rate from unity and the deviation of the numerical solution in momentum space from the analytical result. This trend is also not surprising, because the dephasing rate determines how the field acting on the electron changes in time. An error in the dephasing causes an error in the electron acceleration and, consequently, leads to an error in the electron momentum. Note that errors in the dephasing rate also lead to considerable asymmetry in the transverse momentum *p _{y}* (see the results for Δ

*x*=

*λ*/20 in Figs. 1 and 2), which can result in an unphysical electron drift perpendicular to the direction of the wave propagation.

There is a distinct periodic structure of sharp downwards spikes in the dephasing rate. By comparing the time evolution of the dephasing rate and the electron momentum, we find that the downwards spikes correspond to stopping points. The numerical errors in the dephasing rate change considerably along the electron trajectory and they are the most significant in the vicinities of the stopping points. This is immediately evident from the change in the vertical distance between adjacent data points (dots) in the lower panel of Fig. 1. As stated earlier, the time interval between the adjacent data points is constant and, therefore, the dephasing rate changes at a greater rate around the downwards spikes.

Figure 2 shows the same quantities as Fig. 1, but for a higher wave amplitude of *a*_{0} = 15. The already discussed trends seem to be similar in this case. We deliberately used the same set of grid-sizes and the same ratio *cΔt*/Δ*x* as in the case of *a*_{0} = 10 in order to determine how numerical errors scale with wave amplitude. Let us examine the runs with Δ*x* = *λ*/40 for *a*_{0} = 10 and *a*_{0} = 15. The period of electron oscillations in the laser field increases with *a*_{0} due to the increased *γ*-factor that is primarily associated with the longitudinal motion. This can be seen by comparing the distance between the downwards spikes in the dephasing rate that is also the distance between stopping points. Since the velocity of the longitudinal motion is close to *c*, longer distance translates directly into a longer interval between stopping points and thus a longer period of oscillations. The time-step for both runs at *a*_{0} = 10 and *a*_{0} = 15 is the same, so one might expect that the higher amplitude run would be better resolved and the numerical errors would be reduced. The comparison of the deviation in the dephasing rate from unity indicates that the trend is exactly the opposite. Greater errors in the dephasing rate then lead to greater errors in the electron momentum at higher wave amplitude.

## IV. ANALYSIS OF ERRORS ORIGINATING FROM THE PARTICLE-PUSHER

In general, there is a wide range of factors that contribute to numerical errors in a particle-in-cell simulation (see Refs. 17, 18, 19, 4, and references therein). The results of Sec. III indicate that in our test problem the numerical errors in the dephasing rate tend to significantly increase around specific points of the electron trajectory—the stopping points. This observation serves as a motivation for us to consider the particle pusher separately.

In what follows, we specify the wave field analytically. In the test problem under consideration, this is simply an electromagnetic pulse propagating in a vacuum. This approach allows us to isolate the errors introduced by the particle pusher. We use the standard Boris pusher^{4} to advance the electron momentum and coordinates in time using the analytical solution for the field at the electron location. The following results are for an initially immobile single electron irradiated by a plane electromagnetic wave.

Figure 3 shows the results for three different wave amplitudes and the same time-step, *cΔt*/*λ* = 1/50. The upper panel shows the electron momentum space. Deviation from the analytical solution predicting $px/mec=(py/mec)2/2$ [see Eqs. (15) and (16)] increases with wave amplitude. The lower panel shows the dephasing rate as a function of the longitudinal displacement. The errors in the dephasing rate and the deviation from $\gamma \u2212px/mec=1$ also increase with *a*_{0}. Both trends for the momentum and dephasing are similar to those that emerge from comparing the results for *a*_{0} = 10 and *a*_{0} = 15 in Sec. III (see Figs. 1 and 2). However, the wave fields in Sec. III were calculated numerically using a finite difference scheme for the Maxwell equations. This suggests that the increase in numerical errors is caused by the particle pusher, whereas the errors resulting from numerical integration of the field equations are less significant for the grid-size and time-step used.

In order to see the trend for the numerical errors more clearly, we have performed field amplitude scans for three different time-steps, *cΔt*/*λ* = 1/60, 1/120, and 1/180. The results are shown in Fig. 4 for 5 ≤ *a*_{0} ≤ 40. The upper panel shows the discrepancy in the electron energy gain and the lower panel shows the discrepancy in the dephasing rate. For each field amplitude, we determine the maximum *γ*-factor achieved by the electron in the simulation, denoted as *γ*_{max}. The analytical solution in Sec. II predicts the maximum *γ*-factor to be

where *a*_{0} is the peak amplitude of the pulse. The quantity plotted in the upper panel of Fig. 4 is $(\gamma max\u2212\gamma *)/\gamma *$, which is the relative error in the predicted electron energy gain. For each wave amplitude, we also determine the maximum absolute deviation of the dephasing rate $\gamma \u2212px/mec$ from unity. This is the quantity plotted in the lower panel of Fig. 4.

There are several important trends that become apparent from the scans presented in Fig. 4. The error in the dephasing rate gradually increases with *a*_{0} regardless of the time-step used to integrate the electron equations of motion. In contrast with that, the discrepancy in the electron energy gain exhibits a threshold behavior with the increase of *a*_{0}. The threshold for the onset of large errors in the electron energy gain appears to scale inversely proportional to the time-step used in the particle pusher. If, for example, we define this threshold as $|\gamma max\u2212\gamma *|/\gamma *=0.025$, then we find that the discrepancy exceeds the threshold value at *a*_{0} = 11, 25, and 36 for *cΔt*/*λ* = 1/60, 1/120, and 1/180. This is indeed a threshold, since the discrepancy in the energy gain jumps up considerably at these wave amplitudes. We have listed these numbers in Table I. We have also listed in Table I $|\gamma \u2212px/mec\u22121|max$ above the threshold. The threshold clearly occurs at roughly the same level of discrepancy in the dephasing rate in all three cases. The values from Table I are shown with enlarged markers in Fig. 4.

cΔt/λ
. | a_{0}
. | $|\gamma \u2212px/mec\u22121|max$ . |
---|---|---|

1/60 | 11 | 0.080 |

1/120 | 25 | 0.088 |

1/180 | 36 | 0.083 |

cΔt/λ
. | a_{0}
. | $|\gamma \u2212px/mec\u22121|max$ . |
---|---|---|

1/60 | 11 | 0.080 |

1/120 | 25 | 0.088 |

1/180 | 36 | 0.083 |

We can summarize this section by concluding that there is a general trend for the errors in the dephasing rate to increase with wave amplitude. Once the discrepancy in the dephasing rate approaches 10%, a rapid increase in the discrepancy in the electron energy gain takes place. The corresponding threshold wave amplitude scales inversely proportional to the time-step used to integrate the electron equations of motion.

## V. CRITERION FOR THE TIME-STEP

In Secs. III and IV, we showed that the numerical errors tend to increase with wave amplitude. In this section, we formulate a criterion for the time-step that must be satisfied in order to accurately reproduce the electron dynamics in a strong electromagnetic wave.

It is helpful to begin by reviewing how a standard Boris particle pusher advances the electron momentum in time. This is done in three subsequent stages using given electric and magnetic fields. It first accelerates the electron using the electric field. It then performs a rotation in the magnetic field. The last stage is another push using the electric field. The (half) rotation angle is

where *B* is a given magnetic field and *γ* is the relativistic factor of the electron after the first push by the electric field.

This procedure necessarily requires the rotation angle to be small, which can be easily understood in the context of the test problem that we are considering. As the electron is pushed forward by the wave, its transverse momentum oscillates, while the longitudinal momentum remains positive. Therefore, the electron never performs a full rotation in momentum space. On the other hand, if the rotation angle in the Boris pusher is comparable to *π*, then the particle pusher qualitatively changes the electron motion causing the electron to move backwards.

Our next step is to determine the relation between the rotation angle and the wave amplitude in our test problem. Let us take a pulse that, after some initial ramp-up, has a constant amplitude, with $a=a0\u2009sin(\xi )$. In this case, the electric and magnetic fields acting on the electron are

According to the analytical solution of Sec. II, the *γ*-factor of the electron is

The ratio *B*/*γ* has the largest absolute value for *a* = 0. Thus, the rotation angle for a given time-step Δ*t* has also the largest value for *a* = 0

The requirement $max|\psi |\u226a\pi $ now yields

This condition indicates that the wave amplitude imposes an upper limit on the time-step that can be used by the particle pusher. It should be noted that the discussed criterion is equivalent to the requirement that the smallest time scale in the problem must be resolved. In our case, the restriction is imposed by the wave magnetic field and the corresponding time scale that must be resolved is the gyro-period at *a* = 0.

The points along the electron trajectory where *a* vanishes are stopping points. This is evident from the analytical solution of Sec. II. Equations (15) and (16) give *p _{y}* and

*p*as functions of

_{z}*a*, with both vanishing for

*a*= 0. On the other hand, the electric and magnetic fields of the wave reach their maximum amplitude at a stopping point. Therefore, the vicinity of a stopping point is that part of the electron trajectory where the electron experiences the strongest acceleration and, as a result, the rotation angle has the largest amplitude.

It might seem that the derived criterion for the time-step is unnecessarily restrictive due to the fact that it is imposed by electron dynamics near stopping points. Indeed, even if the criterion is not satisfied, the electron momentum gain that would be calculated incorrectly would still be relatively small compared to both $|py|=a0$ and $px=a02/2$. However, this argument does not take into account the corresponding error in the dephasing rate and its impact on the subsequent electron acceleration by the wave.

In order to show that errors in acceleration near stopping points are indeed critical, let us first estimate how the dephasing rate changes as a result of these errors. We define the dephasing rate as

Let us consider an electron as it starts its motion right at the stopping point with correct initial conditions determined from the analytical solution, so that *p _{x}* =

*p*= 0 and

_{y}*a*= 0. The electron dephasing rate is then

*I*= 1. We now use a standard Boris pusher to advance the electron momentum by one time-step. We assume that the rotation is calculated incorrectly and we want to estimate the resulting error in the dephasing rate. The momentum gain resulting from the acceleration by the wave electric field over a time interval equal to Δ

*t*is roughly $|\Delta p/mec|=2\pi a0c\Delta t/\lambda $. Deliberately assuming the worst case scenario, we set

*p*= 0 and $px=\xb1|\Delta p|$ in Eq. (25). For $|\Delta p/mec|\u226a1$, we have the following estimate for the dephasing after just one time-step:

_{y}An electron with an initial axial momentum $px=\xb1|\Delta p|$ would have the same dephasing rate. An analytical solution for such an electron is given in Ref. 6. The maximum electron *γ*-factor according to Eq. (25) of Ref. 6 is $\gamma max=(1+a02+I2)/2I$. The difference in the maximum *γ*-factor achieved by this electron as compared to the case considered in Sec. II where the electron is initially at rest is

We can now employ this expression to estimate the impact that an error in the dephasing has on the electron energy gain. We use the estimate for the dephasing rate given by Eq. (26) and assume that the error in the dephasing rate is relatively small, $|1\u2212I|\u226a1$. It then follows from Eq. (27) that

This estimate indicates that the relative error in the electron energy gain is of the same order as the error in dephasing rate.

These estimates elucidate the physical basis for the restriction on the time-step given by Eq. (24). If the condition (24) is not satisfied, then the numerical error resulting from numerical integration of the electron equations of motion near a stopping point leads to a considerable error in the dephasing rate. The error in the momentum gain is relatively small at this stage. However, the error in the dephasing rate affects the subsequent electron acceleration by the wave even if no additional errors are introduced, causing a considerable error in the maximum electron energy gain.

## VI. ADAPTIVE SUB-CYCLING

We have so far examined the errors introduced by a standard Boris pusher and determined that the electron dynamics in the vicinity of stopping points imposes a stringent upper limit on the time-step. Guided by these observations, we develop in this section a procedure that allows us to considerably improve the accuracy when calculating the electron dynamics with only a modest increase of the total number of time steps.

A direct way to reduce numerical errors generated by the particle pusher is by decreasing the time-step Δ*t*. However, if the time-step reduction necessary to satisfy the criterion (24) is significant, then this approach would greatly increase the total number of time-steps needed to simulate the same time interval. On the other hand, the increased precision is not helpful and thus unnecessary for the majority of the electron trajectory.

A more efficient way to reduce the numerical errors is by adaptively decreasing the time-step only when necessary to ensure that the criterion (24) is always satisfied. We implement this by introducing a critical rotation angle $\psi *$. We also choose a base value for the time-step Δ*t* = Δ_{0}. At the beginning of each time-step, we estimate the rotation angle for the particle pusher, *ψ _{est}*, using Eq. (20), where the value of

*γ*is taken from the previous time-step and Δ

*t*= Δ

_{0}. If this estimate exceeds $\psi *$, then we reduce Δ

*t*, which reduces

*ψ*, to make

_{est}*ψ*less than $\psi *$. For convenience, we choose the reduced time-step from a list of discrete values $\Delta t=\Delta 0/4k$, where

_{est}*k*= 1,2,3,…. We pick the largest value that yields $\psi est<\psi *$. The algorithm automatically resets Δ

*t*to the base value Δ

_{0}if a reduced time-step is no longer necessary.

The standard Boris pusher updates electron position and momentum at interleaved time points, staggered such that they “leapfrog” over each other. In order to successfully implement the described algorithm, one must synchronize the electron position and momentum at a time point corresponding to the momentum before changing the time-step. Therefore, the particle pushing algorithm should involve the following two steps after the momentum has been updated and the time-step has to be changed. The first step is to advance the electron position by Δ*t*/2, where Δ*t* is the original time-step. The second step is to advance the electron position by Δ*t*/2, where Δ*t* is the new time-step. After this, the standard particle pushing algorithm can be used with a new time-step, starting with a momentum update.

### A. Sub-cycling example

As an example, we re-run the three cases presented in Fig. 3 (see Sec. IV) now using adaptive sub-cycling. The base time-step for all three runs is set to $c\Delta 0/\lambda =1/50$, which was the time-step used to generate Fig. 3. The critical rotation angle is set to $\psi *=0.05$. The resulting electron momentum and electron dephasing are shown in Fig. 5. The deviation of the dephasing rate from unity has been significantly reduced, as compared to the results in the lower panel in Fig. 3. It is less than 10% even for *a*_{0} = 25. As a result, the deviation of the electron momentum from the analytical solution given in Sec. IV is no longer visually detectable. It has been dramatically reduced compared to the results in the upper panel in Fig. 3.

In Table II, we have listed how many times each time-step value is used during the sub-cycling. The relative number of steps using the base time-step value increases with increasing *a*_{0}, so that the number of reduced time-steps is less than 13% for *a*_{0} = 25. Recall that the same value of *a*_{0} had the largest deviation from the analytical solution in Fig. 3. Therefore, the sub-cycling algorithm becomes more efficient at higher wave amplitudes and a significant improvement of the numerical results can be achieved using only a modest increase in the total number of time-steps. This is not surprising, since the electron spends only a small fraction of its time near the stopping points, as the time interval between the stopping points increases as $\gamma max\u221da02$.

a_{0}
. | $\psi *$ . | Δt (%)
. | Δt/4 (%)
. | Δt/16 (%)
. | Δt/64 (%)
. |
---|---|---|---|---|---|

5 | 0.05 | 58 | 30 | 12 | 0 |

15 | 0.05 | 78 | 13 | 8 | 1 |

25 | 0.05 | 87 | 7 | 4 | 2 |

a_{0}
. | $\psi *$ . | Δt (%)
. | Δt/4 (%)
. | Δt/16 (%)
. | Δt/64 (%)
. |
---|---|---|---|---|---|

5 | 0.05 | 58 | 30 | 12 | 0 |

15 | 0.05 | 78 | 13 | 8 | 1 |

25 | 0.05 | 87 | 7 | 4 | 2 |

We have performed wave amplitude scans for three different base time-step values, $c\Delta 0/\lambda =1/60$, 1/120, and 1/180 for $5\u2264a0\u226440$. The critical rotation angle was set to $\psi *=0.05$. Without sub-cycling, these three scans would produce the results shown in Fig. 4. The results using sub-cycling are shown in Fig. 6. As in the case of Fig. 4, the upper panel shows the discrepancy in the electron energy gain and the lower panel shows the discrepancy in the dephasing rate. For each wave amplitude, we determine the maximum *γ*-factor achieved by the electron in the simulation, denoted as $\gamma max$. The analytical solution in Sec. II predicts the maximum *γ*-factor to be $\gamma *=1+a02/2$, where *a*_{0} is the peak amplitude of the pulse. The quantity plotted in the upper panel of Fig. 6 is $(\gamma max\u2212\gamma *)/\gamma *$, which is the relative error in the predicted electron energy gain. For each wave amplitude, we also determine the maximum absolute deviation of the dephasing rate $\gamma \u2212px/mec$ from unity. This is the quantity plotted in the lower panel of Fig. 6.

The trends seen in Fig. 5 are confirmed by the wave amplitude scans in Fig. 6. The deviation from the analytical solution of the dephasing and the electron energy gain has been dramatically reduced. The errors in the dephasing rate remain below 10% for all runs. The errors in the electron energy gain now increase gradually with *a*_{0} and do not exhibit the threshold behavior seen in Fig. 4. This is due to the fact that the amplitude of the rotation angle in the particle pusher always remains smaller than $\psi *=0.05$, which prevents the particle pusher from introducing qualitative changes to the electron dynamics near the turning points.

### B. Efficiency of the sub-cycling algorithm

The key metric of the proposed algorithm is the relative increase in the total number of time steps between two stopping points due to the sub-cycling. In order to evaluate this quantity, we consider a vicinity of a stopping point *ξ* = 0 where the wave amplitude is $a(\xi )=a0\u2009sin(\xi )$. We assume that *a*_{0} ≫ 1 and that the condition given by Eq. (24) is not satisfied for Δ*t* = Δ _{0}, so that the sub-cycling is required.

The fact that the criterion (24) is not satisfied implies that the rotation angle *ψ* defined by Eq. (20) exceeds $\psi *$ at *ξ* = 0. As the electron accelerates from the stopping point, the rotation angle decreases with the increase of *ξ*. We use Eq. (20) together with the expressions for the fields given by Eq. (21) to find that

The condition $\psi (0)=\psi *$ sets the smallest time-step $\Delta tmin$ that the sub-cycling procedure would have to use, which is

The sub-cycling is only needed for $|\xi |<\xi *$, where $\xi *$ is defined by the condition $\psi (\xi *)=\psi *$. Away from the stopping point, where $\psi <\psi *$, the sub-cycling is not necessary. We assume that *ψ* matches $\psi *$ close to *ξ* = 0, with $\xi \u226a\pi $. Otherwise, the base time-step is too large and most of the electron trajectory has to be sub-cycled. The reduction of *ψ* at *ξ* ≪ 1 occurs due to the denominator in Eq. (29) when the term proportional to $sin(\xi )$ starts to dominate. In order to find the value of *ξ* that yields $\psi =\psi *$, we retain only the $sin(\xi )$-term in the denominator and expand the expression in Eq. (29) with respect to *ξ*. To the lowest order in *ξ*, we then have

It follows from Eq. (31) that $\psi =\psi *$ at

The travel time from the stopping point to the point where $\xi =\xi *$ can be calculated using Eq. (17). Taking into account that $\gamma =1+a2(\xi )/2$ and assuming that *ξ* ≪ 1, we find that to the lowest order in $\xi *$ this travel time is

On the other hand, the travel time to the next stopping point is

We found this expression by performing the integration in Eq. (17) from *ξ* = 0 to $\xi =\pi $ and retaining only the leading term, which involves *a*_{0}.

The total number of time-steps with Δ*t* = Δ_{0} between two stopping points is

The total number of time-steps $N*$ with the sub-cycling can be estimated by using $\Delta t=\Delta tmin$ for the time interval that requires the sub-cycling

This gives an upper estimate, since the actual algorithm is adaptive and not all time-steps are as small as $\Delta tmin$. The relative increase in the total number of time-steps is

We next use the expressions for $t*$, *t _{s}*, $\Delta tmin$, and $\xi *$ to obtain that

Therefore, the relative number of extra time-steps decreases at least as fast as $1/a0$, which makes the sub-cycling mechanism more efficient at higher wave amplitudes.

Finally, it is important to point out the dual role of the adjustable parameter $\psi *$. On the one hand, the number of extra steps increases as $1/\psi *5/2$ according to Eq. (38), which means that the efficiency of the sub-cycling decreases with a decrease of the critical rotation angle $\psi *$. On the other hand, the error in the dephasing also decreases with a decrease of $\psi *$. This can be illustrated using the results presented in Fig. 5. In addition to the oscillations in the dephasing between the stopping points, there is also a clear downwards shift introduced by changing the size of the time-step during the sub-cycling procedure. This trend leads to an error of roughly 8% at the end of the run for *a*_{0} = 25. By reducing $\psi *$ from 0.05 to 0.025 and holding other parameters fixed, we have reduced the dephasing error to less than 2%. Therefore, the value of the critical rotation angle $\psi *$ must be carefully chosen such that the desired precision is achieved without sacrificing the efficiency.

## VII. SUMMARY AND DISCUSSION

We have revisited the classic test problem of a single electron irradiated by a high-intensity plane electromagnetic wave to examine the performance of a particle-in-cell code. We found that the numerical accuracy consistently deteriorates with increasing wave amplitude for a fixed time-step. We have separately examined the accuracy of the standard Boris particle pusher and found that the particle pusher introduces significant errors while integrating the electron motion in the vicinities of stopping points. Field amplitude scans reveal that the deviation of the dephasing rate from unity increases with *a*_{0} regardless of the time-step. The errors in the energy gain have a threshold behavior with increasing *a*_{0} caused by increased error in the dephasing rate and the threshold amplitude scales inversely proportional to Δ*t*.

Based on these observations, we have derived a convenient criterion for the time-step used in the particle pusher

Our analysis shows that the stopping points are more prone to numerical errors than other parts of the electron trajectory. Numerical errors from the integration of the electron equation of motion near a stopping point lead to considerable error in the dephasing rate. The error in the momentum gain is relatively small at this stage; however, the error in the dephasing rate affects the subsequent electron acceleration by the wave causing considerable error in the maximum electron energy gain.

We have developed an efficient algorithm that allows us to reduce the numerical errors generated by the particle pusher. The algorithm uses the derived criterion to adaptively decrease the time-step to ensure that the criterion is always satisfied. We have demonstrated that this adaptive sub-cycling actually becomes more efficient at higher wave amplitudes, so that a significant improvement of the numerical results can be achieved using only a modest increase in the total number of time-steps.

In our analysis, we have used a standard Boris pusher to calculate the electron dynamics, since this is the integrator implemented in some of the PIC codes frequently used to model laser-plasma interactions. It has been recently pointed out that the Boris pusher might lead to errors when calculating orbits of relativistic electrons.^{18} This can be important on those segments of the electron trajectory where the electron is moving forward with a high *γ*-factor and the electric and magnetic field contributions in the equation of motion almost cancel each other out. The Boris pusher does not correctly cancel the electric and magnetic field contributions in the case of a relativistic electron, which can lead to a spurious force.^{18} The spurious force and the resulting errors can be eliminated by using the Vay particle pusher^{18} instead of the Boris pusher. In the regime considered in this paper, the relative role of these errors is however not as significant as the role of the errors from integrating the electron equation of motion near the stopping points. This is evident from the dramatic improvement that we have achieved by implementing the sub-cycling algorithm while still using the Boris pusher.

Our discussion has so far been focused on the particle pusher as the main source of numerical errors. The degradation in numerical accuracy observed in PIC simulations under these conditions appears to be primarily from the particle pusher, rather than the field advance. However, another potential source of errors is the numerical dispersion produced by the field solver. It is well known that the frequency *ω* of a wave with wave-length *λ* is less than $\omega *\u22612\pi c/\lambda $

This error can be greatly reduced by keeping the ratio $c\Delta t/\Delta x$ close to unity. The corresponding discussion can be found in Ref. 4.

The spatial grid can also lead to numerical electron heating. The heating is caused by nonphysical instabilities that develop if the grid size exceeds the Debye length (for example, see p. 175 in Ref. 4). The resulting electron energy *ε* can be estimated by noting that the heating stops once the Debye length $\lambda D\u2248\epsilon /4\pi ne2$ becomes comparable to Δ*x*, where *n* is the plasma electron density. We are interested in a regime where *n* is significantly under-dense for a given frequency *ω* of the laser pulse irradiating the plasma. It is then convenient to introduce a critical density *n*_{cr} defined by the condition $4\pi ncre2/me=\omega $. The wave propagation in a significantly under-dense plasma is similar to that in a vacuum and thus we effectively have $\omega \u22482\pi c/\lambda $. Taking this relation and the definition for *n*_{cr} into account, we find that the condition $\lambda D\u2248\Delta x$ yields

Typically, the grid-size is much less than $\lambda /2\pi $, so that the electron energy resulting from numerical heating is smaller than $mec2n/ncr$. This energy is non-relativistic in a significantly under-dense plasma with $n\u226ancr$. On the other hand, the maximum energy of an electron accelerated by a wave with amplitude *a*_{0} is $mec2a02/2$. This energy greatly exceeds *ε*, provided that $a0\u226b1$ and $n\u226ancr$. Moreover, the dephasing rate $I=\gamma \u2212px/mec$ [see Eq. (25)] remains close to unity despite the numerical heating since $\epsilon \u226amec2$. This means that the electron acceleration by the wave remains unaffected. Therefore, it is not necessary to resolve the Debye length when simulating electron acceleration by an ultra-intense laser $(a0\u226b1)$ in a significantly under-dense plasma, since the energy that results from numerical heating is inconsequential.

The errors from the particle pusher that we have examined here should be particularly critical when simulating electron acceleration in an underdense plasma. In a two- and three-dimensional set-up, only a small group of electrons are accelerated directly by the laser pulse. Our approach to adaptive sub-cycling would be well suited in this case to improve the accuracy of the numerical results. Our algorithm would automatically single out energetic electrons accelerated by a high-amplitude field, leaving the time-step for the other electrons unchanged. In order to implement the adaptive sub-cycling in a given particle-in-cell code, one has to orbit-average the current density of the sub-cycled electrons in a way similar to those discussed in Refs. 20 and 21. The orbit-averaged current can then be directly used in the field solver to calculate the fields using the global time-step.

Finally, it is worth emphasizing that the primary role of the proposed algorithm is to ensure that the multi-scale electron dynamics that takes place in the presence of a wave with a large amplitude *a*_{0} is recovered in a simulation by resolving the fundamental time-scale of the electron motion given by Eq. (39). If for given simulation parameters the criterion (39) is not satisfied, then significant errors near stopping points are unavoidable regardless of the numerical scheme. These errors can lead to significant errors in electron energy gain even when the total energy in the system may appear well conserved due to the use of an energy-conserving scheme. In fact, we have also observed the effect reported in Sec. III using the energy-conserving particle advance in the PIC code LSP.^{22} This is the reason why adaptive particle orbit integration was introduced in an implicit energy-conserving scheme described in Ref. 23. These general considerations also indicate that the proposed algorithm should be applied regardless of the plasma density whenever the time-scale set by a conventional PIC code is not sufficient to resolve the time-scale given by Eq. (39). The impact of the discussed errors in the case of a dense plasma still remains to be evaluated and only then can it be determined whether they are significant enough or can be ignored.

## ACKNOWLEDGMENTS

A.V.A. would like to thank Dr. S. P. D. Mangles and Dr. V. N. Khudik for stimulating discussions and constructive comments. Simulations for this paper were performed using the EPOCH code (developed under UK EPSRC Grant Nos. EP/G054940/1, EP/G055165/1, and EP/G056803/1) using HPC resources provided by the Texas Advanced Computing Center at The University of Texas. A.V.A. was supported by AFOSR Contract No. FA9550-14-1-0045, National Nuclear Security Administration Contract No. DE-FC52-08NA28512, and U.S. Department of Energy Contract No. DE-FG02-04ER54742. GEC received support from NNSA Contract No. DE-NA0001976.