2.1. Analytical model

Considering one electron in a vacuum, we can write its motion equation in terms of the electromagnetic vector potential $\boldsymbol {A}$ as

(2.1)\begin{equation} \frac{{\rm d} \boldsymbol{p}}{{\rm d}t} ={-}e \left[ - \frac{\partial \boldsymbol{A}}{\partial t} + \boldsymbol{v} \times \left( \boldsymbol{\nabla} \times \boldsymbol{A}\right) \right], \end{equation}

where $\boldsymbol {p}$ is the electron momentum, $\boldsymbol {v}$ its velocity, $e$ the elementary charge and the scalar potential $\phi = 0$.

For a plane wave propagating along the $x$ axis, the vector potential depends only on $x$ and is polarised in the $yz$ plane. Assuming an electron initially at rest, we thus get to

(2.2)\begin{gather} \frac{{\rm d}p_x}{{\rm d}t}={-} \frac{e^2}{m_e \gamma} \boldsymbol{A} \boldsymbol{\cdot} \frac{\partial \boldsymbol{A}}{\partial x}={-} \frac{e^2}{2m_e \gamma} \frac{\partial A^2}{\partial x}, \end{gather}

(2.3)\begin{gather}\boldsymbol{p}_\perp{=} e \boldsymbol{A}, \end{gather}

with $\boldsymbol {p}_\perp = p_y \boldsymbol {e_y} + p_z \boldsymbol {e_z}$, $m_e$ being the electron mass and $\gamma$ its Lorentz factor. We can recognise in (2.2) the expression for the relativist ponderomotive force in the one-dimensional case.

From the conservation of energy equation for the electron and using the fact that for a plane wave propagating at $c$ along the $x$ axis, $( {\partial }/{\partial t} + c ({\partial }/{\partial }x) ) A = 0$, we can then show that the quantity $\gamma - p_x / m_e c$ is conserved. With the electron being initially at rest, we have $\gamma (t=0)=\gamma _0=1$ and $p_x(t=0)=p_{x,0}=0$ which leads to

(2.4)\begin{equation} \gamma = 1 + \frac{p_x}{m_e c}. \end{equation}

We thus have three equations describing the dynamic of an initially at rest electron in a propagating plane wave:

(2.5)\begin{gather} \frac{p_x}{m_ec} = \frac{1}{2} \left( \frac{e \boldsymbol{A}}{m_e c} \right)^2, \end{gather}

(2.6)\begin{gather}\frac{\boldsymbol{p}_\perp}{m_ec} = \frac{e \boldsymbol{A}}{m_e c}, \end{gather}

(2.7)\begin{gather}\gamma = 1 + \frac{1}{2} \left( \frac{e \boldsymbol{A}}{m_e c} \right)^2 . \end{gather}

2.2. Numerical simulations

All numerical simulations presented in this article were performed using the code Calder (Lefebvre et al. Reference Lefebvre, Cochet, Fritzler, Malka, Aléonard, Chemin, Darbon, Disdier, Faure and Fedotoff2003). Before we delve into the results of those simulations though, we first discuss the computation of electromagnetic fields and their resulting actions in PIC codes in order to introduce and summarise the B-TIS method proposed in Bourgeois & Davoine (Reference Bourgeois and Davoine2020).

2.2.1. Temporal interpolation of the magnetic field $B$ in PIC codes

2.2.1.1. Standard temporal interpolation

Although there have been interesting new developments to create dispersionless Maxwell solvers in recent years, notably most recently (Pukhov Reference Pukhov2020), Calder, as many other PIC codes, still uses the simple and robust finite-difference time-domain Yee scheme (Yee Reference Yee1966) to solve the Maxwell equations. Electric ($E$) and magnetic ($B$) fields are thus calculated using a leap-frog method: they are defined on different grids which are offset spatially by half a cell but also temporally by $\Delta t / 2$.

Let us consider the computation of the transverse Lorentz force applied to a moving particle in our simulation. The particle can move freely in the whole simulated space while the fields are computed only on specific grid points. It is thus necessary to interpolate the electromagnetic fields onto the charged particle position to calculate its motion.

This simple situation, with the field amplitude values known on the different points of the grid and an electron moving freely between those points, is shown in figure 1 for a simple one-dimensional grid. Considering only the longitudinal dimension $x$ to simplify the notations, $A^{n}_{i}$ then denotes the value of the field $A$ at the $n$th time step and the $i$th point of the spatial grid or in other terms $A^{n}_{i}= A (t= n\Delta t,x=i\Delta x )$. Due to the leap-frog nature of the Yee scheme, values of the $E_y$ fields are known at integer time steps and grid points while values for $B_z$ are known at half-integer time steps and grid points. In the remainder of this section, we only consider $E_y$ and $B_z$ – which are required to compute the transverse force $F_y$ – and then omit the $y$ and $z$ indices when it is not absolutely necessary so as not to clutter the notation.

