Pseudospectral Maxwell solvers for an accurate modeling of Doppler harmonic generation on plasma mirrors with particle-in-cell codes

With the advent of petawatt class lasers, the very large laser intensities attainable on target should enable the production of intense high-order Doppler harmonics from relativistic laser-plasma mirror interactions. At present, the modeling of these harmonics with particle-in-cell (PIC) codes is extremely challenging as it implies an accurate description of tens to hundreds of harmonic orders on a broad range of angles. In particular, we show here that due to the numerical dispersion of waves they induce in vacuum, standard finite difference time domain (FDTD) Maxwell solvers employed in most PIC codes can induce a spurious angular deviation of harmonic beams potentially degrading simulation results. This effect was extensively studied and a simple toy model based on the Snell-Descartes law was developed that allows us to finely predict the angular deviation of harmonics depending on the spatiotemporal resolution and the Maxwell solver used in the simulations. Our model demonstrates that the mitigation of this numerical artifact with FDTD solvers mandates very high spatiotemporal resolution preventing realistic three-dimensional (3D) simulations even on the largest computers available at the time of writing. We finally show that nondispersive pseudospectral analytical time domain solvers can considerably reduce the spatiotemporal resolution required to mitigate this spurious deviation and should enable in the near future 3D accurate modeling on supercomputers in a realistic time to solution.

Figure 1

Simulation configuration and effects of FDTD Maxwell solvers on harmonic spectra. For all simulations presented on this figure, an incident laser ($p$-polarized, normalized amplitude $a_0=8$, waist of $6{\lambda }_0$, Gaussian spatial profile, duration ${\tau }_0=10T_0$ with $T_0$ the laser period, hyper-Gaussian temporal profile) impinges with an angle ${\theta }_i={55}^{\circ }$ on the plasma mirror. The plasma mirror has an initial exponential density profile of gradient scale length $L={\lambda }_0/20$, where ${\lambda }_0$ is the laser wavelength, equal to $0.8\mu \text{m}$. The simulation box dimensions are $40{\lambda }_0\times 40{\lambda }_0$ in ($x,z$) directions. The FDTD order 2 Cole-Karkkainen-Cowan (CKC) Maxwell solver was used. (a) Spatial evolution of $B_y$ field and electron density at different time. In this simulation the spatial resolution was $\mathrm{\Delta }x=\mathrm{\Delta }z={\lambda }_0/70$ in ($x,z$) directions. 25 plasma particles per cells were used. The spatial laser amplitude profile ($B_y$ component) is sketched in red at a time $t=24.6T_0$ before reflection. The reflected field spatial amplitude profile ($B_y$ component filtered between harmonic orders 4 to 8) is represented in violet at a time $t=52.7T_0$ after reflection of the laser on the plasma mirror surface. The plasma mirror electron density spatial profile is represented in black at a time $t=39.13T_0$ during the laser-plasma mirror interaction. Panel (d) shows zoom of the plasma mirror surface. (b) Angularly resolved spectrum (log scale) of the reflected field computed from the same simulation as in (a), i.e., with a spatial resolution of $\mathrm{\Delta }x=\mathrm{\Delta }z={\lambda }_0/70$ in ($x,z$) directions. (c) Angularly resolved spectrum (log scale) of the reflected field computed from the same simulation as in (a) but with a higher spatial resolution of $\mathrm{\Delta }x=\mathrm{\Delta }z={\lambda }_0/140$ in ($x,z$) directions. In all cases, the time resolution used is fixed by the Courant condition of the CKC scheme (cf. Table 1).

Figure 2

Refractive index $n_r$ of vacuum for different Maxwell solvers as a function of normalized wave numbers $k\mathrm{\Delta }x$ and propagation angle $\theta =\text{arctan}(k_z/k_x)$ of the electromagnetic waves. For each solver, the time step is given by the Courant condition (cf. Table 1). (a) Index $n_r$ of vacuum for the Yee solver. (b) Index $n_r$ of vacuum for the CKC solver. (c) Index $n_r$ of the PSTD-order 128 solver. (d) Index $n_r$ of the PSATD-order 128.

Figure 3

Model for angular deviation of harmonics. The harmonic beam is modeled as an incident beam generated at the plasma-vacuum interface and entering vacuum with an incidence angle ${\theta }_i$. (a) The physical case is represented. Because the vacuum index is 1, there is no dispersion and the harmonics are generated at the angle of specular reflection. (b) The transmission medium is numerical, therefore dispersive. Each frequency is deflected with an angle ${\theta }_r$, which can be calculated with our simple model based on Snell's law.

Figure 4

Angularly resolved spectra (log scale) for different Maxwell's solvers in the range of ${\omega }_n/{\omega }_0=10$ to 35 for a spatial resolution of ${\lambda }_0/70$ and a time resolution fixed by the Courant condition. Panel (a): Yee solver. Panel (b): CKC solver. Panel (c): PSTD solver. Panel (d): PSATD solver. The laser incidence is ${\theta }_i={55}^{\circ }$. The white solid line is the curve expected by the refraction model. For (a) and (b), $n_r>1$ and the harmonics get closer to the normal at higher frequency with a positive angle of refraction (cf. Fig. 3). For (c), the PSTD solver numerical index $n_r<1$ and the refraction angle is negative. The panel (d) shows a dispersion-free spectrum obtained with the PSATD solver.

Figure 5

Angularly resolved spectrum (log scale) for PSTD [panel (a)] and PSATD [panel (b)]. The white solid line is the curve expected by the refraction model. For ${\omega }_n/{\omega }_0=17$, a numerical cutoff appears for PSTD due to the total reflection on the interface. This cutoff is not physical and is absent with PSATD in the same configuration.

Figure 6

Variations of the coefficients solution of the system (A4) in function of $\theta$, given incident angle. For a stable scheme, the full line curve must be below the dot line curve, i.e., $\eta \le {\eta }_{\mathrm{stable}}$. However, this condition is valid only for $\theta =\pi /4$ and corresponds to the Yee scheme.

Figure 7

Evolution of high harmonic generation efficiency as a function of spatial resolution in the case of the PSATD solver. Numerical and physical parameters of the simulation are the same as in the manuscript. Harmonic generation efficiency for different harmonic range is obtained by integrating the angularly resolved spectrum in angle. The $y$ axis is in log scale. Until at least the $30\text{th–}35\text{th}$ harmonic, this figure shows convergence of the high harmonic spectrum at a resolution of $\lambda /140$ in $2\mathrm{D}$.