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 a0 for a given time-step. We go on to show analytically that the time-step must be significantly less than λ/ca0 to achieve good accuracy. We thus propose adaptive electron sub-cycling as an efficient remedy.

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 density6–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 1021 W/cm2 (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|e|E/mecω1, where E is the amplitude of the wave electric field, ω is the wave frequency, c is the speed of light, and e and me 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 γa02/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 a0 ∼ 1 correctly, one might expect that the relativistic motion of an electron with a0 ≫ 1 would also be correctly reproduced since the period of the oscillations increases with a0. The numerical results presented below (1D-3V) show an exactly opposite trend, with the accuracy quickly deteriorating with the increase of a0 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/a0ω, 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 a0. 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 reaction14,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.

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

a(x,t)=a(ξ)ey,
(1)

where a is only a function of a dimensionless phase variable

ξ2π(ctx)/λ.
(2)

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

E=mec|e|at,
(3)
B=mec2|e|ax,
(4)

where me 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:

ddt(pxmec)=|e|Bγmecpymec,
(5)
ddt(pymec)=|e|Emec+|e|Bγmecpxmec,
(6)
dxdt=cγpxmec,
(7)
dydt=cγpymec,
(8)

where px and py are components of the electron momentum and

γ=1+(px/mec)2+(py/mec)2,
(9)

is the relativistic factor. This system of equations has two integrals of motion

ddt(pymeca)=0,
(10)
ddt(γpxmec)=0.
(11)

We skip the derivation here, which can, for example, be found in Refs. 5 and 6.

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

dξdt=ωγ(γpxmec).
(12)

On the other hand, the proper time that we denote as τ and the time t are related by the expression /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

dξdτ=ω(γpxmec).
(13)

According to the integral of motion (11), the expression on the right-hand side is a constant and, therefore, dξ/dτ 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. IIIVI.

If the electron is at rest (px = py = 0) before the wave arrives (a = 0), then it follows from Eq. (11) that

γpxmec=1.
(14)

Using Eqs. (10) and (14), and the definition of γ, we find that

py/mec=a,
(15)
px/mec=a2/2.
(16)

In momentum space, the electron always moves along a parabola px/mec=(py/mec)2/2, with only the maximum displacement along px and py changing with the wave amplitude.

Equations (15) and (16) give px and py only in terms of ξ and one still needs to integrate Eq. (7) in order to find px(x, t) and py(x, t). One can find from Eq. (7) using the definition of ξ given by Eq. (2) that

ctλ=12π0ξγdξ,
(17)
xλ=12π0ξ(γ1)dξ,
(18)

where γ=1+a2(ξ)/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 γpx/mec remains constant. The maximum electron γ-factor is γ*=1+a02/2 and the corresponding maximum electron energy gain is γ*mec2, where a0 is the maximum value of a.

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=103ncrit, where ncrit 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 a0 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 a0 = 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Δtx = 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λ/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|a0 and 0px/meca02/2. The dephasing rate must remain constant and equal to unity, γpx/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.

FIG. 1.

Electron momentum space (top) and the dephasing rate (bottom) for different values of grid-size, Δx, at a0 = 10. For all three runs, we set cΔt/Δx=0.95.

FIG. 1.

Electron momentum space (top) and the dephasing rate (bottom) for different values of grid-size, Δx, at a0 = 10. For all three runs, we set cΔt/Δx=0.95.

Close modal

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 py (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.

FIG. 2.

Electron momentum space (top) and the dephasing rate (bottom) for different values of grid-size, Δx, at a0 = 15. For all three runs, we set cΔt/Δx=0.95.

FIG. 2.

Electron momentum space (top) and the dephasing rate (bottom) for different values of grid-size, Δx, at a0 = 15. For all three runs, we set cΔt/Δx=0.95.

Close modal

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 a0 = 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Δtx as in the case of a0 = 10 in order to determine how numerical errors scale with wave amplitude. Let us examine the runs with Δx = λ/40 for a0 = 10 and a0 = 15. The period of electron oscillations in the laser field increases with a0 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 a0 = 10 and a0 = 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.

We therefore conclude that the numerical accuracy deteriorates with increasing wave amplitude for a fixed time-step despite the fact that the period of the electron oscillations increases. In Secs. IVVII, we examine the source of the discovered increasing inaccuracy.

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 pusher4 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 γpx/mec=1 also increase with a0. Both trends for the momentum and dephasing are similar to those that emerge from comparing the results for a0 = 10 and a0 = 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.

FIG. 3.

Electron momentum space (top) and the dephasing rate (bottom) for a0 = 5, 15, and 25 calculated using a particle pusher and an analytical field. We set cΔt/λ = 1/50 for all three runs.

FIG. 3.

Electron momentum space (top) and the dephasing rate (bottom) for a0 = 5, 15, and 25 calculated using a particle pusher and an analytical field. We set cΔt/λ = 1/50 for all three runs.

Close modal

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 ≤ a0 ≤ 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

γ*=1+a02/2,
(19)

where a0 is the peak amplitude of the pulse. The quantity plotted in the upper panel of Fig. 4 is (γmaxγ*)/γ*, 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 γpx/mec from unity. This is the quantity plotted in the lower panel of Fig. 4.

FIG. 4.

Field amplitude scans of the relative error in the electron energy gain (top) and the error in the dephasing rate (bottom) using three different time steps: cΔt/λ = 1/60, 1/120, and 1/180. See discussion in text of these figures of merit. The enlarged markers indicate the threshold intensity for the onset of large errors in the electron energy gain for each value of cΔt/λ (see Table I).

FIG. 4.

Field amplitude scans of the relative error in the electron energy gain (top) and the error in the dephasing rate (bottom) using three different time steps: cΔt/λ = 1/60, 1/120, and 1/180. See discussion in text of these figures of merit. The enlarged markers indicate the threshold intensity for the onset of large errors in the electron energy gain for each value of cΔt/λ (see Table I).

Close modal

There are several important trends that become apparent from the scans presented in Fig. 4. The error in the dephasing rate gradually increases with a0 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 a0. 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 |γmaxγ*|/γ*=0.025, then we find that the discrepancy exceeds the threshold value at a0 = 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|γpx/mec1|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.

TABLE I.

Threshold for errors in the energy gain.

cΔta0|γpx/mec1|max
1/60 11 0.080 
1/120 25 0.088 
1/180 36 0.083 
cΔta0|γpx/mec1|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.

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

ψ=12|e|Bλγmec2cΔtλ,
(20)

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=a0sin(ξ). In this case, the electric and magnetic fields acting on the electron are

E=B=mec2a0cos(ξ)2π/λ|e|.
(21)

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

γ=1+a2/2.
(22)

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

max|ψ|=πa0cΔtλ.
(23)

The requirement max|ψ|π now yields

cΔtλ1a0.
(24)

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 py and pz as functions of 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

Iγpxmec=1+(pxmec)2+(pymec)2pxmec.
(25)

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 px = py = 0 and 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 |Δp/mec|=2πa0cΔt/λ. Deliberately assuming the worst case scenario, we set py = 0 and px=±|Δp| in Eq. (25). For |Δp/mec|1, we have the following estimate for the dephasing after just one time-step:

I1±|Δp/mec|1±2πa0cΔt/λ.
(26)

An electron with an initial axial momentum px=±|Δ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 γ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

γmaxγ*γ*=(1I)1+a02II(2+a02).
(27)

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, |1I|1. It then follows from Eq. (27) that

γmaxγ*γ*1I±2πa0cΔtλ.
(28)

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.

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 ψ*. 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 ψ*, then we reduce Δt, which reduces ψest, to make ψest less than ψ*. For convenience, we choose the reduced time-step from a list of discrete values Δt=Δ0/4k, where k = 1,2,3,…. We pick the largest value that yields ψest<ψ*. 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.

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Δ0/λ=1/50, which was the time-step used to generate Fig. 3. The critical rotation angle is set to ψ*=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 a0 = 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.

FIG. 5.

Electron momentum space (top) and the dephasing rate (bottom) for a0 = 5, 15, and 25 calculated using adaptive sub-cycling. The particle pusher uses an analytical field. The base time-step is set at cΔ0/λ = 1/50 for all three runs.

FIG. 5.

Electron momentum space (top) and the dephasing rate (bottom) for a0 = 5, 15, and 25 calculated using adaptive sub-cycling. The particle pusher uses an analytical field. The base time-step is set at cΔ0/λ = 1/50 for all three runs.

Close modal

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 a0, so that the number of reduced time-steps is less than 13% for a0 = 25. Recall that the same value of a0 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 γmaxa02.

TABLE II.

Number of time-steps used for sub-cycling.

a0ψ*Δt (%)Δt/4 (%)Δt/16 (%)Δt/64 (%)
0.05 58 30 12 
15 0.05 78 13 
25 0.05 87 
a0ψ*Δt (%)Δt/4 (%)Δt/16 (%)Δt/64 (%)
0.05 58 30 12 
15 0.05 78 13 
25 0.05 87 

We have performed wave amplitude scans for three different base time-step values, cΔ0/λ=1/60, 1/120, and 1/180 for 5a040. The critical rotation angle was set to ψ*=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 γmax. The analytical solution in Sec. II predicts the maximum γ-factor to be γ*=1+a02/2, where a0 is the peak amplitude of the pulse. The quantity plotted in the upper panel of Fig. 6 is (γmaxγ*)/γ*, 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 γpx/mec from unity. This is the quantity plotted in the lower panel of Fig. 6.

FIG. 6.

Field amplitude scans of the relative error in the electron energy gain (top) and the error in the dephasing rate (bottom) using adaptive sub-cycling for base time-step values cΔ0/λ = 1/60, 1/120, and 1/180.

FIG. 6.

Field amplitude scans of the relative error in the electron energy gain (top) and the error in the dephasing rate (bottom) using adaptive sub-cycling for base time-step values cΔ0/λ = 1/60, 1/120, and 1/180.

Close modal

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 a0 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 ψ*=0.05, which prevents the particle pusher from introducing qualitative changes to the electron dynamics near the turning points.

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(ξ)=a0sin(ξ). We assume that a0 ≫ 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 ψ* 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

ψ(ξ)=πcΔtλa0cos(ξ)1+a02sin2(ξ)/2.
(29)

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

cΔtminλ=ψ*πa0.
(30)

The sub-cycling is only needed for |ξ|<ξ*, where ξ* is defined by the condition ψ(ξ*)=ψ*. Away from the stopping point, where ψ<ψ*, the sub-cycling is not necessary. We assume that ψ matches ψ* close to ξ = 0, with ξπ. 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(ξ) starts to dominate. In order to find the value of ξ that yields ψ=ψ*, we retain only the sin(ξ)-term in the denominator and expand the expression in Eq. (29) with respect to ξ. To the lowest order in ξ, we then have

ψ(ξ)=πcΔ0λ2a0ξ2.
(31)

It follows from Eq. (31) that ψ=ψ* at

ξξ*1a02πψ*cΔ0λ.
(32)

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

ctλ=ct*λ12πa022ξ*33.
(33)

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

ctλ=ctsλa028.
(34)

We found this expression by performing the integration in Eq. (17) from ξ = 0 to ξ=π and retaining only the leading term, which involves a0.

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

N0=ts/Δ0.
(35)

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

N*=ts2t*Δ0+2t*Δtmin.
(36)

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

N*N0N02t*tsΔ0Δtmin.
(37)

We next use the expressions for t*, ts, Δtmin, and ξ* to obtain that

N*N0N023πa02ξ*5=23π1a0(2πψ*cΔ0λ)5/2.
(38)

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 ψ*. On the one hand, the number of extra steps increases as 1/ψ*5/2 according to Eq. (38), which means that the efficiency of the sub-cycling decreases with a decrease of the critical rotation angle ψ*. On the other hand, the error in the dephasing also decreases with a decrease of ψ*. 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 a0 = 25. By reducing ψ* 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 ψ* must be carefully chosen such that the desired precision is achieved without sacrificing the efficiency.

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 a0 regardless of the time-step. The errors in the energy gain have a threshold behavior with increasing a0 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

cΔt/λ1/a0.
(39)

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 pusher18 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 ω*2πc/λ

ωω*ω*π26(Δxλ)2[ 1(cΔtΔx)2 ].
(40)

This error can be greatly reduced by keeping the ratio cΔt/Δ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 λDε/4π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 ncr defined by the condition 4πncre2/me=ω. The wave propagation in a significantly under-dense plasma is similar to that in a vacuum and thus we effectively have ω2πc/λ. Taking this relation and the definition for ncr into account, we find that the condition λDΔx yields

εmec2nncr(2πΔxλ)2.
(41)

Typically, the grid-size is much less than λ/2π, 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 nncr. On the other hand, the maximum energy of an electron accelerated by a wave with amplitude a0 is mec2a02/2. This energy greatly exceeds ε, provided that a01 and nncr. Moreover, the dephasing rate I=γpx/mec [see Eq. (25)] remains close to unity despite the numerical heating since εmec2. 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 (a01) 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 a0 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.

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.

1.
S. P. D.
Mangles
,
C. D.
Murphy
,
Z.
Najmudin
,
A. G. R.
Thomas
,
J. L.
Collier
,
A. E.
Dangor
,
E. J.
Divall
,
P. S.
Foster
,
J. G.
Gallacher
,
C. J.
Hooker
,
D. A.
Jaroszynski
,
A. J.
Langley
,
W. B.
Mori
,
P. A.
Norreys
,
F. S.
Tsung
,
R.
Viskup
,
B. R.
Walton
, and
K.
Krushelnick
,
Nature
431
,
535
(
2004
).
2.
C. G. R.
Geddes
,
Cs.
Toth
,
J.
van Tilborg
,
E.
Esarey
,
C. B.
Schroeder
,
D.
Bruhwiler
,
C.
Nieter
,
J.
Cary
, and
W. P.
Leemans
,
Nature
431
,
538
(
2004
).
3.
J.
Faure
,
Y.
Glinec
,
A.
Pukhov
,
S.
Kiselev
,
S.
Gordienko
,
E.
Lefebvre
,
J.-P.
Rousseau
,
F.
Burgy
, and
V.
Malka
,
Nature
431
,
541
(
2004
).
4.
C. K.
Birdsall
and
A. B.
Langdon
,
Plasma Physics Via Computer Simulation
(
IOP Publishing Ltd.
,
2002
).
5.
T. J. M.
Boyd
and
J. J.
Sanderson
,
The Physics of Plasmas
(
Cambridge University Press
,
2003
), p.
38
.
6.
A. V.
Arefiev
,
V. N.
Khudik
, and
M.
Schollmeier
,
Phys. Plasmas
21
,
033104
(
2014
).
7.
A. V.
Arefiev
,
B. N.
Breizman
,
M.
Schollmeier
, and
V. N.
Khudik
,
Phys. Rev. Lett.
108
,
145004
(
2012
).
8.
A. P. L.
Robinson
,
A. V.
Arefiev
, and
D.
Neely
,
Phys. Rev. Lett.
111
,
065002
(
2013
).
9.
S.-W.
Bahk
,
P.
Rousseau
,
T.
Planchon
,
V.
Chvykov
,
G.
Kalintchenko
,
A.
Maksimchuk
,
G.
Mourou
, and
V.
Yanovsky
,
Opt. Lett.
29
,
2837
(
2004
).
10.
See www.eli-laser.eu for Extreme Light Infrastructure European Project.
11.
S. A.
Gaillard
,
T.
Kluge
,
K. A.
Flippo
,
M.
Bussmann
,
B.
Gall
,
T.
Lockard
,
M.
Geissel
,
D. T.
Offermann
,
M.
Schollmeier
,
Y.
Sentoku
, and
T. E.
Cowan
,
Phys. Plasmas
18
,
056710
(
2011
).
12.
T.
Kluge
,
S. A.
Gaillard
,
K. A.
Flippo
,
T.
Burris-Mog
,
W.
Enghardt
,
B.
Gall
,
M.
Geissel
,
A.
Helm
,
S. D.
Kraft
,
T.
Lockard
,
J.
Metzkes
,
D. T.
Offermann
,
M.
Schollmeier
,
U.
Schramm
,
K.
Zeil
,
M.
Bussmann
, and
T. E.
Cowan
,
New J. Phys.
14
,
023038
(
2012
).
13.
A. G.
Krygier
,
D. W.
Schumacher
, and
R. R.
Freeman
,
Phys. Plasmas
21
,
023112
(
2014
).
14.
A.
Zhidkov
,
J.
Koga
,
A.
Sasaki
, and
M.
Uesaka
,
Phys. Rev. Lett.
88
,
185002
(
2002
).
15.
L. L.
Ji
,
A.
Pukhov
,
I. Yu.
Kostyukov
,
B. F.
Shen
, and
K.
Akli
,
Phys. Rev. Lett.
112
,
145003
(
2014
).
16.
A. M.
Fedotov
,
N. B.
Narozhny
,
G.
Mourou
, and
G.
Korn
,
Phys. Rev. Lett.
105
,
080402
(
2010
).
17.
B. B.
Godfrey
and
J.-L.
Vay
,
J. Comput. Phys.
248
,
33
(
2013
).
18.
J.-L.
Vay
,
Phys. Plasmas
15
,
056701
(
2008
).
19.
T.
Esirkepov
,
Comput. Phys. Commun.
135
,
144
(
2001
).
20.
B.
Cohen
,
Multiple Time Scales
(
Academic Press
,
1985
), p.
311
.
21.
G.
Chen
and
L.
Chacon
,
Comput. Phys. Commun.
185
,
2391
(
2014
).
22.
D. R.
Welch
,
D. V.
Rose
,
M. E.
Cuneo
,
R. B.
Campbell
, and
T. A.
Mehlhorn
,
Phys. Plasmas
13
,
063105
(
2006
).
23.
G.
Chen
,
L.
Chacon
, and
D. C.
Barnes
,
J. Comput. Phys.
230
,
7018
(
2011
).