Figure 1. Usual interpolation process of electromagnetic fields on a one-dimensional spatial grid at the $n$th time step in order to get the transverse force $F_y^n$ applied on the particle represented by the yellow circle. (a) The initial configuration where the field values are known at different times and spatial points. The temporal interpolation to get $\tilde {B}^n$ is shown with dashed arrows. (b) The subsequent spatial interpolation of $E^n$ and $\tilde {B}^n$ on the particle (again with dashed arrows and first order).

As the fields are initially not known at the same time step, a temporal interpolation is necessary to get $\tilde {B}^n$ before the spatial interpolation and the computation of the force which is done at the integer time step $n$. The simplest and most common method is the linear time interpolation (LTI) $\tilde {B}^n = ( B^{n-{1}/{2}} + B^{n+{1}/{2}} )/2$ as is shown in figure 1(a) but a quadratic time interpolation (QTI), with

(2.8)\begin{equation} \tilde{\tilde{B}}_{i+{1}/{2}}^{n} = \frac{3}{8} B_{i+{1}/{2}}^{n+{1}/{2}} + \frac{3}{4} B_{i+{1}/{2}}^{n-{1}/{2}} - \frac{1}{8} B_{i+{1}/{2}}^{n-{3}/{2}}, \end{equation}

can sometimes be used instead for higher accuracy (Lehe et al. Reference Lehe, Thaury, Guillaume, Lifschitz and Malka2014). Once both fields $E$ and $B$ are known at the same time step, they are interpolated spatially onto the particles. This step is shown in figure 1(b). Once again different kinds of spatial interpolation can be used with varying orders of interpolation. The one depicted here is, for simplicity's sake, a first-order method (linear interpolation).

2.2.1.2. Bypassing the need for a temporal interpolation: B-TIS

The idea behind our new scheme proposed in Bourgeois & Davoine (Reference Bourgeois and Davoine2020) is to simply use the available value of the magnetic field as the value we need. That is, interpolating the fields at the particle position using ($E_i^{n}$, $\hat {B}_{i}^{n}= B_{i+{1}/{2}}^{n+{1}/{2}}$) instead of ($E_i^{n}$, $\tilde {B}_{i+{1}/{2}}^{n}$) as previously shown. Of course, doing so is only physically sound if $B_{i}^{n} \approx B_{i+{1}/{2}}^{n+{1}/{2}}$. Thankfully this condition means for a wave propagating at $c$ that $B(x,t) \approx B(x+\Delta x/2, t + \Delta t/2)$ which is true as long as $\Delta x \approx c \Delta t$. Although this last condition can be limited in the simulation by the Courant–Friedrichs–Lewy condition, it is usually verified in PIC simulations as having a ratio $\Delta x / \Delta t$ close to $c$ improves the quality of the results.

This new method effectively results in a translation of the calculated $B$ field before the spatial interpolation, hence its name: B-translated interpolation scheme. This process is described in figure 2.

Figure 2. Modified interpolation process of electromagnetic fields: B-TIS on a one-dimensional grid at the $n$th time step. (a) The same initial configuration as figure 1(a) but instead of the temporal interpolation of $B$, the relation $\hat {B}^n_i = B^{n+{1}/{2}}_{i+{1}/{2}}$ is used. (b) The result of this translation and the subsequent spatial interpolation of $E^n$ and $\hat {B}^n_i$ on a particle.

2.2.1.3. Computational error introduced

Both approaches introduce an error on the computed values of the magnetic field compared with an ideal situation where the magnetic field can be computed at integer time steps.

Let us consider a given magnetic field $B(\varphi ) = B_0 \cos (\varphi )$. The computation of the numerical value of that field in the simulation introduces an error $\varepsilon$ which depends on the method used:

(2.9)\begin{gather} \varepsilon_{\text{LTI}} ={-} B_0 \cos (\varphi) \left(1 - \cos(\delta\varphi) \right), \end{gather}

(2.10)\begin{gather}\varepsilon_{\text{QTI}} ={-} B_0 \cos (\varphi) \left(1 - \cos(\delta\varphi) \left(1 + \frac{1}{2} \sin^2(\delta\varphi) \right) \right) - \frac{1}{2} B_0 \sin(\varphi) \sin^3(\delta\varphi), \end{gather}

(2.11)\begin{gather}\varepsilon_{\text{B-TIS}} ={-} B_0 \cos (\varphi) \left(1 - \cos(\widehat{\delta \varphi}) \right) - B_0 \sin(\varphi) \sin(\widehat{\delta \varphi}), \end{gather}

with $\varphi = k_x x + k_y y + k_z z - \omega t$, $\delta \varphi = {\omega \Delta t}/{2}$ and $\widehat {\delta \varphi } = (k_x \Delta x - \omega \Delta t)/2$. Here $\omega$ is the field's angular frequency and $k_x$, $k_y$, $k_z$ its wavevector components.

