A solver based on pseudo-spectral analytical time-domain method for the two-fluid plasma model

A number of physical processes in laser-plasma interaction can be described with the two-fluid plasma model. We report on a solver for the three-dimensional two-fluid plasma model equations. This solver is particularly suited for simulating the interaction between short laser pulses with plasmas. The fluid solver relies on two-step Lax–Wendroff split with a fourth-order Runge–Kutta scheme, and we use the Pseudo-Spectral Analytical Time-Domain (PSATD) method to solve Maxwell’s curl equations. Overall, this method is only based on finite difference schemes and fast Fourier transforms and does not require any grid staggering. The Pseudo-Spectral Analytical Time-Domain method removes the numerical dispersion for transverse electromagnetic wave propagation in the absence of current that is conventionally observed for other Maxwell solvers. The full algorithm is validated by conservation of energy and momentum when an electromagnetic pulse is launched onto a plasma ramp and by quantitative agreement with wave conversion of p-polarized electromagnetic wave onto a plasma ramp.

www.nature.com/scientificreports/ propagation in the absence of current (see Fig. 1 of Vay et al. 24 ). This method is also particularly well adapted for laser pulse propagation. The algorithm is tested with laser-plasma interaction problems with intensities around 10 14 to 10 15 W/cm 2 , since it is intended for the study of electron-hole plasma dynamics in solids. We model a plasma embedded in a medium of background relative permittivity ε r .
Here, we build a two-fluid plasma solver based on PseudoSpectral Analytical Time-Domain PSATD for solving Maxwell's equations. A schematic representation of our solver is given in Fig. 1. The integration of Maxwell's equations is performed by using the PSATD method. The electromagnetic fields are transmitted to the fluid equations as a Lorentz force source term. The fluid equations are integrated by using a Strang splitting 25 . In the splitting, the homogeneous system is solved via a Lax Wendroff (LW) scheme, while the source terms are integrated with a fourth-order Runge-Kutta scheme (RK4). The updated fluid variables are used to calculate the current density, which is injected in Maxwell's equations.
This paper is divided in four main parts. We will first recall the two-fluid plasma model equations before summarizing the numerical integration. We will then validate the solver and demonstrate its benefits in terms of numerical dispersion and in the reduction of the constraint imposed on time-steps with the solver of reference 21 .

