A simple solver for the two-fluid plasma model based on PseudoSpectral Time-Domain algorithm
AA simple solver for the two-fluid plasma model based on PseudoSpectralTime-Domain algorithm
B. Morel, R. Giust, K. Ardaneh and F. Courvoisier ∗ Institut FEMTO-ST, Universite Bourgogne Franche-Comte,15B Avenue des Montboucons, 25030 Besancon cedex, France
We present a solver of 3D two-fluid plasma model for the simulation of short-pulse laser inter-actions with plasma. This solver resolves the equations of the two-fluid plasma model with idealgas closure. We also include the Bhatnagar-Gross-Krook collision model. Our solver is based onPseudoSpectral Time-Domain (PSTD) method to solve Maxwell’s curl equations. We use a Strangsplitting to integrate Euler equations with source term: while Euler equations are solved with acomposite scheme mixing Lax-Friedrichs and Lax-Wendroff schemes, the source term is integratedwith a fourth-order Runge-Kutta scheme. This two-fluid plasma model solver is simple to implementbecause it only relies on finite difference schemes and Fast Fourier Transforms. It does not requirespatially staggered grids. The solver was tested against several well-known problems of plasmaphysics. Numerical simulations gave results in excellent agreement with analytical solutions or withprevious results from the literature. a r X i v : . [ phy s i c s . p l a s m - ph ] O c t LWLFn
Fluid equations
PSTD RK4
Maxwell equations
Strang splitting
FIG. 1. General overview of the algorithm.
INTRODUCTION
Since the early 1960s, the interaction of laser light with a plasma has revealed to be an extremely rich topic :laser-plasma accelerators, inertial fusion, X-Ray sources, nonlinear plasmonics or laser materials processing [1]. Anumber of plasma effects have been characterized in laser-plasma experiments, but many challenging problems remain[2]. Simulations in laser-plasma interaction domain are essential for the understanding of the physical phenomena inhigh-intensity laser interactions. However, this requires efficient and precise numerical simulations. The hydrodynamicapproach is particularly interesting to save computational effort when each of the species can be assumed in localthermodynamic equilibrium [3].The most fundamental approach for hydrodynamic code is the set of multi-fluid plasma equations [4]. Theseequations can be derived by taking moments of the Vlasov equation with respect to velocity and truncating momentequations by making additional assumptions. Here, we will consider a two-fluid system of electrons and ions. Theclosure of the system of equations is performed via ideal gas closure. This model describes the evolution in spaceand time of the density, mean velocity and pressure of electrons and ions fluids. The two-fluid plasma equationstherefore consist of two sets of Euler equations with source term, as well as Maxwell’s equations. More conventionalplasma models, e.g., two-temperature plasma equations, single-fluid equations and MagnetoHydroDynamic (MHD)are derived from the two-fluid plasma model but require additional assumptions.Until now, only few numerical codes solve the complete two-fluid plasma equations [5]. However, their implemen-tation is often complex for non-specialist groups. ANTHEM code was the first two-fluid plasma algorithm [6, 7]. Itis based on a flux-corrected transport (FCT) algorithm for fluids equations and a Finite-Difference Time Domain(FDTD) algorithm for the evolution of the fields. This code was used to simulate high-density plasmas out of reachto Particle-In-Cell (PIC) codes. It was although difficult to extend to arbitrary geometries because of the staggeredgrid of the FDTD scheme in the Yee algorithm [8]. More recently, U. Shumlak et al. [9] presented an algorithmbased on Roe-type Riemann solver [10] for the two-fluid plasma model. Time updates are carried out by treating thesource term implicitly and the flux terms explicitly. This code was developed to observe two-fluid effects which areinaccessible to magnetohydrodynamic models. The same group added the high-order discontinuous Galerkin methodto improve the result’s accuracy [11–14]. A. Alvarez Laguna et al. [15] introduce another algorithm where the two-fluid plasma equations are solved with a fully-implicit finite volume method for unstructured meshes combined witha modified version of Rusanov scheme [16] for Maxwell equations. However, these solvers are relying on complexmathematical methods requiring specialized skills. Other approaches have been developed relying on finite differencealgorithm: S. Baboolal [17] reports a one-dimensional scheme for shock and soliton simulations. H. Kumar et al. [18]describe a scheme combining entropy conservative fluxes and suitable numerical diffusion operators. However, thefirst is limited to one dimension and the second requires heavy finite difference scheme to avoid entropy dissipationfor shocks.These codes are usually developed for complex problems as nuclear fusion and astrophysics, and also tested onvery complex problems as the Brio and Wu shock tube test [19]. Our code is tested with laser/plasma interactionproblems with intensities around 10 to 10 W/cm , since this code is intended for the study of electron/hole plasmadynamics in solids.In this framework, we present a new simple method to solve equations of the adiabatic inviscid five-moment two-fluid plasma model with collisions. It only relies on finite-difference schemes and Fast Fourier Transforms (FFT)over a spatially non-staggered grid. This solver is intended to be used in absence of discontinuities. Our approach issummarised in Fig. 1. Maxwell’s equations are solved using PseudoSpectral Time-Domain (PSTD) algorithm. Theelectromagnetic fields are passed to the fluid equations as source terms. We solve the fluid equations using a Strangsplitting. In the splitting, the Euler equations (without source terms) are solved with a composite scheme, mixingLax-Friedrichs and Lax-Wendroff schemes, while the source terms are integrated via fourth order Runge-Kutta scheme(RK4). The computed currents are then used as sources for Maxwell’s equations.This method has successfully been tested against classical plasma physics problems covering: i/ the three wavepropagation modes in warm, fully ionized, isotropic plasmas; ii/ mode conversion and iii/ a nonlinear laser-plasmainteraction regime.This paper is organised as follows. In sections 2, 3, 4 and 5 we briefly recall the models and methods on whichour approach is based: the two-fluid plasma model, the collision model, the PseudoSpectral Time-Domain (PSTD)algorithm and the fluid equation solver. This will provide the reader a consistent way of writing the equations. Then,in section 6 we describe our algorithm and the stability conditions. In section 7, we present the 5 different testsvalidating the solver. TWO-FLUID PLASMA EQUATIONS
As for any hydrodynamic code, we start from the two-fluid plasma equations [4]. We consider here a plasmacomposed of an electron fluid and an ion fluid where we assume no heat flow and no viscosity. The inviscid five-moment two-fluid plasma model equations can be derived from the Vlasov-Maxwell equations [20]. The two-fluidplasma model equations listed below consists of Euler equations with source term for each fluid, as well as Maxwell’sequations. We first write the collisionless equations and we will see in the next section how they can be introducedas an additional source term. Following Ref. [20], the fluid equations are written as balance laws, i.e. a conservativeform on the left side and a source term on the right one. This approach offers access to a large number of numericalmethods. ∂∂t U + ∇ · [ F ( U )] = S ( U , E , B ) (1)where U is the fluid variables vector: U ≡ ρ e ρ e u e (cid:15) e ρ i ρ i u i (cid:15) i (2) F ( U ) is the flux tensor: F ( U ) ≡ ρ i u e ρ e u e ⊗ u e + p e I ( (cid:15) e + p e ) u e ρ i u i ρ i u i ⊗ u i + p i I ( (cid:15) i + p i ) u i (3)and S ( U , E , B ) is the source term: S ( U , E , B ) ≡ ρ e q e m e ( E + u e × B ) ρ e q e m e u e · E ρ i q i m i ( E + u i × B ) ρ i q i m i u i · E (4)Here, e and i are indexes related respectively to the electron fluid and to the ion fluid. I is the identity matrix and ⊗ is tensor product. q k is the k specie charge, m k the specie mass, ρ k the mass density, u k the mean velocity, p k the pressure, (cid:15) k the fluid energy density, E the electric field and B the magnetic field. The system of equations (1)corresponds to continuity equations, motion equations and energy transport equations for electron and ion fluids. Thesystem of equations is closed with ideal gas closure for each fluid k [20] : (cid:15) k ≡ p k γ − (cid:124) (cid:123)(cid:122) (cid:125) Thermal energy + 12 ρ k || u k || (cid:124) (cid:123)(cid:122) (cid:125) Kinetic energy (5)where γ is the adiabatic index.The electric and magnetic fields appearing in the source term of the Euler equations are determined by the Maxwellequations: ∇ · (cid:15) r E = 1 (cid:15) (cid:20) q e m e ρ e + q i m i ρ i (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) Total charge density (6) ∇ · B = 0 (7) ∇× E = − ∂ B ∂t (8) ∇× B = µ (cid:20) q e m e ρ e u e + q i m i ρ i u i (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) Current density J + 1 c ∂(cid:15) r E ∂t (9)where (cid:15) and µ are respectively the vacuum permittivity and permeability. (cid:15) r is the relative permittivity of themedium, and c = ( (cid:15) µ ) − / is the speed of light. COLLISIONS IN THE TWO-FLUID PLASMA MODEL
The addition of a basic collisional model is straightforward in the two-fluid plasma model. One of the simplestmodel is the BGK (Bhatnagar-Gross-Krook) model [21]. In this model, it is assumed that the effect of collisions isto restore a situation of local equilibrium. It assumes that the local equilibrium is reached exponentially with time.The inclusion of the BGK collision model in the two-fluid plasma equations is handled through an additional sourceterm given by [22] : S coll ( U ) ≡ − ρ e m e µν ei ( u e − u i ) − κν ei (cid:104) p e γ − − ρ e m i ρ i m e p i γ − + ρ e m e (cid:0) m e u e2 − m i u i2 + ( m i − m e ) u e · u i (cid:1)(cid:105) − ρ e m e µν ei ( u i − u e )+ κν ei (cid:104) p e γ − − ρ e m i ρ i m e p i γ − + ρ e m e (cid:0) m e u e2 − m i u i2 + ( m i − m e ) u e · u i (cid:1)(cid:105) (10)where µ = m e m i m e + m i is the reduced mass, κ = µm e + m i the energy transfer coefficient, and ν ei the collision frequency formomentum transfer between electrons and ions. This collision term preserves momentum and energy conservation. THE MAXWELL SOLVER
Numerical methods to solve Maxwell equations have been widely studied [8, 23–27]. The PSTD algorithm [23] is apseudo-spectral method which is simple to implement and does not need a spatially staggered grid, in contrast withthe Yee FDTD algorithm [8]. The PSTD is therefore more flexible to be coupled with another algorithm without toomany interpolations. Importantly, PSTD is low-dispersive since it provides accurate dispersion relation with only twocells per wavelength. The relative error for the PSTD algorithm with 2 cells per wavelength is smaller than that forthe FDTD method with 16 cells per wavelength [23]. The PSTD algorithm imposes periodicity because it is based onthe Fast Fourier Transform (FFT), but this constraint can be removed by using Berenger’s perfectly matched layers(PML) [28]. Furthermore, the PSTD scheme allows to implement non-uniform background permittivity (cid:15) r . It is alsopossible to model non-linear polarization term as shown in T.-W. Lee et al. [29]. Therefore, PSTD algorithm is anexcellent candidate to construct a simple and reliable hydrodynamic solver.We note that solving fluid equations and Maxwell’s curl equations enforces the conservation of divergence propertiesof the fields [30]. Hence, if fields divergence properties are satisfied at an initial time, they will be satisfied at a latertime. It is therefore not necessary to solve equations (6) and (7).In PSTD, spatial derivatives are approximated using an FFT algorithm [23]. Moreover, central finite difference isused for temporal derivatives. Thus, the spatial grid is unstaggered but the temporal grid is staggered. For example,if E is defined at t = n ∆ t , B is defined at t = (cid:0) n + (cid:1) ∆ t , where ∆ t is the time step of the PSTD algorithm.We present here 3D PSTD algorithm valid for a uniform medium, to which current density is added. The spatialderivatives are performed in Fourier space. Eq. (11) shows the partial derivative of E y with respect to x : (cid:18) ∂E y ∂x (cid:19) n = iFFT x (cid:2) ik x FFT x (cid:2) E ny (cid:3)(cid:3) (11)where iFFT x denotes to the inverse FFT along the x -axis, and k x is the spatial frequency along the x -axis. All thespatial derivatives will be done in Fourier space as shown in Eq. (11).Knowing the magnetic field at t n − / and spatial derivatives of electric field at t n , magnetic field B n+1 / can becalculated: B n+1 / = B n − / + ∆ t (cid:20)(cid:18) ∂E y ∂z (cid:19) n − (cid:18) ∂E z ∂y (cid:19) n (cid:21) (12) B n+1 / = B n − / + ∆ t (cid:20)(cid:18) ∂E z ∂x (cid:19) n − (cid:18) ∂E x ∂z (cid:19) n (cid:21) (13) B n+1 / = B n − / + ∆ t (cid:20)(cid:18) ∂E x ∂y (cid:19) n − (cid:18) ∂E y ∂x (cid:19) n (cid:21) (14)The spatial derivatives of B n+1 / are carried out as the same way than for electric field (see Eq. (11)). Finally, E n+1 is updated using E n and spatial derivatives of the magnetic field: E n+1x = E nx + c ∆ t (cid:34)(cid:18) ∂B z ∂y (cid:19) n+1 / − (cid:18) ∂B y ∂z (cid:19) n+1 / (cid:35) − ∆ t(cid:15) J x (15) E n+1y = E ny + c ∆ t (cid:34)(cid:18) ∂B x ∂z (cid:19) n+1 / − (cid:18) ∂B z ∂x (cid:19) n+1 / (cid:35) − ∆ t(cid:15) J y (16) E n+1z = E nz + c ∆ t (cid:34)(cid:18) ∂B y ∂x (cid:19) n+1 / − (cid:18) ∂B x ∂y (cid:19) n+1 / (cid:35) − ∆ t(cid:15) J z (17)where J x , J y and J z are the current density. Current density will be provided later on from fluid variables.We finish this section by recalling the PSTD stability criterion in D dimensions [23] :∆ t ≤ π ∆ xc √ D (18)where ∆ x is the spatial-step. THE FLUID SOLVER
In this section, we first discuss the numerical method chosen for homogeneous 3D Euler equations. We thenintroduce the method used to integrate the source term.
Numerical resolution of Euler equations
In the vectorial Eq. (1), fluid equations are written under conservative form. The conservative form of Eulerequations emphasizes that mass density, momentum and energy are conserved. Computationally, there are someadvantages in using the conservative form and consequently it gives access to a large class of numerical methodscalled conservative methods [31]. Here, the goal is to numerically solve Eq. (1) with a simple and reliable scheme.As mentioned in the introduction, extremely reliable schemes have been developed in the past, that can even handlelarge discontinuities and shocks, but are excessively complex for problems without such strong requirements. Veryattractive schemes are the composite ones introduced by Liska and Wendroff [32]. Theses schemes are simple andhigh-resolution [33]. They consist of a temporal composition of basic classical finite difference schemes. Here we havechosen a composite scheme based on two-step Lax-Wendroff and two-step Lax-Friedrichs schemes.
The LWLFn composite scheme
The Euler equations are numerically solved using an LWLFn composite scheme [32]. It consists of applying ( n − n can be tuned depending on the problem. If the solution contains strong gradients, it is betterto use a small number n (typically LWLF4). On the other hand, when the solution is smooth, it is possible to useexclusively the two-step LW scheme. When shape of solution is initially completely unknown, it is advisable to beginnumerical simulation with a small number n and to increase it progressively. The LWLFn scheme offers a wider fieldof applications than LW or LF scheme alone, without additional complexity. Two-step Lax-Friedrichs scheme
The classical (single-step) LF scheme is a basic method to solve Hyperbolic Partial Differential Equations (HPDE)[36]. The two-step LF is a variant of the classical LF scheme which allows better damping of high-frequencies, andless damping of low and middle frequencies than the classical method [37]. It is a better method to filter high-frequency oscillations generate by the two-step LW scheme. This scheme is stable under CFL (Courant-Friedrichs-Lewy) conditions [37] but, its use alone is often limited because it is only first-order accurate.The two-step LF can be extended to multidimensional in a simple way with the symmetrized dimensionally-splitscheme (SYS) technique [35]. The dimensionally-split schemes are created by successively approximating solutionsof the 3D conservation laws in solutions of three 1D problems: ∂∂t U + ∂∂x F x ( U ) = , ∂∂t U + ∂∂y F y ( U ) = and ∂∂t U + ∂∂z F z ( U ) = .First, we define L x the operator for the two-step LF step along x : L x ( U nj , l , m ) = 12 (cid:104) U n+1 / / , l , m + U n+1 / − / , l , m (cid:105) − ∆ t ∆ x [ F x n+1 / / , l , m − F x n+1 / − / , l , m ] (19)with U n+1 / / , l , m = 12 (cid:2) U nj+1 , l , m + U nj , l , m (cid:3) − ∆ t x (cid:2) F x nj+1 , l , m − F x nj , l , m (cid:3) (20) U n+1 / − / , l , m = 12 (cid:2) U nj , l , m + U nj − , l , m (cid:3) − ∆ t x (cid:2) F x nj , l , m − F x nj − , l , m (cid:3) (21)where j , l and m are respectively indexes for x , y and z positions. To make the notation less cluttered, we have set F x nj , l , m ≡ F x (cid:16) U nj , l , m (cid:17) .Similarly, we define L y the operator for the two-step LF step along y : L y ( U nj , l , m ) = 12 (cid:104) U n+1 / , l+1 / , m + U n+1 / , l − / , m (cid:105) − ∆ t ∆ y [ F y n+1 / , l+1 / , m − F y n+1 / , l − / , m ] (22)with U n+1 / , l+1 / , m = 12 (cid:2) U nj , l+1 , m + U nj , l , m (cid:3) − ∆ t y (cid:104) F y nj , l+1 , m − F y nj , l , m (cid:105) (23) U n+1 / , l − / , m = 12 (cid:2) U nj , l , m + U nj , l − , m (cid:3) − ∆ t y (cid:104) F y nj , l , m − F y nj , l − , m (cid:105) (24)Similarly, we define L z the operator for the two-step LF step along z : L z ( U nj , l , m ) = 12 (cid:104) U n+1 / , l , m+1 / + U n+1 / , l , m − / (cid:105) − ∆ t ∆ z [ F z n+1 / , l , m+1 / − F z n+1 / , l , m − / ] (25)with U n+1 / , l , m+1 / = 12 (cid:2) U nj , l , m+1 + U nj , l , m (cid:3) − ∆ t z (cid:2) F z nj , l , m+1 − F z nj , l , m (cid:3) (26) U n+1 / , l , m − / = 12 (cid:2) U nj , l , m + U nj , l , m − (cid:3) − ∆ t z (cid:2) F z nj , l , m − F z nj , l , m − (cid:3) (27)The symmetrized dimensionally-split scheme (SYS) allows to obtain the value U n+1j , l , m from U nj , l , m by averaging allpossible simple dimensionally-split schemes [35]: U n+1j , l , m = 16 ( L x L y L z + L x L z L y + L y L x L z + L y L z L x + L z L x L y + L z L y L x ) U nj , l , m (28) Two-step Lax-Wendroff scheme
The LW scheme is a classical numerical method for HPDE, based on finite differences [38]. But for non-linearHPDE, the classical Lax-Wendroff scheme requires evaluating Jacobian matrices. Richtmyer [39] proposed a two-stepLW scheme in order to avoid calculating Jacobian matrices. The stability of this scheme is constrained by CFLconditions [39].The 3D two-step LW is obtained just by replacing equations (19), (22) and (25) respectively by : L x ( U nj , l , m ) = U nj , l , m − ∆ t ∆ x (cid:104) F x n+1 / / , l , m − F x n+1 / − / , l , m (cid:105) (29) L y ( U nj , l , m ) = U nj , l , m − ∆ t ∆ y (cid:104) F y n+1 / , l+1 / , m − F y n+1 / , l − / , m (cid:105) (30)and L z ( U nj , l , m ) = U nj , l , m − ∆ t ∆ z (cid:104) F z n+1 / , l , m+1 / − F z n+1 / , l , m − / (cid:105) (31) t n Source termHomogeneoussystemMaxwellequations U n E n B n U ** E n+1/2 B n+1/2 E n U n+1 E n+1/2 B n+1/2 U * E n B n U n+1 E n+1/2 B n+1/2 U ** E n B n U n+1 E n+1 B n+1 t n+1/2 t n+1 t RK4 Δ t/2PSTD Δ t/2 RK4 Δ t/2PSTD Δ t/2LWLFn Δ t Initial state U * E n B n U ** B n U ** E n+1/2 B n+1/2 Final stateIntermediate state
FIG. 2. Description of the algorithm to solve the two-fluid plasma equations. An RK4 scheme is used to integrate the sourceterm, an LWLFn scheme is used to integrate the homogeneous system and a PSTD algorithm is used to advance the fields.The algorithm involves five steps shown in blue to advance fields and fluid variables from a time step n to n + 1. Source term integration
In this subsection, we select a numerical technique to integrate the source term in Eq. (1), which is a seriousnumerical challenge according to [31]. The most widely used method is the splitting method which splits the systeminto two sub-problems, homogeneous system and the source system. Each of these systems are solved separately.Here, we employ the Strang splitting technique as presented in [40] and [20]. This splitting is second-order accuratewhereas a classical Godunov splitting is only first-order accurate [36]. In the Strang splitting, the system (1) is splitas follows: ∂∂t U + ∇ · [ F ( U )] = (32) ∂∂t U = S ( U , E , B ) (33)Eq. (32) is a HPDE system while Eq. (33) is a Ordinary Differential Equation (ODE) system. The HPDE system issolved using LWLFn scheme and a fourth-order Runge-Kutta (RK4) scheme [41] is used to integrate the ODE system.In the Strang splitting, we first solved Eq. (33) with an RK4 scheme over ∆ t/
2, then the homogeneous Eq. (32) isintegrated over ∆ t with an LWLFn scheme, and finally Eq. (33) is again solved with an RK4 over time step ∆ t/ SOLVER STRUCTURE
As presented in previous sections, we solve Maxwell’s curl equations with a PSTD algorithm and the combinationof LWLFn and RK4 solves the fluid equations based on a Strang splitting. The aim of this section is to describe thecoupling of these two algorithms. This coupling is performed via the current density for Maxwell’s equations and viathe electromagnetic fields for the the fluid equations.The full algorithm is detailed in Fig. 2. This figure shows how to advance the values of fluid variables U , electricfield E and magnetic field B from time t n to t n+1 .The initial state is defined as: • fluid variables vector U n is known at time t n . • The electric field is known at time t n . Due to the staggered temporal grid of PSTD, the magnetic field B isknown at time t n − / and at t n+1 / . We deduce B n with a linear interpolation between B n − / and B n+1 / .Then algorithm follows five main steps: • Step 1: Integrate the source term (33) with an RK4 scheme over a temporal step ∆ t/ E n , B n and U n to obtain the intermediate value of fluid variables U ∗ . • Step 2: Integrate the homogeneous system (32) with an LWLFn scheme over a temporal step ∆ t using fluidvariables vector U ∗ to obtain a new intermediate value U ∗∗ . • Step 3: Calculate the current density J with densities and velocities from U ∗∗ . Then, carry out a PSTD stepwith J to calculate E n+1 / and B n+3 / . Calculate B n+1 / using a linear interpolation between B n+1 / and B n+3 / . • Step 4: Integrate the source term (33) with an RK4 algorithm over a temporal step ∆ t/ U ∗∗ , E n+1 / and B n+1 / to obtain the final value of fluid variables U n+1 . • Step 5: Calculate the current density J with densities and velocities from U n+1 . Then, carry out a PSTD stepwith J to calculate E n+1 and B n+5 / . Finally, calculate B n+1 using a linear interpolation between B n+3 / and B n+5 / .In order to capture the physics of the system, it is necessary to resolve the maximal frequency of the system [12],which is either the laser pulsation ω or the hybrid pulsation [20]: ω max = (cid:113) ω + ω (34)where ω pe = (cid:112) ρ e q /m (cid:15) is the electron plasma pulsation and ω ce = q e B/m e is the cyclotron plasma pulsation.In practice, when transverse electromagnetic waves are expected or injected, we use the stability criterion (Eq.(18)) to ensure the stability of the PSTD algorithm, which is usually the most demanding constraint. In absence oftransverse electromagnetic waves, the sampling is adapted to resolve both the plasma thermal velocity and the hybridpulsation. VALIDATION OF THE SOLVER
The implementation of the LWLFn algorithm was first validated against the Noh test problem in 2D and 3D asin Refs. [35, 42], which is known as a demanding test. We obtained quantitatively identical results. We note thatfor simpler tests such as the smooth problem of [35], the symmetrization of Eq. (28) is not absolutely necessary. Abasic spatial splitting such as U n+1j , l , m = L x L y L z ( U nj , l , m ) is enough to obtain reliable results while significantly savingcomputation time.Now, we validate our full solver against several classical plasma physics problems in which analytical solutions areavailable. In a warm fully ionized isotropic plasma, there are three modes of wave propagation [30]. The goal ofthe first two tests is to retrieve with the solver the dispersion relations of these three modes. In the third test, wecheck mode conversion between an electromagnetic mode and a longitudinal electron plasma mode. The fourth testvalidates the integration of the collisional model. Finally in the fifth test, we demonstrate that our solver is alsovalid in the nonlinear regime such as wakefield generation driven by ponderomotive force. We remark that in allthe following tests, we have used the general framework of the LWLFn with n=5000 to keep the generality of ourapproach. Although in our case the LW scheme could be sufficient, problems with higher gradients may require toincrease the number of Lax-Friedrichs steps in order to obtain convergence of the numerical solution (for example inthe 3D Noh test of reference [35], as we used to validate the fluid solver of our code). For instance, the dispersioncurves shown in Fig. 4 are identical whether the pure LW or LWLF5000 are used.The simulation setups and the boundary conditions (BC) for the five tests are summarized in table I. For all thesetests, we set m e = 9 . × − kg, q e = − . × − C, q i = − q e , (cid:15) = 8 . × − F/m, (cid:15) r = 1, c = 3 . × m/sand the Boltzmann constant is k B = 1 . × − J/K.0
TABLE I. Summary of different variables for the five tests.Parameter Test 1 Test 2 Test 3 Test 4 Test 5∆ x = ∆ y = ∆ z
10 nm 0 . t
14 as 6 .
34 as 96 as 28 as 14 asLWLFn LWLF5000 LWLF5000 LWLF5000 LWLF5000 LWLF5000 N x N y N z x PSTD BC Periodic Periodic PML Periodic Periodic y PSTD BC Periodic Periodic Periodic Periodic Periodic z PSTD BC PML PML PML PML PML x LWLFn BC Periodic Periodic Open Periodic Periodic y LWLFn BC Periodic Periodic Periodic Periodic Periodic z LWLFn BC Open Open Open Open OpenInitial u e and u i m i /m e p e and p i ν ei (fs − ) 0 0 0 variable 0.5 γ Test 1: Transverse electromagnetic wave propagation in an unmagnetized plasma
The dispersion relation for transverse electromagnetic plasma mode in unmagnetized plasma is given by [30]: ω = ω + ω + (cid:18) πλ p (cid:19) c (35)where ω pe and ω pi are respectively the electron and ion plasma frequencies. Thus the wavelength of the electromagneticwave in the plasma is: λ p = 2 πc (cid:113) ω − ω − ω (36)For numerical simulations, we consider the following initial densities and fields: • The initial densities are uniform. Several simulations are performed for different densities ranging between0 . × cm − and 1 . × cm − . • The electric and magnetic fields are initially set to zero everywhere. A monochromatic plane wave propagatingin positive z -direction is injected with a transient on ∼
50 fs. The wave is polarized along the x -direction witha wavelength λ = 0 . µ m and an amplitude in vacuum E = 10 V/m.The final time of the simulations is t final = 150 fs. For each simulation, a FFT z is performed on E x field to extractthe wavelength in the plasma. The wavelengths from the numerical simulations are plotted as a function of plasmadensity with blue points in Fig. 3. The error bar is defined by the spatial frequency sampling. The numerical resultsare in excellent agreement with the theoretical curve from Eq. (36). Test 2: Longitudinal modes in a warm plasma
Isotropic plasmas support two longitudinal plasma waves: longitudinal electron plasma waves and longitudinal ionplasma waves. Linearization of the two-fluid plasma equations for small amplitude perturbations gives the followingdispersion relation for longitudinal electron plasma mode [30]: ω = ω + ω + γ k B T e m e k (37)1 Electron density (cm -3 ) p ( µ m ) Numerical resultsTheoretical curve
FIG. 3. Wavelength of the transverse electromagnetic waves in the plasma as function of electron density. Numerical resultsare shown in blue points while the red curve is plotted from the theoretical relation 36. where T e is the electron temperature. The dispersion relation of the longitudinal ion plasma wave in low-frequencylimit (cid:16) ω (cid:28) ω (cid:104) T i T e (cid:105)(cid:17) reads [30]: ω = γ k B ( T i + T e ) m i k (38)where T i is the ion temperature. In the high-frequency limit (cid:16) ω (cid:29) ω (cid:104) T i T e (cid:105)(cid:17) , the dispersion relation of thelongitudinal ion plasma wave is given by [30]: ω = γ k B T i m i k (39)The goal of this second test is to reproduce relations (37), (38) and (39). For this purpose, we consider the followinginitial conditions: • The initial ion density is uniform: ρ i /m i = 5 × cm − . • The initial electron density is uniform with a small perturbation in the vicinity of z = 0 : ρ e = ρ e0 (cid:104) d zW z exp (cid:0) − z /W (cid:1)(cid:105) where ρ e /m e = 5 × cm − , ∆ d = 10 − and W z = 1 nm. • Based on Gauss’s law, the E z component associated to this perturbation is E z = − ρ e0 q e W z ∆ d m e (cid:15) exp (cid:0) − z /W (cid:1) .The others components of electric and magnetic fields are initially set to zero. • The initial electron and ion temperatures are: T e = 5 × K and T i = 2 × K respectively. The initialpressures are linked to temperature via ideal gas law ( p = ρm k B T ). • The ion-to-electron mass ratio m i /m e is fixed to 25 so as to observe both electron and ion waves in the samesimulation.We plot in log scale in Fig. 4 the dispersion diagram ( ω , k z ) of the E z field after 130 fs of simulation. The dispersionrelation of the longitudinal electron plasma wave (Eq. 37), and dispersion relation of the longitudinal ion plasmawaves in two asymptotic limits (Eqs. (38) and (39)), are shown as black dashed lines. The numerical results show avery good agreement with the analytical ones. We remark that in this case, since only small gradients are involved,the first order LF scheme alone, the second order LW scheme alone or the composite scheme LWLFn provide thesame result.2 k z (µm -1 ) ( f s - ) Longitudinal electron plasma waveLongitudinal ion plasma wave
Low-frequency asymptote High-frequency asymptote
FIG. 4. Dispersion diagram ( ω , k z ) of longitudinal waves in a warm plasma. Black dashed lines correspond to the theoreticaldispersion relations of longitudinal electron waves and longitudinal ion waves in two asymptotic limits (low and high frequency).The colorbar is in logarithmic scale. -3 U em U Kinetic U Thermal E Pulse (a) (b) -10 -5 0 5 10 z (µm) t (fs) n / n c r E ne r g y den s i t y ( J / m )
500 100 150 200 -10-50510 x ( µ m ) T u r n i ng s u r f a c e C r i t i c a l s u r f a c e k in k out V FIG. 5. (a) Concept of the simulation. A p-polarized laser field impinges on an inhomogeneous plasma slab (see main text).The blue line shows the plasma density profile. The white dashed rectangle shows the contour of the integration volume V .(b) Temporal evolution of the linear densities of: electromagnetic energy (dashed blue line), kinetic energy (dashed-dotted redline), thermal energy (dotted green line), and total energy (solid purple line) integrated over the volume V for the laser pulsedescribed in (a) with an incidence angle of 15 ◦ (see main text). The initial pulse linear density of energy E Pulse is plotted as adotted black line.
Test 3: Wave conversion on plasma density ramp
Here we test wave conversion in a cold, unmagnetized, collisionless plasma. Wave conversion occurs when anelectromagnetic wave, p-polarized, is obliquely incident onto a inhomogeneous plasma [43]. At the critical surface,the local plasma frequency equals the frequency of the electromagnetic wave and since the incident electric field hasa component along the density gradient, a plasma wave is driven resonantly. Therefore, part of the energy of theincident light wave is transferred to the electron plasma wave. Analytical solutions to this difficult problem usuallyrequire a number of approximations. This field of research has been extensively investigated and we have selectedout for comparison, the results from the literature considered as the most accurate. Speziale et al. described theasymptotic behaviors [43], Forslund et al. obtained numerical results based on PIC simulations [44], and the analyticalresults of Hinkel-Lipsker et al. match well with numerical simulations and experimental results [2, 45].Our simulation setup is shown in Fig. 5(a). A spatially Gaussian laser beam with waist w = 4 µ m is obliquelyincident with an angle θ on an inhomogeneous plasma slab. The density profile is invariant in directions x , y and3linearly increases in z -direction. Its profile is shown as a solid blue line. The laser pulse is p-polarized and itsamplitude is temporally described by sin (cid:0) π tT (cid:1) for t ≤ T and is zero for t > T . We set T = 40 fs in the amplitudeprofile which is equivalent to a Full Width at Half Maximum (FWHM) of 12.7 fs in the intensity profile. The beam isinvariant along y -direction. For reference, we display as black dashed lines the turning and critical surfaces at whichthe plasma density reaches respectively n cr cos θ and n cr , which is the critical density.We use the following initial densities and fields: • The initial plasma densities ρ i /m i and ρ e /m e have the following profile (in cm − ): z < µm . × ( z −
1) 1 µm ≤ z ≤ µm . × z > µm (40)This profile corresponds to a plasma density scale length of L = 3 . µ m, which is the length in z -directionat which the critical density is reached. A weak uniform background density of 10 cm − is added to avoiddivisions by zero ( e.g. when calculating velocity from momentum and density) and too strong discontinuity atthe ramp onset. • Initial electric and magnetic fields are uniformly zero. The electromagnetic pulse is injected at z = − µ m.The free-space wavelength is λ = 0 . µ m and the amplitude in free-space is E = 10 V/m.As we will further on investigate the conversion of energy between its different types in a finite volume, we definea volume V - shown as white dashed line in Fig. 5(a) - comprised between z = − µ m and z = 9 µ m and whichcontains all the numerical points in x - and y -directions. The local electromagnetic energy density is [46]: U em = 12 (cid:20) (cid:15) || E || + 1 µ || B || (cid:21) (41)The linear density of electromagnetic energy Σ U em is calculated by integrating the energy density U em over the volume V and dividing by the box thickness in y -direction.We also define the thermal energy density as (cid:80) k p k γ − and the kinetic energy density as (cid:80) k 12 ρ k || u k || . The corre-sponding linear density of thermal energy and kinetic energy are respectively written Σ U Thermal and Σ U Kinetic . Thelinear density of total energy is Σ (cid:15) = Σ U em + Σ U Thermal + Σ U Kinetic .In Fig. 5(b) we show the result of the simulation for an incidence angle θ = 15 ◦ . We plot the different types ofenergy defined above as a function of time. The linear density of energy of the input pulse in the x − z plane is givenby : E P ulse = E (cid:113) (cid:15) µ T w (cid:112) π = 10 mJ/m and is shown as a black dotted line.The electromagnetic energy increases as the pulse enters into the volume V between t = 0 and t <
45 fs. At t = 45 fs, the pulse is completely contained in the volume and has not yet interacted with the plasma ramp. Wesee that the pulse energy corresponds to the predicted value of 10 mJ/m. In the temporal window 45-130 fs, energyexchange with the plasma takes place (the kinetic energy of the plasma increases). Then, between 130 fs and 160 fs,the reflected pulse leaves the integration volume and the electromagnetic energy decreases. We see that the a fractionof energy remains in the volume. This residual energy corresponds to the one transferred from the electromagneticwave into the electron plasma wave. The conversion factor for this simulation is 46%. Since the collisions frequencyis taken to zero, the thermal energy linear density stays zero (dashed green line). We also see in the same graph thatthe total numerically integrated energy density, shown as solid purple line, actually corresponds to the incident pulseenergy when the pulse is fully inside the volume V . This shows that our algorithm preserves energy.We repeat the same numerical simulation with different angles of incidence. The conversion factor is shown in Fig. 6as a function of τ = (cid:0) πLλ (cid:1) / sin θ . The results of the other references are shown for comparison: the theoreticalasymptotes of Speziale (solid blue) [43] et al. , the analytical results from Hinkel-Lipsker et al. (green dotted) [45] andthe numerical results of Forslund et al. (orange dashed) [44] (the latter were obtained for a relatively cold plasmawith k B T e / ( m e c ) = 0 . et al. [44]. Our numerical results, shown as blue crosses, are in excellent agreement with all these previous results.4 C on v e r s i on f a c t o r Numerical pointsSpeziale asymptotesForslund (numerical)Hinkel-Lipsker (analytical)
FIG. 6. Mode conversion factor as a function of τ = (cid:0) πLλ (cid:1) / sin θ . Blue crosses are our numerical results. Results fromother references are shown for comparison: Forslund PIC simulations (dashed orange line) [44], Hinkel-Lipsker analytical model(dotted green line) [45] and Speziale asymptotes (solid blue lines) [43]. -5 0 5 10 15 z (µm) -1-0.500.51 E x ( V / m ) Numerical fieldNumerical envelope ei (fs -1 ) z / e ( µ m ) (a) (b) Theory : = 0.3 10 cm -3 Numerical : = 0.3 10 cm -3 Theory : =0.9 10 cm -3 Numerical : = 0.9 10 cm -3 FIG. 7. (a) Spatial profile of electromagnetic wave amplitude propagating into a homogeneous collisional plasma of density0 . × cm − and a collision time of 2 fs − . The vacuum/plasma interface is at z = 0 µ m. (b) Numerical and theoreticalcharacteristic decay length z / e as function of collision frequency ν ei . Results are plotted for two different densities: 0 . × cm − (blue line, orange crosses) and 0 . × cm − (yellow line, purple crosses). Test 4: Collisional plasmas
Here we test the integration of the BGK collisional model. The Drude model for a cold uniform plasma of heavyions provides the complex permittivity [48]: (cid:15) = 1 − ω τ ω τ + i ω τ c ω (1 + ω τ ) (42)We link the collision time τ c of Drude model to the electron-ion collision frequency of BGK model via τ c = ν ei .We check our model through the collisional damping of a monochromatic plane wave impinging on a plasma slab ofuniform density. The distance at which the wave amplitude is reduced by a factor of 1/e is given by: z /e = c (cid:61) ( √ (cid:15) ) ω (43)5Fig. 7 shows the results obtained for the following initial conditions: • Initial plasma density is uniform and set to 0 . × cm − and to 0 . × cm − in a second series ofsimulations. The plasma is only in region z ≥ µ m. • The collision frequency ν ei is varied in the range 0.25-20 fs − . • The electric and magnetic fields are initially set to zero. A monochromatic electromagnetic plane wave propa-gating along the z -axis is progressively branched with normal incidence on the plasma. It is polarized along the x -axis with λ = 0 . µ m and amplitude of E = 10 V/m.In Fig. 7(a), we plot the amplitude of the electric field as a function of propagation distance z . The exponentialdamping is observed as expected. In Fig. 7(b), we report the values of the decay length z / e which we compare to theanalytical results from Eq. (43). The error bars correspond to the graphical measurement. An excellent agreementis found. In parallel, we measured the evolution of wavelength in same conditions, which we compared with thetheoretical value 2 πc/ [ ω (cid:60) ( √ (cid:15) )] (not shown). The discrepancy between numerical and theoretical results were lessthan 1%.Separately, we also test the conservation of energy within the collisional plasma. The same simulation as forFig. 5(b) is performed in the collisional plasma with a collision frequency ν ei = 0 .
05 fs − (see Fig. 8). We observe asimilar evolution as in section , except that now, due to the collisions, the energy of waves is converted to thermalenergy with time. The total density of energy is preserved as one can see from the purple line in the temporal window[45-130] fs, when the pulse is entirely in the integration volume V . It matches the initial pulse energy density. t (fs)0246810 E ne r g y den s i t y ( J / m )
500 100 150 200 -3 U em U Kinetic U Thermal E Pulse
FIG. 8. Temporal evolution of the linear densities of: electromagnetic energy (dashed blue line), kinetic energy (dashed-dottedred line), thermal energy (dotted green line), and total energy (solid purple line) integrated over the volume V for the laserpulse described in section and a collision frequency ν ei = 0 .
05 fs − . The initial pulse linear density of energy E Pulse is plottedas a dashed black line.
Overall, these results validate the integration of the collisional model.
Test 5: Wakefield generation
Here, we test the algorithm in a nonlinear regime, where transport effects are non-negligible. For this, we considerthe interference between two contrapropagative laser pulses in a sub-critical plasma. The interference generates aperiodic pattern of low and high intensities along the beams axis. The ponderomotive force swiftly pushes the electronsaway from high-intensity regions [49], which creates a net space charge and consequently a wakefield. The use of astanding wave therefore allows the observation of strong ponderomotive effects without using relativistic laser pulses[50].We reproduce this physical situation in this test. We consider the following initial conditions:6 -1 -0.5 0 0.5 1 z (µm) E l e c t r on den s i t y ( c m - ) -202 E ( V / m ) z FIG. 9. Electron density profile (dashed blue line) and wakefield E z (red line) at the time where the two interfering pulsescross at z = 0. • The sub-critical plasma has initial uniform densities ρ i /m i = ρ e /m e = 5 × cm − and the collision frequencyis ν ei = 0 . − . • The electric and magnetic fields are initially set to zero in the plasma. Two incident electromagnetic pulses areinjected to generate a standing plane wave. The first one is injected in z = − µ m and is propagating along z -direction. The second one is injected in z = 20 µ m and is propagating in the opposite direction. Both of themare polarized along x axis with λ = 0 . µ m and amplitude of E = 1 × V/m. The intensity profile of bothpulses is a (sin ) with FWHM 47.6 fs.The simulation is run until the pulses peaks cross at z = 0. Fig. 9 shows, for that time, the electron density profile(dashed blue line) and the longitudinal component of the electric field E z as a function of the propagation distance z .Our results fit qualitatively very well with reference [4]. On a quantitative point of view, reference P.-W Smorenburg et al. [50] provide an estimate of wakefields of ∼ V/m for standing wave field strength of the order of ∼ W/cm .Thus, we performed the same simulation as for Fig. 9 but in a collionless plasma and with a field which correspondsto a standing wave with an intensity of 10 W/cm . We obtained a wakefield E z = 5 × V/m. Our results aretherefore in good quantitative agreement with Ref. [50]. In addition, we note that very strong modifications of thedensity profile (up to 50%) are reached without compromising the stability of the simulation.
CONCLUSIONS
In conclusion, we have developed a simple two-fluid plasma model solver. Our algorithm is composed of a PSTDsolver for Maxwell’s curl equations and of a combination of LWLFn and RK4 algorithms via a Strang splitting forthe fluid equations. We have also integrated a simple collisional model. We have validated the solver against severalwell-known problems of laser-plasma interaction. In all cases, very good agreement has been achieved.This solver is simple and robust. Its implementation is relatively straightforward since, in contrast with complexitiescontained in most of the two-fluid plasma solvers, it relies only on finite difference schemes and FFT. This type ofsolver is a good compromise for physical problems that do not involve strong discontinuities. It can represent a simplefirst step before investing time and efforts into more elaborate solvers. We anticipate that this approach will opennew perspectives in the fields of nonlinear plasmonics and laser-plasma interaction in solid.
ACKNOWLEDGMENTS
Numerical simulations have been performed using the M´esocentre de Calcul de Franche-Comt´e. The research leadingto these results has received funding from the European Research Council (ERC) under the European Union’s Horizon72020 research and innovation program (grant agreement No 682032-PULSAR), R´egion Bourgogne Franche-Comt´e andthe EIPHI Graduate School (ANR-17-EURE-0002). ∗ [email protected][1] S. Eliezer and K. Mima, Applications of Laser-Plasma Interactions (CRC Press, 2008).[2] W. Kruer,
The Physics Of Laser Plasma Interactions (Westview Press, 2003).[3] P. McKenna, D. Neely, R. Bingham, and D. Jaroszynski,
Laser-Plasma Interactions and Applications (Springer, 2013).[4] P. Gibbon,
Short Pulse Laser Interactions with Matter: An Introduction (Imperial College Press, 2005).[5] R. Abgrall and H. Kumar, J. Sci. Comput. , 584 (2014).[6] R. J. Mason, J. Comput. Phys. , 429 (1987).[7] R. J. Mason and C. W. Cranfill, IEEE T. Plasma Sci. , 45 (1986).[8] K. Yee, IEEE Trans. Antennas Propag. , 302 (1966).[9] U. Shumlak and J. Loverich, J. Comput. Phys. , 620 (2003).[10] P. L. Roe, J. Comput. Phys. , 357 (1981).[11] J. Loverich and U. Shumlak, Comput. Phys. Commun. , 251 (2005).[12] J. Loverich, A. Hakim, and U. Shumlak, Commun. Comput. Phys. , 240 (2011).[13] B. Srinivasan and U. Shumlak, Phys. Plasmas , 092 (2011).[14] E. M. Sousa and U. Shumlak, J. Comput. Phys. , 56 (2016).[15] A. Alvarez Laguna, N. Ozak, A. Lani, H. Deconinck, and S. Poedts, Comput. Phys. Commun. , 31 (2018).[16] V. V. Rusanov, J. Comput. Math. Phys. , 267 (1961).[17] S. Baboolal and R. Bharuthram, Math. Comput. Simul. , 3 (2007).[18] H. Kumar and S. Mishra, J. Sci. Comput. , 401 (2012).[19] M. Brio and C. C. Wu, J. Comput. Phys. , 400 (1988).[20] A. Hakim, J. Loverich, and U. Shumlak, J. Comput. Phys. , 418 (2006).[21] P. L. Bhatnagar, E. P. Gross, and M. Krook, Phys. Rev. , 511 (1954).[22] V. E. Golant, A. P. Zhilinsky, and I. E. Sakharov, Fundamentals of Plasma Physics (John Wiley & Sons, 1980).[23] Q. H. Liu, Microw. Opt. Technol. Lett. , 158 (1997).[24] H. Vincenti and J.-L. Vay, Comput. Phys. Commun. , 22 (2018).[25] J. Shang and R. M. Fithen, J. Comput. Phys. , 378 (1996).[26] C. D. Munz, P. Ommes, and R. Schneider, Comput. Phys. Commun. , 83 (2000).[27] A. Monorchio, A. R. Bretones, R. Mittra, G. Manara, and R. G. Martin, IEEE Trans. Antennas Propag. , 2666 (2004).[28] J.-P. Berenger, J. Comput. Phys. , 185 (1994).[29] T.-W. Lee and S. C. Hagness, J. Opt. Soc. Am. B , 330 (2004).[30] J. A. Bittencourt, Fundamentals of Plasma Physics (Springer, 2013).[31] E. F. Toro,
Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction (Springer, 2013).[32] R. Liska and B. Wendroff, SIAM J. Numer. Anal. , 2250 (1998).[33] T. R. Hagen, M. O. Henriksen, J. M. Hjelmervik, and K.-A. Lie, How to Solve Systems of Conservation Laws NumericallyUsing the Graphics Processor as a High-Performance Computational Engine (Springer, 2007).[34] R. Liska and B. Wendroff,
Composite Centered Schemes for Multidimensional Conservation Laws (Birkhauser, 1999).[35] M. Kucharik, R. Liska, S. Steinberg, and B. Wendroff, Appl. Numer. Math. , 589 (2006).[36] R. J. Leveque, Finite Volume Methods for Hyperbolic Problems (Cambridge University Press, 2002).[37] L. F. Shampine, Appl. Math. Lett. , 1134 (2005).[38] P. Lax and B. Wendroff, Comm. Pure Appl. Math. , 217 (1960).[39] D. Richtmyer, A Survey of Difference Methods for Non-Steady Fluid Dynamics (Technical Note NCAR/TN-63-2+STR,1962).[40] G. Strang, SIAM J. Numer. Anal. , 506 (1968).[41] J. C. Butcher, Numerical Methods for Ordinary Differential Equations (John Wiley & Sons, 2004).[42] R. Liska and B. Wendroff,
Comparison of Several Difference Schemes for the Euler Equations in 1D and 2D (Springer,2003) pp. 831–840.[43] T. Speziale and P. J. Catto, Phys. Fluids , 990 (1977).[44] D. W. Forslund, J. M. Kindel, K. Lee, E. L. Lindman, and R. L. Morse, Phys. Rev. A , 679 (1975).[45] D. E. Hinkel-Lipsker, B. D. Fried, and G. J. Morales, Phys. Rev. Lett. , 2680 (1989).[46] D. J. Griffiths, Introduction to Electrodynamics (Prentice Hall, 1999).[47] G. J. Pert, Phys. Plasmas , 175 (1978).[48] N. W. Ashcroft and N. D. Mermin, Solid State Physics (Cengage Learning, 2011).[49] P. Mulser and D. Bauer,
High Power Laser-Matter Interaction (Springer, 2010).[50] P. W. Smorenburg, J. H. M. Kanters, A. Lassise, L. P. J. Brussaard, and O. J. Luiten, Phys. Rev. A83