For longitudinally propagating waves with $\omega /k_x \approx c$ and $c \Delta t \approx \Delta x$, $\widehat {\delta \varphi }$ can become very small and the gain in accuracy should then be significant. This improvement of the computation of the magnetic field is what enables B-TIS to effectively suppress the effects of NCR. Indeed, the more accurate computation result in the induced Lorentz force, from the interaction of the NCR with particles, being negligible. In effect, the spurious radiation is still generated and present in the simulation but has no impact on the simulated particles.

2.2.2. Comparison of the different temporal interpolation methods

Three simulations were originally performed, each with a different temporal interpolation method introduced earlier: LTI, QTI and B-TIS. To reproduce in a two-dimensional simulation the situation presented above, we use a plane wave propagating along the $x$ axis, infinite in width along the $y$ axis. Periodical boundary conditions are used on the longitudinal edges of the simulation box to recreate the infinite transverse property of the wave.

We chose the polarisation direction in the simulation plane (along the $y$ axis). The laser wavelength is $\lambda _0 = 2 {\rm \pi}c / \omega _0$, with $\omega _0$ the laser angular frequency, and the wave amplitude is characterised by $a_0 = e A_0 / m_e c = 5$. Simulations were performed in a moving window following the wave propagation of $960 \times 200$ cells with $\Delta x = 0.15 c/\omega _0$, $\Delta y = 10 \ c/\omega _0$ and a time step $\Delta t = 0.149 \ 1/\omega _0$.

The propagating wave's temporal profile has a finite extension with a trapezoidal envelope as shown in figure 3(a). Each ramp and the plateau of the trapezoidal profile are $4 \lambda _0$ long. We take care that $\int E_y(t)\,{\rm d}t = 0$, otherwise $A_y$, and thus $p_y$, is non-zero even once the laser has passed. To make the analysis easier we also chose the temporal profile so that $\int _{t_1}^{t_2} E_y(t) \,{\rm d}t = 0$, with $t_1$ and $t_2$ being the time boundaries of the profile ramps or plateau. This ensures that the transverse impulsion $p_y$ is zero on average during the whole passing of the laser intensity plateau by the electron.

Figure 3. (a) Temporal profile (normalised vector potential $a$) of the incident wave. The dashed line describes the trapezoidal envelope. (b) Initial situation of the simulation. The incident wave is propagating to the right.

Only one particle is considered here. It is initially at rest at the centre of the simulation box. This single particle is given a very small statistical weight in the simulation so that the charge density (computed within each cell of the simulation) around the particle is negligible and does not affect the particle dynamic. Figure 3(b) shows the initial state of the simulation. It is an electron density map (normalised to the plasma critical density $n_c = \epsilon _0 m_e \omega _0^2 / e^2$) showing the initial position of the particle superimposed on a map of the wave transverse electric field $E_y$.

Considering a propagating plane wave along $x$, polarised along the $y$ axis and with a temporal envelope $f$ such as $\boldsymbol {A} = A_0 f(t - x/c) \sin (\omega _0 (t - x/c) ) \boldsymbol {e_y}$, then (2.5), (2.6) and (2.7) may be rewritten as

(2.12)\begin{gather} \frac{p_x}{m_ec} = \frac{1}{2} a_0^2 f^2(t - x/c) \sin^2 \! \left(\omega_0 (t - x/c) \right), \end{gather}

(2.13)\begin{gather}\frac{p_y}{m_ec} = a_0 f(t - x/c) \sin \! \left(\omega_0 (t - x/c) \right), \end{gather}

(2.14)\begin{gather}\gamma = 1 + \frac{1}{2} a_0^2 f^2(t - x/c) \sin^2 \! \left(\omega_0 (t - x/c) \right). \end{gather}

Figure 4 shows results for each simulation by comparing the particle momenta $p_x$ and $p_y$ with those expected according to theoretical computations based on (2.12) and (2.13). These figures use normalised values, meaning $t$ is expressed as multiples of $1/\omega _0$, $x$ as multiples of $c/\omega _0$ and the momentum $p$ as multiples of $m_e c$.

Figure 4. Evolution of the electron normalised momenta $p_x$ (a) and $p_y$ (b) with respect to the normalised coordinate $ct-x$ for the considered different methods. The red dashed line shows the theoretical $p_x$ (respectively $p_y$) given by (2.12) (respectively (2.13)).

The QTI and B-TIS methods give values of $p_y$ closer to the theoretical value than LTI but, overall, all three simulations reproduce fairly well the transverse behaviour of the particle with small differences.

Still, notable differences between analytical and simulation results though are present at the beginning and end of the trapezoid ramps ($t-x \approx 20$ and $t-x \approx 90$) in figure 4(b). Those are due to a small difference in how the laser temporal profile is defined: in the theoretical case we considered a trapezoidal envelope for the vector potential $A$, whereas we considered a trapezoidal envelope for the electric and magnetic fields $E$ and $B$ in the simulation, as it is easier to initialise the value of those fields. This definition gives rise to vector potential $A$ with a slightly different shape when $A$ varies significantly within a single wavelength, as is the case at the boundaries of the temporal profile of the wave.