Results
Two-fluid plasma equations. The two-fluid plasma model equations consists of Euler equations with source term for each fluid, as well as Maxwell's equations. This system of equations corresponds to continuity equations, motion equations and energy transport equations for electron and ion fluids. Following reference 26 , the fluid equations can be presented under the following form: where U is the fluid variables vector, F(U) is the flux tensor and S(U, E, B) is the Lorentz force source term. In this paper, i and e are indexes related respectively to the ion fluid and to the electron fluid. q is the charge, m the mass, ρ the mass density, u the mean velocity, p the pressure, ǫ the fluid energy density, E the electric field and B the magnetic field. ⊗ is tensor product and I is the identity matrix.
One more equation is required to close the system of equation, an ideal gas closure for each fluid k is used 26 : where γ is the adiabatic index. Electric and magnetic fields in the source term S(U, E, B) are determined by the Maxwell equations. It is known that solving fluid equations and Maxwell's curl equations together enforces the conservation of divergence properties of the fields 27 . Therefore, as long Maxwell's divergence equations are satisfied at the initial time, they continue to be satisfied during the whole simulation because of the combined numerical resolution of fluid equations and Maxwell's curl equations. Furthermore, Maxwell curl's equations can be written by expressing the current density J as function of fluid variables:  www.nature.com/scientificreports/ Here, ε 0 and µ 0 are respectively the vacuum permittivity and permeability. c = (ε 0 µ 0 ) −1/2 is the speed of light. ε r is the relative permittivity of the background medium: we do the assumption that this quantity is time and space independent. In this model, the plasma is also contained inside a medium of relative permittivity ε r .
The numerical integration. The Maxwell solver. As mentioned in the introduction, the PSATD method 24 is used for solving Maxwell curl's equations. This method is simple to implement and does not need the staggering of spatial and temporal grids. This is in contrast with the Finite Difference Time Domain (FDTD) method 28 which requires spatially and temporal staggered grids or in contrast with the PseudoSpectral Time-Domain (PSTD) 22 which requires a temporally staggered grid. The PSATD is therefore more flexible to be coupled with another algorithm without interpolations. Moreover, in absence of current, PSATD induces zero numerical dispersion in contrast with FDTD or PSTD. An additional strong benefit is that the PSATD is not subject to a Courant condition for transverse electromagnetic field propagation in the absence of current. The PSATD algorithm is inherently periodic because it is based on FFT, but open systems can be modeled by using Perfectly Matched Layers (PML) as in Shapoval et al. 29 .
The PSATD algorithm provides the fields in the Fourier space 24 : where ã is the Fourier transform of the quantity a . Here C 0 = cos (kv�t) , S 0 = sin (kv�t) , κ = k/k and v = c The two main assumptions made in the PSATD method are: (1) the time-step t is enough small to assume that the current density is constant over a time-step (2) the background permittivity ε r is uniform. We note that the plasma permittivity is not limited by this constraint because it is taken into account via charge currents.
The fluid solver. For the fluid equations solver, we consider a similar solver as in reference 21 . Here, we simplified the solver by restricting ourselves to problems without discontinuities such that it becomes unnecessary to introduce numerical dissipation to make gradient smoother. Instead of using a composite scheme LWLFn as in reference 21 , we will use a simple two-step Lax-Wendroff (LW) scheme 30 which is second order accurate and introduces less numerical dissipation than the two-step Lax-Friedrichs (LF) scheme 31 . The LW scheme solves the homogeneous part of Eq. (1), as we recall below. First, we set L x the operator for the two-step LW along x direction: with where j, l and m are respectively indexes for x, y and z directions. Similar operations are done in y and z directions as reference 21 to obtain L y (U n j,l,m ) and L z (U n j,l,m ) . A basic spatially dimensionally-split scheme is used to obtain the value U n+1 j,l,m from U n j,l,m 32 : For the numerical integration of the source term S(U, E, B) of Eq. (1), we use the Strang splitting technique presented by Strang 25 . The Strang splitting allows an estimation of current density J n+1/2 at a half time step of PSATD. The concept of Strang splitting is shown on the steps 1, 2 and 4 in Fig. 2. We first integrate the source term with an RK4 scheme over �t/2 , then the homogeneous system is integrated over t with an LW scheme, and finally source term is again integrated with an RK4 over time step �t/2.
Full two-fluid plasma solver algorithm. The full algorithm for the two-fluid plasma model is described in Fig. 2 and can be decomposed in 4 main steps: www.nature.com/scientificreports/ 1. Integration of the source term with an RK4 scheme over a temporal step �t/2 by using E n , B n and U n to obtain the intermediate value of fluid variables U * . 2. Integration of the homogeneous system with an LW scheme over a temporal step t using fluid variables vector U * to obtain a new intermediate value U * * . 3. Computation of the current density J n+1/2 with densities and velocities from U * * . Then, carry out a PSATD step with J n+1/2 to calculate E n+1 and B n+1 . 4. Integration of the source term with an RK4 algorithm over a temporal step �t/2 using U * * , E n+1 and B n+1 to obtain the final value of fluid variables U n+1 .
The PSATD naturally represents all field values at the nodes of a grid, it also avoids temporal interpolation of the magnetic field that was necessary in reference 21 .
For the PSATD algorithm alone without currents, the sampling is in principle only limited by Nyquist theorem. However, in order to derive Eqs. (5) and (6), we make the assumption that the current is constant over the temporal step t . Therefore, the temporal step is chosen small enough to make this assumption valid. The spatial step x is simply chosen to resolve the both plasma and electromagnetic waves.
Validation of the numerical solver. S-polarized electromagnetic wave over a plasma ramp. In this first test, we check the conservation of momentum and energy during reflection of a s-polarized electromagnetic wave over a plasma ramp. The numerical setup is shown in Fig. 3. A laser pulse is propagating toward an overcritical plasma ramp with an angle of incidence θ = 15 • . The initial plasma density profile is invariant in y and z directions, and the following initial density profiles, for electron and ion fluids, are used in x direction: where n e ≡ ρ e /m e and n i ≡ ρ i /m i are given in cm −3 . The length in x direction at which the critical density n c = 1.75 × 10 21 cm −3 is reached is L = 3.08 µm . We add a weak uniform background density of 10 17 cm −3 to avoid divisions by zero in the algorithm and too strong discontinuity at the ramp onset. In this test, the uniform background is vacuum: ε r = 1 . For the plasma, we take m e = 9.11 × 10 −31 kg , m i = 1837m e and γ = 5/3 . The initial mean velocities and pressure are zero.
The laser pulse is a spatially Gaussian beam with a waist w 0 = 4 µm and is described temporally by a single period of a sin 2 function (period T = 40 fs ). The free-space wavelength is = 0.8 µm and the amplitude in freespace is E 0 = 4.3 × 10 10 V/m . We choose this electric field amplitude to demonstrate the possibility of working with high field amplitudes with this algorithm. Note that the beam is invariant along z direction.
The number of points in x and y directions is N x = N y = 512 and N z = 2 is z direction. PML (resp. open) boundary conditions in x and y directions for PSATD (resp. fluid algorithm) are implemented. For fluids and fields, periodic boundary condition are used in z direction. The spatial step is x = y = z = 60 nm and the temporal step is = 96 as.
In Fig. 4a, we plot the different momenta in x direction as function of the simulation time. These momenta are normalized to the absolute value of the x momentum P 0 of the incident pulse. To measure the normalization factor P 0 , we performed beforehand the simulation without the plasma, and we measure the x momentum of the pulse defined by the x component of the Eq. (12) integrated over the simulation window.  www.nature.com/scientificreports/ The density of electromagnetic momentum is defined by 33 : The dashed red curve of Fig. 4a corresponds to the normalized electromagnetic x momentum density integrated in the simulation window. The dashed dotted blue curve corresponds to the normalized fluids x momentum integrated in the simulation window. The momentum of fluids is defined by: The black line of Fig. 4a is the sum of the electromagnetic and fluids x momentum. We observe three main sequences in Fig. 4a: • 1: Since the laser pulse goes from the right to the left, the electromagnetic momentum along x (red dashed curve) decreases as the pulse enters into the simulation window (between t = 0 and t < 50 fs ). At t = 50 fs , the pulse is completely contained in the simulation window and has not yet interacted with the plasma ramp. We see that the electromagnetic x momentum corresponds to the incident pulse x momentum −P 0 . • 2: In the temporal window 70-130 fs, momentum exchange with the plasma takes place: the fluids momentum decreases until −2P 0 , whereas the electromagnetic x momentum increases until reaches +P 0 . This is the signature that the laser pulse transfers twice its initial momentum to the plasma during its reflection, as can be expected.  We see that the momentum is preserved over the temporal window over which the pulse is fully enclosed within the simulation window. The error on the conservation of the total momentum is around 1%. It is reasonable in view of the chosen spatial and temporal steps. The numerical algorithm also preserves the conservation of momentum with a good accuracy. In Fig. 4b, we plot the linear density of energy as function of simulation time. The dashed red curve correspond to the electromagnetic energy density 33 that we have integrated over the x − y plane. The dashed dotted blue (resp. dashed green) corresponds to the electron fluid (resp. ion fluid) energy density given by Eq. (2) and then integrated over the x − y plane. The total energy plotted in black line is defined as the sum of electromagnetic, electron fluid and ion fluid energies. The linear density of energy of the input pulse in the x − y plane can be calculated analytically and is given by: E Laser =  Fig. 4b.
We observe the three main sequences in Fig. 4b: • 1: The electromagnetic energy increases as the pulse enters into the numerical window between t = 0 and t < 50 fs . At t = 50 fs , the pulse is completely contained in the numerical window and has not yet interacted with the plasma ramp. The electromagnetic energy corresponds to the predicted analytical value E Laser . • 2: In the temporal window 70-130 fs, energy exchange with the plasma takes place (the electron energy increases). • 3: Between 170 fs and 210 fs, the reflected pulse leaves the integration volume and the electromagnetic energy decreases. No electromagnetic energy remains in the simulation window. This is expected for s-polarized wave.
We remark the conservation of the energy when the pulse is fully in the simulation window. The error on the conservation the total energy is around 0.1%. The numerical algorithm also preserves the energy conservation with a good accuracy. Additionally, we show in Fig. 5 the error induced by the algorithm on the conservation of energy and momentum as a function of time step t . The calculations of errors are performed from simulations of s-polarized electromagnetic wave over a plasma ramp. In Fig. 5, we remark for �t < 150 as that the numerical error remains below 0.4% for energy conservation (blue curve), and below 2.8% for momentum conservation (red curve). Furthermore, we noted a loss of stability of PSATD/Hydro code for simulations if t exceeds 150 as with an input intensity of 2.5 × 10 14 W cm −2 .
Wave conversion on plasma density ramp. In this second test, we consider the same numerical setup as shown in Fig. 3, but we inject a p-polarized laser pulse.
The energy of the system is plotted as a function of time in Fig. 6a. In the central white area, the error on the conservation the total energy is around 0.1%. Furthermore, we observe for the sequence 4 that a fraction of the input energy remains in the simulation box while the laser pulse has left. This is due to the phenomenon of wave conversion onto a inhomogeneous plasma, i.e. conversion of an electromagnetic wave into a plasma wave which occur only for p-polarization 34 .
The conversion factor depends in particular on the plasma density gradient and the angle of incidence. Obtaining analytical solutions to this difficult problem usually requires a number of approximations. We www.nature.com/scientificreports/ performed a series of simulations with different angles of incidence, and we plot (blue circles) in Fig. 6b the factor of the energy conversion as function of the quantity τ 2 = 2πL 2/3 sin 2 θ. We compare conversion factors obtained with the PSATD/Hydrodynamic code (this work) and results of the literature. The PSATD/Hydro conversion factor curve is quantitatively superimposable to the one obtained the PSTD/hydro simulation (blue crosses) obtained in reference 21 . Our results are also in good agreement with the analytical results of Speziale et al. who described only the asymptotic behaviors 34 for τ → 0 and τ → ∞ , with the results of Hinkel-Lipsker et al. 35 obtained for any value of τ and also with those of Forslund et al. who have used Particle-In-Cell (PIC) simulations 36 . The fact that we injected a short pulse (polychromatic) gaussian beam instead a monochromatic plane wave can explain the tenuous differences. In addition, the analytical results of references Speziale et al. 34 and Hinkel-Lipsker et al. 35 , have carried out assumptions that are not exactly fulfilled in our numerical test, such as the linearization of fluid equations. But overall, the results we obtained with the present solver are in good agreement with the state of the art.

