Overcoming timestep limitations in boosted-frame particle-in-cell simulations of plasma-based acceleration

Explicit electromagnetic particle-in-cell (PIC) codes are typically limited by the Courant-Friedrichs-Lewy (CFL) condition, which implies that the timestep multiplied by the speed of light must be smaller than the smallest cell size. In the case of boosted-frame PIC simulations of plasma-based acceleration, this limitation can be a major hindrance, as the cells are often very elongated along the longitudinal direction and the timestep is thus limited by the small, transverse cell size. This entails many small-timestep PIC iterations and can limit the potential speed-up of the boosted-frame technique. Here, by using a CFL-free analytical spectral solver, and by mitigating additional numerical instabilities that arise at large timestep, we show that it is possible to overcome traditional limitations on the timestep and thereby realize the full potential of the boosted-frame technique over a much wider range of parameters.

Figure 1

NCI growth rate of the Galilean PSATD scheme. Normalized NCI growth rate $\text{Im}\left(\omega \right)/{\omega }_{p,r}$ in spectral ($k_x,k_z$) space, calculated from the analytical stability analysis for different Galilean velocities $v_{\text{gal}}=v_0$ (a), (b) and $v_{\text{gal}}=0.99v_0$ (c), (d), and for different timesteps $c\mathrm{\Delta }t=\mathrm{\Delta }x$ (left) and $c\mathrm{\Delta }t=6\mathrm{\Delta }x$ (right). The solid and dotted lines correspond to well-known aliased NCI resonant modes, with alias number $(m_z,n)$, as given by Eq. (1). In this simulation, a uniform plasma drifts at a velocity $v_0=c{(1-1/{\gamma }_b^2)}^{1/2}$ with ${\gamma }_b=130$, and the transverse and longitudinal cell sizes are $\mathrm{\Delta }x=6.4\times {10}^{-2}{k}_{p,r}^{-1}$ and $\mathrm{\Delta }z=6\mathrm{\Delta }x$, respectively (where $k_{p,r}^2=n_0{e}^2/\left(m_e{\epsilon }_0{c}^2{\gamma }_0\right)$, and where $n_0$ is the plasma density).

Figure 2

Illustration of the field push and particle momenta push, in the averaged Galilean PSATD algorithm. The quantities represented in solid frames are the ones that are known just before the field and particle push. Note that this includes the deposited charge and current $\widehat{\rho }$ and $\widehat{\mathcal{J}}$. The quantities in dashed and dotted frames are the ones that are being computed during the field and particle push. As part of the field push, the regular fields $\widehat{\mathcal{E}}$ and $\widehat{\mathcal{B}}$ are updated, and the averaged fields $\langle \widehat{\mathcal{E}}\rangle$ and $\langle \widehat{\mathcal{B}}\rangle$ are calculated and transformed to the spatial grid. As part of the particle push, the averaged fields $\langle E\rangle$ and $\langle B\rangle$ are gathered onto the macroparticles, to update the macroparticles' momenta. The rest of the PIC cycle (e.g., charge and current deposition, particle position push) is not shown here but is identical to the standard Galilean PSATD [16, 17]. See the appendices for more details on the PIC loop and exact definitions of the notation.

Figure 3

NCI growth rate: Galilean PSATD vs. averaged Galilean PSATD schemes. Normalized NCI growth rate $\text{Im}\left(\omega \right)/{\omega }_{p,r}$ in spectral ($k_x,k_z$) space, calculated from the analytical stability analysis (a), (b) and from WarpX simulation results (c), (d), obtained using the Galilean PSATD (a), (c) and averaged Galilean PSATD (b), (d) schemes, at infinite spectral order, with large timestep $c\mathrm{\Delta }t=\mathrm{\Delta }z=6\mathrm{\Delta }x$ and slightly detuned Galilean velocity $v_{\text{gal}}=0.99v_0$.

Figure 4

Instability mitigation in a 2D laser-wakefield simulation with large timestep. The upper half of each subplot shows the $E_z$ field seen by the macroparticles (i.e., the regular $E_z$ field for the standard Galilean PSATD, and the averaged $\langle E_z\rangle$ field for the averaged Galilean PSATD). The lower half of each subplot shows the longitudinal current density $J_z$ of the plasma. All the fields are shown in the boosted frame. The different subplots correspond to the standard Galilean PSATD (a), (c), (e) and averaged Galilean PSATD (b), (d), (f), with $c\mathrm{\Delta }t=0.1\mathrm{\Delta }z=\mathrm{\Delta }x$ (a), (b), $c\mathrm{\Delta }t=0.5\mathrm{\Delta }z=5\mathrm{\Delta }x$ (c), (d), and $c\mathrm{\Delta }t=\mathrm{\Delta }z=10\mathrm{\Delta }x$ (e), (f).

Figure 5

Instability mitigation in a 3D laser-wakefield simulation with large timestep. The upper half of each subplot shows the $E_z$ field seen by the macroparticles, i.e., the regular $E_z$ field for the standard Galilean PSATD (a), and the averaged $\langle E_z\rangle$ field for the averaged Galilean PSATD (b). The lower half of each subplot shows the longitudinal current density $J_z$ of the plasma. All the fields are shown in the boosted frame with $c\mathrm{\Delta }t=\mathrm{\Delta }z=10\mathrm{\Delta }x$.

Figure 6

Algorithm comparison for a 3D plasma-wakefield simulation. (a), (b) Snapshot of the $E_z$ field ($\langle E_z\rangle$ in the case of the averaged Galilean PSATD) in the laboratory frame, for the small-timestep $(c\mathrm{\Delta }t=\mathrm{\Delta }x)$ standard Galilean PSATD simulation (a) and large-timestep $(c\mathrm{\Delta }t=\mathrm{\Delta }z)$ averaged Galilean PSATD simulation (b) with fixed $\mathrm{\Delta }z=6.4\mathrm{\Delta }x$. The green dots are representative random samples of the macroparticles in the beam driver. (c), (d) Evolution of the emittance and relative energy spread of the beam, in the laboratory frame.