The longitudinal momentum is more problematic and only the LTI method reproduces qualitatively the expected behaviour of the particle, with a notable error on the maximal value of $p_x$ which is underestimated. Both QTI and B-TIS lead to an overestimation of $p_x$ with an error increasing over time (artificial increase in the energy of the particle). This error leads to a non-physical residual momentum after the wave has gone by.

In the next section we investigate the source of this error and present a possible solution for the B-TIS method.

2.3. Inaccurate longitudinal momentum with QTI and B-TIS

The laser wave propagation in Calder is computed using the fields $E$ and $B$ not the vector potential $A$ so we will use here expressions with respect to the electromagnetic fields. The electron motion equation, for the longitudinal momentum, is thus

(2.15)\begin{equation} \frac{{\rm d}p_x}{{\rm d}t} ={-}e \left( E_x + v_y B_z - v_z B_y \right). \end{equation}

Considering a plane wave propagating along $x$ and polarised along the $y$ axis such that $E_y = E_0 \cos ( \varphi )$ with $\varphi = \omega _0 (t-x/c)$, then $B_z = (E_0 / c) \cos ( \varphi )$ and $E_x = 0$, $B_y =0$. The previous equation thus becomes

(2.16)\begin{equation} \frac{{\rm d}p_x}{{\rm d}t} ={-} \frac{e E_0}{c} v_y \cos ( \varphi ). \end{equation}

From (2.6), we get that $v_y \propto A(\varphi ) \propto \sin ( \varphi )$, which leads to

(2.17)\begin{equation} p_x(t) = \int_0^t \frac{{\rm d}p_x}{{\rm d}t} dt \propto \int_0^t \sin ( \varphi ) \cos ( \varphi ) \,{\rm d}\varphi. \end{equation}

This means that for an electron initially at rest, after an integer number of wave periods, the momentum $p_x$ should be $0$ again as $\int _0^{2{\rm \pi} } \sin ( \varphi ) \cos ( \varphi ) \,{\rm d}\varphi = 0$. This is indeed what is observed in figure 4 for the theoretical curve and LTI but not for QTI and B-TIS. This error comes from the computational error introduced during the magnetic field interpolation step.

Looking at (2.9)–(2.11), for all three methods, we can see that the error depends on a term that is proportional to $\cos (\varphi )$ but in the case of QTI and B-TIS there is a second term proportional to $\sin (\varphi )$. It is this last part which is problematic here. Taking a look at the numerically computed momentum $p_x^{{\rm num}}$ in our simulation, we can express it as

(2.18)\begin{equation} p_x^{{\rm num}}(t) = p_x(t) + \int_0^t - \frac{e E_0}{c} v_y \, \varepsilon \, {\rm d}t = p_x(t) + \varepsilon_p, \end{equation}

where $\varepsilon _p$ is the numerical error on the momentum $p_x$. Thus, with LTI, we get that

(2.19)\begin{equation} \varepsilon_p^{\text{LTI}} = K \left( 1 - \cos ( \delta \varphi ) \right) \int_0^t \sin ( \varphi ) \cos ( \varphi ) \,{\rm d}\varphi ,\end{equation}

with $K$ a constant depending on the wave amplitude. This error will fluctuate with time but it will periodically cancel out, whereas with QTI and B-TIS we get

(2.20)\begin{gather} \varepsilon_p^{\text{QTI}} = K' \int_0^t \sin ( \varphi ) \cos ( \varphi ) \,{\rm d}\varphi + \frac{1}{2} K \sin^3 ( \delta \varphi ) \int_0^t \sin^2 ( \varphi ) \,{\rm d}\varphi, \end{gather}

(2.21)\begin{gather}\varepsilon_p^{\text{B-TIS}} = K'' \int_0^t \sin ( \varphi ) \cos ( \varphi ) \,{\rm d}\varphi + K \sin ( \widehat{\delta \varphi} ) \int_0^t \sin^2 ( \varphi ) \,{\rm d}\varphi, \end{gather}

with $K' = K ( 1 - \cos ( \delta \varphi ) ( 1 + \tfrac {1}{2} \sin ^2 ( \delta \varphi ) ) )$ and $K'' = K ( 1 - \cos ( \widehat {\delta \varphi } ) )$. The second term in these expressions leads to an error that grows continuously and monotonically with time instead of periodically cancelling out. This explains the divergence of results observed in figure 4 though the initial error is really small. Note also that, with our values for the parameters $\Delta x$ and $\Delta t$, $\tfrac {1}{2} \sin ^3 ( \delta \varphi ) \approx 2\times 10^{-4}$ whereas $\sin ( \widehat {\delta \varphi } ) \approx 5\times 10^{-4}$. This explains the difference in growth rate for the error with QTI and B-TIS.