Discussion
In this section, we compare the benefits and drawbacks of PSATD/Hydro solver compared to PSTD/Hydro solver. We numerically simulate a single cycle pulse plane wave in normal incidence onto a plasma ramp. The laser wavelength and plasma parameters are identical to the ones of Fig. 3. The pulse amplitude is E 0 = 4.3 × 10 10 V/m . The computation is performed in 3D, with the same numerical sampling parameters as in Fig. 3. We use the periodic boundary conditions in y and z directions. For the PSTD/Hydro solver, the time-step is fixed to t = 50 as since we are constraint by Courant Friedrichs Lewy (CFL) conditions 22 . In contrast, the time-step for the PSATD/ Hydro solver is set to t = 200 as , as it is only constrained by the sampling of the laser and the plasma wave frequencies.
In Fig. 7a, we show a snapshot during propagation of the laser pulse with the different solvers. The snapshot is taken when the pulse has propagated through vacuum and just reaches the onset of the plasma ramp. We see that the PSATD/Hydro solver result (solid blue line) is precisely superimposed on the analytical solution in black dashed line. In contrast, the PSTD/Hydro solver (red dashed-dotted line) exhibits distortion of the laser pulse. Indeed, pre-pulses are generated by numerical dispersion of the PSTD algorithm. The amplitude of the last artifact pre-pulse (located at x ≈ −5 µm ) reaches around 15% of the amplitude of the main peaks.
In Fig. 7b, we plot the velocity component v z of the electron fluid at the same time of the snapshot of Fig. 7a. We observe that the laser pulse has not yet interacted with the plasma in the PSATD/Hydro simulation (blue line). However, in the PSTD/Hydro simulation (red dashed-dotted line), the artifact pre-pulses already interact with the plasma and accelerate the electrons to velocities around 10 5 m/s . This effect is obviously undesirable, particularly in the case of the simulation of few cycle laser pulses with plasmas 37 .
We also obtained better results in the PSATD/Hydro simulation whereas the time-step t was 4 times greater than those in PSTD/Hydro simulation.
The PSATD/Hydro solver is well suited to pulse propagation. Specifically, the fact that PSATD is not constrained by the CFL condition, releases the strong numerical link between spatial and temporal sampling. The computational gain is therefore particularly significant in the case where high spatial resolution is required together with less demanding temporal resolution. We finish this section by reminding that the PSATD method requires that the background medium permittivity is uniform.
As a conclusion, we have developed a solver for the two-fluid plasma model based on a relatively simple technique, which does not necessitate staggered grids and which benefits of fundamentally having no numerical We have shown that PSATD/Hydro solver has two main advantages compared to the PSTD/Hydro solver: the pre-pulses generated by numerical dispersion are removed and the time-step is not constraint by CFL conditions. For simulations which require low temporal resolution and high spatial resolution, the gain in terms of computational resources with PSATD/Hydro solver can be really significant. The PSATD/Hydro solver is a computationally inexpensive but powerful tool for the study of laser-plasma interaction. An interesting future development would be the extension of the PSATD algorithm to take into account the dispersion of background permittivity ǫ r with frequency.

Methods
Simulations were performed with Nvidia Tesla K40 GPU card. This card has 12 GB memory size, 2880 CUDA cores and 745 MHz processor core clock. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.