The obvious solution would be to use a finer grid so as to minimise $\Delta t$ and $\Delta x - c \Delta t$ – thus minimising $\delta \varphi$, $\widehat {\delta \varphi }$ and the induced error – but this is highly detrimental to computation performance. Restricting oneself to using only small values of $a_0$ is another way to circumvent this problem – as the error is proportional to $a_0^2$ – but it is problematic in the case of the study of VLA or DLA. In any case, both of those approaches only mitigate the problem as the error will still be present and, though small initially, it will compound and increase with time.

It is, however, possible to slightly modify B-TIS in a way that drastically reduces this source of error while still benefiting from all of the advantages of this method.

2.4. B-TIS modification

As we saw earlier, B-TIS replaces the following temporal interpolation of the magnetic field:

(2.22)\begin{equation} \tilde{B}_{i+{1}/{2}}^{n} = \frac{1}{2} \left(B_{i+{1}/{2}}^{n+{1}/{2}} + B_{i+{1}/{2}}^{n-{1}/{2}} \right) \end{equation}

with the following translation (as shown in figure 2a):

(2.23)\begin{equation} \hat{B}_i^n = B_{i+{1}/{2}}^{n+{1}/{2}}. \end{equation}

This is based on the assumption that the wave propagates at $c$ along the $x$ axis and that $c \Delta t \approx \Delta x$.

However, if the value of the magnetic field is conserved going forward in time and space, it is also conserved when going backward, and within those assumptions, that translation is completely equivalent to

(2.24)\begin{equation} \hat{B}_i^n = B_{i-{1}/{2}}^{n-{1}/{2}}. \end{equation}

Both possible translations are shown in figure 2(b) where the assumption that the magnetic field $B$ propagates at $c$ and that $c \Delta t \approx \Delta x$ means that $B$ is constant along the diagonals of the grid, hence the equivalence of the two translations.

Let us call rB-TIS the ‘reversed’ version of B-TIS using $\hat {B}_i^n = B_{i-{1}/{2}}^{n-{1}/{2}}$. This method then introduces the following error on the magnetic field:

(2.25)\begin{align} \varepsilon_{\text{rB-TIS}} & ={-} B_0 \cos (\varphi) \left(1 - \cos(-\widehat{\delta\varphi}) \right) - B_0 \sin(\varphi) \sin(-\widehat{\delta\varphi}), \end{align}

(2.26)\begin{align} & ={-} B_0 \cos (\varphi) \left(1 - \cos(\widehat{\delta\varphi}) \right) + B_0 \sin(\varphi) \sin(\widehat{\delta\varphi}). \end{align}

With the exception of the sign of the second term, this is the same expression as the error for the ‘classical’ B-TIS. As such, we should be able to eliminate that second term by combining both B-TIS and rB-TIS. We can indeed choose to use both translations in conjunction at each time step through a simple average such as

(2.27)\begin{equation} \hat{B}_i^n = \frac{1}{2} \left( B_{i+{1}/{2}}^{n+{1}/{2}} + B_{i-{1}/{2}}^{n-{1}/{2}}\right). \end{equation}

It is important to note, though, that it is not a simple temporal averaging as is used in LTI: both the time index and the spatial index are different. This new method, which we call B-TIS3, then introduces the following error on $B$:

(2.28)\begin{align} \varepsilon_{\text{B-TIS3}} & = \frac{1}{2} \varepsilon_{\text{B-TIS}} + \frac{1}{2} \varepsilon_{\text{rB-TIS}}, \end{align}

(2.29)\begin{align} & ={-} B_0 \cos (\varphi) \left(1 - \cos(\widehat{\delta\varphi}) \right). \end{align}

The second term involving $\sin (\varphi )$ has indeed disappeared as we expected. The expression of this error is similar to the one introduced by LTI with the important difference though that we have $\widehat {\delta \varphi } < \delta \varphi$, which should point to an improvement of accuracy.

To illustrate this point, let us consider a forward-propagating wave in a vacuum with $\omega = \omega _0$ and $k_x = k = {\omega _0}/{c}$. We thus have $\delta \varphi = {\omega _0 \Delta t}/{2}$ while $\widehat {\delta \varphi } = \tfrac {1}{2} ({\omega _0}/{c}) ( \Delta x - c \Delta t )$.

Then, per (2.9), (2.10), (2.11) and (2.29), with our numerical parameters we get

(2.30)\begin{gather} \varepsilon^{\text{las}}_{\text{LTI}} ={-} B_0 \cos (\varphi) \times 2.77\times10^{{-}3}, \end{gather}

(2.31)\begin{gather}\varepsilon^{\text{las}}_{\text{LTI}} ={-} B_0 \cos (\varphi) \times 1.15\times 10^{{-}5} - B_0 \sin (\varphi) \times 2.06\times 10^{{-}4}, \end{gather}

(2.32)\begin{gather}\varepsilon^{\text{las}}_{\text{B-TIS}} ={-} B_0 \cos (\varphi) \times 1.25\times10^{{-}7} - B_0 \sin (\varphi) \times 5\times 10^{{-}5}, \end{gather}

(2.33)\begin{gather}\varepsilon^{\text{las}}_{\text{B-TIS3}} ={-} B_0 \cos (\varphi) \times 1.25\times10^{{-}7} . \end{gather}

On the other hand, let us consider a wave, propagating forward in the simulation at the Nyquist frequency so that $\omega = {{\rm \pi} }/{\Delta t}$ and $k_x = k = {{\rm \pi} }/{\Delta x}$. Then, as $\delta \varphi = {{\rm \pi} }/{2}$ and $\widehat {\delta \varphi } = 0$:

(2.34)\begin{gather} \varepsilon^{\text{Nyq}}_{\text{LTI}} ={-} B_0 \cos (\varphi), \end{gather}

(2.35)\begin{gather}\varepsilon^{\text{las}}_{\text{LTI}} ={-} B_0 \cos (\varphi) - B_0 \sin (\varphi), \end{gather}

(2.36)\begin{gather}\varepsilon^{\text{las}}_{\text{B-TIS}} = 0, \end{gather}

(2.37)\begin{gather}\varepsilon^{\text{Nyq}}_{\text{B-TIS3}} = 0. \end{gather}

Both of these simple cases show the immense improvement in accuracy that B-TIS3 brings to our simulations regarding the magnetic field over all the previously considered methods.

A comparison of the results from the B-TIS3 method and the previous one is presented in figure 5. We can see that this improved B-TIS3 eliminates the troublesome compounding error impacting B-TIS1 and leads to computed values of $p_x$ even closer to the theoretical ones than those obtained using LTI. Results on transverse momentum $p_y$ are just as good as before, the improvement margin being already very slim even with LTI.

Figure 5. Evolution of the electron normalised momentum $p_x$ (a) and $p_y$ (b) with respect to the normalised coordinate $ct-x$ for the considered different methods. The red dashed line shows the theoretical $p_x$ (respectively $p_y$) given by (2.12) (respectively (2.13)).

At this point, the reader should note that we have focused on forward-propagating waves, as they are those that are most commonly encountered in VLA and LWFA and for which this scheme was designed. Counter-propagating radiation, however, can be present in a few situations, because of backscattering for instance or even in the case of colliding laser pulses. It is thus of interest to check the robustness of our new scheme in such scenarios.

The accuracy of both LTI and QTI is not impacted for a counter-propagating wave, as we still have $\delta \varphi = {\omega \Delta t}/{2}$. For B-TIS and B-TIS3, however, with $k_x = -k$, we now get $\widehat {\delta \varphi } = ( k \Delta x + \omega \Delta t )/2$. With this, for a counter-propagating laser pulse we get

(2.38)\begin{equation} \varepsilon^{\text{back las}}_{\text{B-TIS3}} ={-} B_0 \cos (\varphi) \times 1.11\times10^{{-}2}. \end{equation}

And for a wave propagating backward at the Nyquist frequency:

(2.39)\begin{equation} \varepsilon^{\text{back Nyq}}_{\text{B-TIS3}} ={-} 2 B_0 \cos (\varphi). \end{equation}

While these results appear much worse than in the case of forward propagation, it should be noted that they are not that different from the accuracy level of LTI, with a factor of 2 to 4 between them. Therefore, the presence of counter-propagating radiation in the simulation, although not the ideal use case for B-TIS3, should not prove physics-breaking and the overall gains might be well worth this trade-off.

From these results, it appears that we have successfully solved the initial problem encountered in the computation of the longitudinal momentum $p_x$ and that B-TIS3 appears well suited to reproduce the physics of VLA contrary to B-TIS1 and QTI.

As LTI is by far the most common method used in PIC codes, it seems an appropriate method with which to compare. We then only present simulation results for both of these methods. A summary of results for LTI and B-TIS3 is presented in figure 6. The better reproduction of the momentum with B-TIS3 shown in figure 6(a,b) leads to an electron trajectory more faithful to the theoretical predictions as can be seen in figure 6(c,d).

Figure 6. Evolution of the electron normalised momentum (a) $p_x$ and (b) $p_y$ and normalised position (c) $x$ and (d) $y$ with respect to the normalised coordinate $ct-x$ for LTI and B-TIS3. The red dashed line shows the theoretical values expected from the analytical model introduced earlier.

2.5. Simulations with initial longitudinal momentum

We have so far investigated only the case of an electron initially at rest but, in practice, electrons are far more likely to have a non-zero initial momentum. Going back to the analytical model developed in § 2.1, with $p_{x,0} \neq 0$ and $\gamma _0 \neq 1$ we now get

(2.40)\begin{equation} \gamma = \gamma_0 + \frac{p_x}{m_e c} - \frac{p_{x,0}}{m_ec}, \end{equation}

which leads to the three following equations now describing the electron motion:

(2.41)\begin{gather} \frac{p_x}{m_ec} = \frac{p_{x,0}}{m_ec} + \frac{1}{2} \left( \gamma_0 + \frac{p_{x,0}}{m_ec} \right) \left( \frac{e \boldsymbol{A}}{m_e c} \right)^2, \end{gather}

(2.42)\begin{gather}\frac{\boldsymbol{p}_\perp}{m_ec} = \frac{e \boldsymbol{A}}{m_e c}, \end{gather}

(2.43)\begin{gather}\gamma = \gamma_0 + \frac{1}{2} \left( \gamma_0 + \frac{p_{x,0}}{m_ec} \right) \left( \frac{e \boldsymbol{A}}{m_e c} \right)^2. \end{gather}

Figure 7 shows results from two simulations where the electron had an initial longitudinal momentum $p_0$: the first one with $p_0 = m_e c$ and the second one with $p_0 = 5 m_e c$. All other numerical and physical parameters are the same as previously. The two methods result in much more noticeable differences than previously.

Figure 7. Evolution of the electron normalised momentum (a) (respectively c) $p_x$ and (b) (respectively d) $p_y$ with respect to the normalised coordinate $ct-x$ for LTI and B-TIS3 with $p_0 = m_e\, c$ (respectively $p_0 = 5\, m_e\, c$). The red dashed line shows the theoretical values expected from the analytical model introduced earlier.

It is apparent that the greater the initial velocity of the electron, the more LTI underestimates the longitudinal momentum. The maximum longitudinal momentum $p_{x,max}$ being underestimated by a factor of at least $3$ compared with the theoretical value when using LTI with an initial momentum $p_0 = 5\, m_e\, c$. B-TIS3, on the other hand, slightly overestimates $p_x$ but overall produces results much closer to the predictions of the analytical model.

Looking now at the transverse momentum $p_y$, both methods reproduce fairly well the theoretical results, though both tend to overestimate them. Once again, the B-TIS3 results appear closer to those of the model and, though the error increases when the initial momentum $p_0$ increases, this increase appears much smaller for B-TIS3 than for LTI.

Note that all of the presented simulations have been realised using the Boris pusher (Boris Reference Boris1970). We found very negligible differences between simulations using the Boris scheme and simulations using the scheme proposed in Vay (Reference Vay2008) to correct the Boris pusher inaccuracies at relativistic speed, and we observed the same inaccuracies with LTI and the same improvements with B-TIS3, irrespective of the choice of pusher.

3.1. Theoretical description

We keep the same trapezoidal temporal profile for our laser pulse so as to make comparisons with earlier results easier; however, its transverse size is now finite and given by a Gaussian profile with a width (full width at half maximum of the intensity profile) of $160 \omega _0/c$. Figure 8 shows this initial situation.

Figure 8. Initial situation with a spatially Gaussian wave.

We can come back to the equations introduced in § 2.1 but considering now the more general case of a wave polarised along the $y$ axis and propagating along the $x$ axis such that $\boldsymbol {A}(x,y,t) = A(x,y,t) \boldsymbol {e}_y$. The laser waist is assumed here much larger than the wavelength and as such we can neglect the existence of a vector potential component along $x$ (we discuss this approximation later). The energy conservation equation then leads to

(3.1)\begin{equation} \frac{{\rm d}}{{\rm d}t}\left( \gamma m c^2 \right) = e v_y \frac{\partial A}{\partial t}. \end{equation}

And the motion equation (2.1) gives us

(3.2)\begin{gather} \frac{{\rm d}p_x}{{\rm d}t} ={-}e v_y \frac{\partial A}{\partial x}, \end{gather}

(3.3)\begin{gather}\frac{{\rm d}}{{\rm d}t}\left(p_y-e A\right) ={-} e v_y \frac{\partial A}{\partial y}. \end{gather}

Using the fact that ${{\rm d} \boldsymbol {A}}/{{\rm d}t} = {\partial \boldsymbol {A}}/{\partial t} + ( \boldsymbol {v} \boldsymbol {\cdot } \boldsymbol {\nabla } ) \boldsymbol {A}$, we finally get to

(3.4)\begin{gather} \frac{{\rm d}p_x}{{\rm d}t} ={-}e v_y \frac{\partial A}{\partial x}, \end{gather}

(3.5)\begin{gather}\frac{{\rm d}p_y}{{\rm d}t} = e \left(\frac{\partial A}{\partial t} + v_x \frac{\partial A}{\partial x} \right),\end{gather}

(3.6)\begin{gather}\frac{{\rm d}}{{\rm d}t}\left( \gamma m c^2 \right) = e v_y \frac{\partial A}{\partial t}. \end{gather}

This system is much more complex than the previous one and we do not attempt to solve it analytically here. We can still gather though that the ponderomotive force of the laser will push aside the electron and that it will gain a residual momentum, both longitudinal and transverse. For an exact solution of the motion and acceleration of electrons by the laser fields, we would have to consider the longitudinal component of $\boldsymbol {A}$ ($\boldsymbol {A}=A_x\boldsymbol {e_x}+A_y\boldsymbol {e_y}$) which would modify (3.1) from (3.6) above. For simplicity's sake, however, we do not consider the component $A_x$ which is negligible in our case. More comprehensive descriptions and computations on this subject can be found in Mora & Antonsen (Reference Mora and Antonsen1997) and Quesnel & Mora (Reference Quesnel and Mora1998).

3.2. Numerical results

Figures 9 and 10 sum up the simulation results for an electron initially close to the propagation axis and, respectively, with no initial longitudinal momentum ($p_0 = 0$) or with $p_0 = m_e c$.

Figure 9. Evolution of the electron normalised momentum (a) $p_x$ and (b) $p_y$ and normalised position (c) $x$ and (d) $y$ with respect to the normalised coordinate $ct-x$ for LTI and B-TIS3 with a transversely Gaussian laser beam and an electron initially at rest ($p_0 = 0$).

Figure 10. Evolution of the electron normalised momentum (a) $p_x$ and (b) $p_y$ and normalised position (c) $x$ and (d) $y$ with respect to the normalised coordinate $ct-x$ for LTI and B-TIS3 with a transversely Gaussian laser beam and an electron with an initial longitudinal momentum ($p_0 = m_e c$).

The behaviour observed in the simulations is in accordance with our expectations: the electron drifts transversely as it is being accelerated and thus ends up with residual non-zero longitudinal and transverse momentum, $p_x$ and $p_y$, even after it has exited the laser beam. This residual extra momentum is quite small when the electron is initially at rest but it quickly becomes more important when the electron has some initial velocity as we can see in figure 10(a,b).

Both methods predict overall a similar behaviour for the electron though we can observe notable differences in final momentum and positions. Indeed as the electron is evidently more accelerated in the B-TIS3 case, it travels a much longer distance while inside the laser pulse. This difference is especially apparent in the case with a non-zero initial momentum, which is coherent with our previous observations.

As we cannot compare with an analytical model here, we use a numerical method to estimate the theoretical result. When reducing both the spatial grid size and the time step duration, both numerical methods should converge on a unique solution corresponding to the physical solution.

We look at the results of three pairs of simulations where we chose $\Delta x = 0.15 c/\omega _0$ and $c \Delta t = 0.149 c/\omega _0$ for the first one, $\Delta x = 0.10 c/\omega _0$ and $c \Delta t = 0.099 c/\omega _0$ for the second and $\Delta x = 0.05 c/\omega _0$ and $c \Delta t = 0.049 c/\omega _0$ for the third. Observed momenta are shown in figure 11. Those results clearly show that both methods are converging on the same solution which appears fairly close to the one obtained initially through B-TIS3. We thus observe once again that B-TIS3 tends to slightly overestimate the electron momentum while LTI noticeably underestimates it.

Figure 11. Evolution of the electron normalised momentum (a) $p_x$ and (b) $p_y$ with respect to the normalised coordinate $ct-x$ for LTI and B-TIS3, with $p_0 = \, m_e c$ and varying cell sizes and time steps.

It is important to note, however, that even when dividing by 3 both $\Delta t$ and $\Delta x$ (thus increasing the computation time by a factor of 9), LTI still does not give better results than B-TIS3 with a worse resolution. B-TIS3 appears, in that light, as a fairly good improvement on the standard LTI.

3.3. Influence of the initial position

The ponderomotive force acting on a particle depends on the laser pulse intensity gradient which varies widely depending on the transverse position.

We show in this section that we can indeed observe this behaviour in our simulations. We present here in figure 12 the momenta of electrons with different initial positions $y_0$.

Figure 12. Evolution of the electron normalised momentum (a) $p_x$ and (b) $p_y$ with respect to the normalised coordinate $ct-x$ for LTI and B-TIS3, with $p_0 = \, m_e c$ and varying initial transverse position $y_0$ of the electron. Note that the transverse profile of the laser is centred on $y=0$.

Both methods tend to predict similar final momentum for the electron as long as it is outside of the stronger field region of the laser pulse. Indeed differences in results between LTI and B-TIS3 are less important the farther is the electron initially from the propagation axis.

As we showed that the error introduced by LTI was greater when the electron was strongly accelerated, it thus seems logical that its predictions are better – thus more in accordance with B-TIS3 – when the electron is far from the high-intensity region near the propagation axis and less accelerated.

Though final momentum may appear similar, the small differences can have a great impact on the particle trajectory. Furthermore, when trying to study VLA we are of course more interested in the strongly accelerated electrons. Thus the usage of the B-TIS3 method appears better suited to these simulations than the standard LTI.