A solver based on pseudo-spectral analytical time-domain method for the two-fluid plasma model
AA solver based on pseudo-spectral analytical time-domain method for the two-fluid plasma model
B. Morel ∗ , R. Giust, K. Ardaneh and F Courvoisier ∗ FEMTO-ST institute, Univ. Bourgogne Franche-Comt´e, CNRS,15B avenue des Montboucons, 25030, Besanc¸on Cedex, France ∗ Corresponding authors: [email protected] and [email protected] is a pre-peer-review version of an article published in
Scientific Reports . The final version is available online at: https://doi.org/10.1038/s41598-021-82173-9
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 particularlysuited 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 PseudoSpectral AnalyticalTime-Domain (PSATD) method to solve Maxwell’s curl equations. Overall, this method is only based on finitedifference schemes and fast Fourier transforms and does not require any grid staggering. The PseudoSpectralAnalytical Time-Domain method removes the numerical dispersion for transverse electromagnetic wave propa-gation in the absence of current that is conventionally observed for other Maxwell solvers. The full algorithm isvalidated by conservation of energy and momentum when an electromagnetic pulse is launched onto a plasmaramp and by quantitative agreement with wave conversion of p-polarized electromagnetic wave onto a plasmaramp.
INTRODUCTION
Chirped pulse amplification technology in 1985 [1] hasmade possible the generation of extremely powerful laserpulses [2]. When a solid, liquid or a gas is irradiated by sucha powerful pulse, the ionization phenomena swiftly create aplasma at the surface of the material or within the gas. Thedevelopment of applications such as inertial fusion [3], laser-plasma accelerators [4], laser materials processing [2], X-Raylasers [5] or nonlinear plasmonics at lower intensities[6] re-quires laser plasma interactions modeling.Hydrodynamic models are particularly useful to describeshort pulse interaction with plasmas when each of the speciescan be assumed in local thermodynamic equilibrium [7]. Thetwo-fluid plasma equations is the starting point of the hy-drodynamic models [8]. This model describes the spatio-temporal evolution of the density, mean velocity and pres-sure of electrons and ions fluids. The two-fluid plasma equa-tions therefore consist of two sets of Euler equations withsource term, as well as Maxwell’s equations. The fluid de-scription involves the assumptions of local thermodynamicequilibrium for each species (electrons, ions). The conven-tional hydrodynamic models, e.g. , two-temperature plasmaequations, single-fluid equations and MagnetoHydroDynamic(MHD) can be derived from the two-fluid plasma model bymeans of additional assumptions. At present, solving the complete two-fluid plasma equa-tions is a difficult challenge[9]. Their implementation is oftencomplex for non-specialist groups since most of these codesare developed to be particularly robust for shock’s problems.A good example is given in Shumlak et al. [10] which presentan algorithm based on Roe-type Riemann solver [11] for thetwo-fluid plasma model. The same group added the high-order discontinuous Galerkin method to improve the result’saccuracy [12–15]. References [16–20] describe numericalmethods well adapted for shock’s problems. In contrast, forproblems without strong shocks, our group has recently pro-posed a relatively simple approach [21], based on finite differ-ence schemes and Fast Fourier Transform (FFT). This solveris based on the PseudoSpectral Time-Domain method (PSTD)[22] to solve Maxwell equations and a composite scheme [23]to solve the fluid equations. However, the PSTD is based ona temporally staggered grid, which requires temporal inter-polations for the coupling with the fluid solver. In addition,the PSTD algorithm is numerically dispersive. It emits wavecomponents that are faster than light.J.-L. Vay et al. [24] proposed the PseudoSpectral Analyt-ical Time-Domain (PSATD), and its application to pseudo-spectral Particle-In-Cell (PIC) simulations. In PSTD, the tem-poral integration is performed via finite differences, while inPSATD, the integration is analytical (except for the integra-tion of current). Thus, unlike PSTD, the PSATD requires no a r X i v : . [ phy s i c s . p l a s m - ph ] F e b temporal staggered grid and is free of numerical dispersionfor transverse electromagnetic propagation in the absence ofcurrent (see Fig. 1 of J.L. Vay et al. [24]). This method is alsoparticularly well adapted for laser pulse propagation. The al-gorithm is tested with laser/plasma interaction problems withintensities around to W/cm , since it is intended forthe study of electron/hole plasma dynamics in solids.Here, we build a two-fluid plasma solver based onPseudoSpectral Analytical
Time-Domain PS A TD for solv-ing Maxwell’s equations. A schematic representation of oursolver is given in Fig. 1. The integration of Maxwell’sequations is performed by using the PSATD method. Theelectromagnetic fields are transmitted to the fluid equationsas a Lorentz force source term. The fluid equations are in-tegrated 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 vari-ables are used to calculate the current density, which is in-jected in Maxwell’s equations. LW Fluid equations
PSATD RK4
Maxwell equations
Strang splitting
FIG. 1. Schematic representation of the solver’s structure.
This paper is divided in four main parts. We will first recallthe two-fluid plasma model equations before summarizing thenumerical integration. We will then validate the solver anddemonstrate its benefits in terms of numerical dispersion andin the reduction of the constraint imposed on time-steps withthe solver of reference [21].
RESULTSTwo-fluid plasma equations
The two-fluid plasma model equations consists of Eu-ler equations with source term for each fluid, as well asMaxwell’s equations. This system of equations correspondsto continuity equations, motion equations and energy trans-port equations for electron and ion fluids. Following reference[26], the fluid equations can be presented under the following form: ∂∂t ρ e ρ e u e (cid:15) e ρ i ρ i u i (cid:15) i (cid:124) (cid:123)(cid:122) (cid:125) ≡ 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 (cid:124) (cid:123)(cid:122) (cid:125) ≡ F ( U ) = ρ 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 (cid:124) (cid:123)(cid:122) (cid:125) ≡ S ( U , E , B ) (1)where U is the fluid variables vector, F ( U ) is the flux ten-sor and S ( U , E , B ) is the Lorentz force source term. In thispaper, i and e are indexes related respectively to the ion fluidand to the electron fluid. q is the charge, m the mass, ρ themass density, u the mean velocity, p the pressure, (cid:15) the fluidenergy 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 equa-tion, an ideal gas closure for each fluid k is used [26]: (cid:15) k ≡ p k γ − ρ k u k2 (2)where γ is the adiabatic index.Electric and magnetic fields in the source term S ( U , E , B ) are determined by the Maxwell equations. Since solving fluidequations and Maxwell’s curl equations enforces the conser-vation of divergence properties of the fields [27], it is notnecessary to solve Maxwell’s divergence equations. Howeverthey must be satisfied at an initial time. Furthermore, Maxwellcurl’s equations can be written by expressing the current den-sity J as function of fluid variables: ∇× E = − ∂ B ∂t (3) ∇× 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) ≡ J + ε r c ∂ E ∂t (4)Here, ε and µ are respectively the vacuum permittivity andpermeability. c = ( ε µ ) − / is the speed of light. ε r is therelative permittivity of the background medium: we do theassumption that this quantity is time and space independent.In this model, the plasma is also contained inside a mediumof 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 methodis simple to implement and does not need the staggering ofspatial and temporal grids. This is in contrast with the FiniteDifference Time Domain (FDTD) method [28] which requiresspatially and temporal staggered grids or in contrast with thePseudoSpectral Time-Domain (PSTD) [22] which requires atemporally staggered grid. The PSATD is therefore more flex-ible to be coupled with another algorithm without interpola-tions. Moreover, in absence of current, PSATD induces zeronumerical dispersion in contrast with FDTD or PSTD. An ad-ditional strong benefit is that the PSATD is not subject to aCourant condition for transverse electromagnetic field propa-gation in the absence of current. The PSATD algorithm is in-herently periodic because it is based on FFT, but open systemscan be modeled by using Perfectly Matched Layers (PML) asin O. Shapoval et al. [29].The PSATD algorithm provides the fields in the Fourierspace [24]: ˜E n+1 = C ˜E n + ivS κ × ˜B n − ε ε r S kv ˜J n+1 / +(1 − C ) κ · ( κ · ˜E n ) + 1 ε ε r (cid:18) S kv − ∆ t (cid:19) κ · ( κ · ˜J n+1 / ) (5) ˜B n+1 = C ˜B n − i S v κ × ˜E n + iµ − C k κ × ˜J n+1 / (6)where ˜a is the Fourier transform of the quantity a . Here C = cos ( kv ∆ t ) , S = sin ( kv ∆ t ) , κ = k /k and v = cε / r .The two main assumptions made in the PSATD method are:1) the time-step ∆ t is enough small to assume that the currentdensity is constant over a time-step 2) the background permit-tivity ε r is uniform. The fluid solver
For the fluid equations solver, we consider a similar solveras in reference [21]. Here, we simplified the solver by re-stricting ourselves to problems without discontinuities suchthat it becomes unnecessary to introduce numerical dissipa-tion to make gradient smoother. Instead of using a compos-ite scheme LWLFn as in reference [21], we will use a simpletwo-step Lax-Wendroff (LW) scheme [30] which is second or-der accurate and introduces less numerical dissipation than thetwo-step Lax-Friedrichs (LF) scheme [31]. The LW schemesolves 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: L x ( U nj , l , m ) = U nj , l , m − ∆ t ∆ x × (cid:104) F x (cid:16) U n+1 / / , l , m (cid:17) − F x (cid:16) U n+1 / − / , l , m (cid:17)(cid:105) (7)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 (cid:0) U nj+1 , l , m (cid:1) − F x (cid:0) U nj , l , m (cid:1)(cid:3) (8) U n+1 / − / , l , m = 12 (cid:2) U nj , l , m + U nj − , l , m (cid:3) − ∆ t x (cid:2) F x (cid:0) U nj , l , m (cid:1) − F x (cid:0) U nj − , l , m (cid:1)(cid:3) (9)where j , l and m are respectively indexes for x , y and z di-rections. Similar operations are done in y and z directionsas reference [21] to obtain L y ( U nj , l , m ) and L z ( U nj , l , m ) . A ba-sic spatially dimensionally-split scheme is used to obtain thevalue U n+1j , l , m from U nj , l , m [32]: U n+1j , l , m = L x L y L z ( U nj , l , m ) (10)For the numerical integration of the source term S ( U , E , B ) of Eq. (1), we use the Strang splitting techniquepresented by G. Strang[25]. The Strang splitting allows anestimation of current density J n+1 / at a half time step ofPSATD. The concept of Strang splitting is shown on the steps1, 2 and 4 in Fig. 2. We first integrate the source term with anRK4 scheme over ∆ t/ , then the homogeneous system is in-tegrated over ∆ t with an LW scheme, and finally source termis again integrated with an RK4 over time step ∆ t/ . Full two-fluid plasma solver algorithm
The full algorithm for the two-fluid plasma model is de-scribed in Fig. 2 and can be decomposed in 4 main steps:1. Integration of the source term with an RK4 scheme overa temporal step ∆ t/ by using E n , B n and U n to obtainthe intermediate value of fluid variables U ∗ .2. Integration of the homogeneous system with an LWscheme over a temporal step ∆ t using fluid variablesvector U ∗ to obtain a new intermediate value U ∗∗ .3. Computation of the current density J n+1 / with densi-ties and velocities from U ∗∗ . Then, carry out a PSATDstep with J n+1 / to calculate E n+1 and B n+1 .4. Integration of the source term with an RK4 algorithmover a temporal step ∆ t/ 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 nodesof a grid, it also avoids temporal interpolation of the magneticfield that was necessary in reference [21].For the PSATD algorithm alone without currents, the sam-pling is in principle only limited by Nyquist theorem. How-ever, in order to derive Eqs. (5) and (6), we make the assump-tion that the current is constant over the temporal step ∆ t .Therefore, the temporal step is chosen small enough to makethis assumption valid. The spatial step ∆ x is simply chosento resolve the both plasma and electromagnetic waves. t n U n E n B n E n U * E n B n U n+1 E n+1 B n+1 U ** E n B n U ** E n+1 B n+1 t n+1/2 t n+1 t RK4 Δ t/2 RK4 Δ t/2LW Δ t U * E n B n U ** B n U ** E n+1 B n+1 PSATD Δ t Source termHomogeneous system Maxwellequations Initial state Intermediate state Final state
FIG. 2. Schematic representation of the algorithm. The source term is integrated with an RK4 scheme while the homogeneous system isintegrated with an LW scheme. Maxwell’s equations are solved with the PSATD method. The algorithm requires four steps shown in red toadvance fluid variables and fields from a time step n to n + 1 . Validation of the numerical solver
S-polarized electromagnetic wave over a plasma ramp
In this first test, we check the conservation of momentumand energy during reflection of a s-polarized electromagneticwave over a plasma ramp. The numerical setup is shown inFig. 3. A laser pulse is propagating toward an overcriticalplasma ramp with an angle of incidence θ = 15 ◦ . The ini-tial plasma density profile is invariant in y and z directions,and the following initial density profiles, for electron and ionfluids, are used in x direction: n e = n i = for x > − µm . × (5 − x ) for − µm ≤ x ≤ − µm . × for x < − µm (11)where n e ≡ ρ e /m e and n i ≡ ρ i /m i are given in cm − .The length in x direction at which the critical density n c =1 . × cm − is reached is L = 3 . µ m. We add a weakuniform background density of cm − to avoid divisionsby zero in the algorithm and too strong discontinuity at theramp onset. In this test, the uniform background is vacuum: ε r = 1 . For the plasma, we take m e = 9 . × − kg, m i = 1837 m e and γ = 5 / . The initial mean velocities andpressure are zero.The laser pulse is a spatially Gaussian beam with a waist w = 4 µ m and is described temporally by a single period ofa sin function (period T = 40 fs). The free-space wave-length is λ = 0 . µ m and the amplitude in free-space is E = 4 . × V/m. We choose this electric field amplitudeto demonstrate the possibility of working with high field am-plitudes with this algorithm. Note that the beam is invariantalong 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) bound-ary conditions in x and y directions for PSATD (resp. fluidalgorithm) are implemented. For fluids and fields, periodicboundary condition are used in z direction. The spatial stepis ∆ x = ∆ y = ∆ z = 60 nm and the temporal step is ∆ = 96 as.In Fig. 4(a), we plot the different momenta in x directionas function of the simulation time. These momenta are nor-malized to the absolute value of the x momentum P of theincident pulse. To measure the normalization factor P , weperformed beforehand the simulation without the plasma, andwe measure the x momentum of the pulse defined by the x component of the Eq. (12) integrated over the simulation win-dow.The density of electromagnetic momentum is defined by[33]: P em = ε E × B (12)The dashed red curve of Fig. 4(a) corresponds to the nor-malized electromagnetic x momentum density integrated in the simulation window. The dashed dotted blue curve corre-sponds to the normalized fluids x momentum integrated in thesimulation window. The momentum of fluids is defined by: P f = ρ e u e + ρ i u i (13)The black line of Fig. 4(a) is the sum of the electromagneticand fluids x momentum.We observe three main sequences in Fig. 4(a):• : Since the laser pulse goes from the right to the left,the electromagnetic momentum along x (red dashedcurve) decreases as the pulse enters into the simulationwindow (between t = 0 and t < fs). At t = 50 fs,the pulse is completely contained in the simulation win-dow and has not yet interacted with the plasma ramp.We see that the electromagnetic x momentum corre-sponds to the incident pulse x momentum − P .• : In the temporal window 70-130 fs, momentum ex-change with the plasma takes place: the fluids momen-tum decreases until − P , whereas the electromagnetic x momentum increases until reaches + P . This is thesignature that the laser pulse transfers twice its initialmomentum to the plasma during its reflection, as canbe expected.• : Between fs and fs, the reflected pulse leavesthe simulation window, thus the electromagnetic mo-mentum goes back to zero.We see that the momentum is preserved over the temporalwindow over which the pulse is fully enclosed within the sim-ulation window. The error on the conservation of the totalmomentum is slightly less than 1%. It is reasonable in viewof the chosen spatial and temporal steps. The numerical al-gorithm also preserves the conservation of momentum with agood accuracy.In Fig. 4(b), we plot the linear density of energy as functionof simulation time. The dashed red curve correspond to theelectromagnetic energy density [33] U em = 12 (cid:20) ε E + 1 µ B (cid:21) (14)that we have integrated over the x − y plane. The dasheddotted blue (resp. dashed green) corresponds to the electronfluid (resp. ion fluid) energy density given by Eq. (2) and thenintegrated over the x − y plane. The total energy plotted inblack line is defined as the sum of electromagnetic, electronfluid and ion fluid energies. The linear density of energy of theinput pulse in the x − y plane can be calculated analyticallyand is given by: E Laser = E (cid:113) ε µ T w (cid:112) π = 0 . J/m.This analytical linear density of energy is shown as a blackdotted line in Fig. 4(b).We observe the three main sequences in Fig. 4(b):• : The electromagnetic energy increases as the pulseenters into the numerical window between t = 0 and k E B θ L a s e r p u l s e Electron/ion Plasma P M L / O p e n b o un d a r y c o n d i t i o n PML/Open boundary conditionPML/Open boundary condition P M L / O p e n b o un d a r y c o n d i t i o n Density x Criticaldensity n c xyz P u l s e i n j e c t o r n c cos² θ FIG. 3. Numerical setup: A laser pulse in oblique incidence is injecting toward an overcritical plasma ramp. The laser pulse is reflected at theturning surface such as density n = n c cos θ . t (fs) a) b) t (fs) -2-1.5-1-0.500.51 P x / P E n e r g y d e n s i t y ( J / m ) FluidsElectromagneticTotal ElectromagneticElectron fl uidIon fl uidTotalAnalytical Pulse fully in the window Pulse fully in the window
FIG. 4. (a) Conservation of momentum and (b) conservation of energy during s-polarized pulse reflection over the plasma ramp. In the centralwhite area, the pulse is fully in the numerical window and conservation of momentum and energy are preserved. t < fs. At t = 50 fs, the pulse is completelycontained in the numerical window and has not yet in-teracted with the plasma ramp. The electromagneticenergy corresponds to the predicted analytical value E Laser .• : In the temporal window 70-130 fs, energy exchangewith the plasma takes place (the electron energy in-creases).• : Between fs and fs, the reflected pulse leavesthe integration volume and the electromagnetic energydecreases. No electromagnetic energy remains in thesimulation window. This is expected for s-polarizedwave.We remark the conservation of the energy when the pulse isfully in the simulation window. The error on the conservationthe total energy is around 0.1%. The numerical algorithm alsopreserves the energy conservation with a good accuracy.Figs. 4(a) and 4(b) demonstrated conservation of momen-tum and energy during s-polarized reflection over a plasmaramp. Wave conversion on plasma density ramp
In this second test, we consider the same numerical setupas shown in Fig. 3, but we inject a p-polarized laser pulse.The energy of the system is plotted as a function of timein Fig. 5(a). In the central white area, the error on the con-servation the total energy is around 0.1%. Furthermore, weobserve for the sequence n ◦ i.e. conversion of an electromagneticwave into a plasma wave which occur only for p-polarization[34].The conversion factor depends in particular on the plasmadensity gradient and the angle of incidence. Obtaining analyt-ical solutions to this difficult problem usually requires a num-ber of approximations. We performed a series of simulationswith different angles of incidence, and we plot (blue circles)in Fig. 5(b) the factor of the energy conversion as function ofthe quantity τ = (cid:0) πLλ (cid:1) / sin θ .We compare conversion factors obtained with thePSATD/Hydrodynamic code (this work) and results of the lit-erature. The PSATD/Hydro conversion factor curve is quan-titatively superimposable to the one obtained the PSTD/hydrosimulation (blue crosses) obtained in reference [21]. Our re-sults are also in agreement with the analytical results of D.E.Hinkel-Lipsker et al. [35], with T. Speziale et al. who de-scribed the asymptotic behaviors [34] and also with those ofD.W. Forslund et al. who have used Particle-In-Cell (PIC)simulations [36]. The fact that we injected a short pulse (poly-chromatic) gaussian beam instead a monochromatic plane wave can explain the tenuous differences. In addition, the an-alytical results of references T. Speziale et al. [34] and D.E.Hinkel-Lipsker et al. [35], have carried out assumptions thatare not exactly fulfilled in our numerical test. This can explainthe minor discrepancies observed. But overall, the results weobtained with the present solver are in good agreement withthe state of the art. DISCUSSION
In this section, we compare the benefits and drawbacks ofPSATD/Hydro solver compared to PSTD/Hydro solver. Wenumerically simulate a single cycle pulse plane wave in nor-mal incidence onto a plasma ramp. The laser wavelength andplasma parameters are identical to the ones of Fig. 3. Thepulse amplitude is E = 4 . × V/m. The computation isperformed in 3D, with the same numerical sampling param-eters as in Fig. 3. We use the periodic boundary conditionsin y and z directions. For the PSTD/Hydro solver, the time-step is fixed to ∆ t = 50 as since we are constraint by CourantFriedrichs Lewy (CFL) conditions [22]. In contrast, the time-step for the PSATD/Hydro solver is set to ∆ t = 200 as, as it isonly constrained by the sampling of the laser and the plasmawave frequencies.In Fig. 6(a), we show a snapshot during propagation ofthe laser pulse with the different solvers. The snapshot istaken when the pulse has propagated through vacuum andjust reaches the onset of the plasma ramp. We see that thePSATD/Hydro solver result (solid blue line) is precisely su-perimposed on the analytical solution in black dashed line. Incontrast, the PSTD/Hydro solver (red dashed-dotted line) ex-hibits distortion of the laser pulse. Indeed, pre-pulses are gen-erated by numerical dispersion of the PSTD algorithm. Theamplitude of the last artifact pre-pulse (located at x ≈ − µ m)reaches around 15% of the amplitude of the main peaks.In Fig. 6(b), we plot the velocity component v z of the elec-tron fluid at the same time of the snapshot of Fig. 6(a). We ob-serve that the laser pulse has not yet interacted with the plasmain the PSATD/Hydro simulation (blue line). However, in thePSTD/Hydro simulation (red dashed-dotted line), the artifactpre-pulses already interact with the plasma and accelerate theelectrons to velocities around m/s. This effect is obvi-ously undesirable, particularly in the case of the simulation offew cycle laser pulses with plasmas [37].We also obtained better results in the PSATD/Hydro simu-lation whereas the time-step ∆ t was 4 times greater than thosein PSTD/Hydro simulation.The PSATD/Hydro solver is well suited to pulse propaga-tion. Specifically, the fact that PSATD is not constraint bythe CFL condition, releases the strong numerical link betweenspatial and temporal sampling. The computational gain istherefore particularly significant in the case where high spa-tial resolution is required together with less demanding tem-poral resolution. We finish this section by reminding that thePSATD method requires that the background medium permit- t (fs) E n e r g y d e n s i t y ( J / m ) Pulse fully in the window a) b)
ElectromagneticElectron fl uidIon fl uidTotalAnalytical [36][35] [34][21]) FIG. 5. (a) Evolution of electromagnetic energy, fluid energy, and total energy as function of time for p-polarized incident pulse. In the centralwhite area, the pulse is fully enclosed within the numerical window and the conservation of the total energy is preserved. (b) Mode conversionfactor as a function of τ = (cid:0) πLλ (cid:1) / sin θ . Green circles are our numerical results with PSATD/Hydrodynamic code. Results from othersreferences are shown for comparison: Hinkel-Lipsker analytical model (dotted black line) [35], Speziale asymptotes (solid blue lines) [34],Forslund PIC simulations[36] (dashed orange line) and PSTD/Hydrodynamic (Blue crosses) [21]. tivity 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 ben-efits of fundamentally having no numerical dispersion forthe propagation of electromagnetic waves in absence of cur-rent. The algorithm relies on the PseudoSpectral Analyti-cal Time-Domain (PSATD) technique which is a powerfulmethod for propagating laser pulses, and on a combination oftwo-step Lax-Wendroff (LW) and fourth-order Runge-Kutta(RK4) for the fluid equations. We have demonstrated that thePSATD/Hydro solver preserves momentum and energy dur-ing a test with s-polarized laser pulse incident over a plasmaramp. The tests of wave conversion on plasma ramps havedemonstrated an excellent quantitative agreement with numer-ical and analytical results of the state of the art. We haveshown that PSATD/Hydro solver has two main advantagescompared to the PSTD/Hydro solver: the pre-pulses gener-ated by numerical dispersion are removed and the time-step isnot constraint by CFL conditions. For simulations which re-quire low temporal resolution and high spatial resolution, thegain in terms of computational resources with PSATD/Hydrosolver can be really significant. The PSATD/Hydro solver isa computationally inexpensive but powerful tool for the studyof laser-plasma interaction. METHODS
Simulations were performed with Nvidia Tesla K40 GPUcard. This card has 12 GB memory size, 2880 CUDA coresand 745 MHz processor core clock.
ACKNOWLEDGEMENTS
Numerical simulations have been performed using theM´esocentre de Calcul de Franche-Comt´e. The research lead-ing to these results has received funding from the EuropeanResearch Council (ERC) under the European Union’s Horizon2020 research and innovation program (grant agreement No682032-PULSAR), R´egion Bourgogne Franche-Comt´e, theEIPHI Graduate School (ANR-17-EURE-0002) and I-SITEBFC project (ANR-15-IDEX-0003).
AUTHOR CONTRIBUTIONS STATEMENT
B.M. conceived the algorithm structure, implemented thefluid solver and performed the simulations, R.G. implementedthe PSATD solver, K.A. and F.C. analysed the results. B.M.and F.C. prepared the manuscript which was reviewed by allauthors. [1] Strickland, D. & Mourou, G. Compression of ampli-fied chirped optical pulses.
Opt. Commun. , 447–449 (1985). URL .[2] Eliezer, S. & Mima, K. Applications of Laser-Plasma Interac-tions (CRC Press, 2008).[3] Basov, N. & Krokhin, O. Condition for heating up of a plasmaby the radiation from an optical generator.
J. Exp. Theor. Phys , 123–125 (1964).[4] Esarey, E., Sprangle, P., Krall, J. & Ting, A. Overview ofplasma-based accelerator concepts. IEEE Plasma Sci , 252–288 (1996).[5] Daido, H. Review of soft x-ray laser researches and develop-ments. Rep. Prog. Phys. , 1513–1576 (2002).[6] Kauranen, M. & Zayats, A. Nonlinear plasmonics. Nat. Pho-tonics , 737–748 (2012).[7] McKenna, P., Neely, D., Bingham, R. & Jaroszynski, D. Laser-Plasma Interactions and Applications (Springer, 2013).[8] Gibbon, P.
Short Pulse Laser Interactions with Matter: An In-troduction (Imperial College Press, 2005).[9] Abgrall, R. & Kumar, H. Robust Finite Volume Schemes forTwo-Fluid Plasma Equations.
J. Sci. Comput. , 584–611(2014).[10] Shumlak, U. & Loverich, J. Approximate Riemann solver forthe two-fluid plasma model. J. Comput. Phys. , 620–638(2003).[11] Roe, P. L. Approximate Riemann solvers, parameter vectors,and difference schemes.
J. Comput. Phys. , 357–372 (1981).[12] Loverich, J. & Shumlak, U. A discontinuous Galerkin methodfor the full two-fluid plasma model. Comput. Phys. Commun. , 251–255 (2005).[13] Loverich, J., Hakim, A. & Shumlak, U. A DiscontinuousGalerkin Method for Ideal Two-Fluid Plasma Equations.
Com-mun. Comput. Phys. , 240–268 (2011).[14] Srinivasan, B. & Shumlak, U. Analytical and computationalstudy of the ideal full two-fluid plasma model and asymptoticapproximations for Hall-magnetohydrodynamics. Phys. Plas-mas , 092–113 (2011).[15] Sousa, E. M. & Shumlak, U. A blended continuous-discontinuous finite element method for solving the multi-fluidplasma model. J. Comput. Phys. , 56–75 (2016).[16] Alvarez Laguna, A., Ozak, N., Lani, A., Deconinck, H. &Poedts, S. Fully-implicit finite volume method for the idealtwo-fluid plasma model.
Comput. Phys. Commun. , 31–44(2018).[17] Mason, R. J. An electromagnetic field algorithm for 2D implicitplasma simulation.
J. Comput. Phys. (1987).[18] Mason, R. J. & Cranfill, C. W. Hybrid Two-Dimensional Elec-tron Transport in Self-Consistent Electromagnetic Fields. IEEET. Plasma Sci. , 45–52 (1986).[19] Baboolal, S. & Bharuthram, R. Two-scale Numerical Solutionof the Electromagnetic Two-fluid plasma-Maxwell Equations:Shock and Soliton Simulation. Math. Comput. Simul. , 3–7(2007).[20] Kumar, H. & Mishra, S. Entropy Stable Numerical Schemesfor Two-Fluid Plasma Equations. J. Sci. Comput. , 401–425(2012).[21] Morel, B., Giust, R., Ardaneh, K. & Courvoisier, F. A simplesolver for the two-fluid plasma model based on PseudoSpectraltime-domain algorithm. Comm. Comput. Phys. , 955–978(2021). a) b) -7 -6 -5 -4 -3 -2 x (µm) -1-0.500.51 E z / E PSTD/Hydro
Propagation direction Interaction of artifactpre-pulses with plasmaPlasma Vacuum Artifactpre-pulses
PSATD/HydroAnalytical -7 -6 -5 -4 -3 -2 x (µm) -2-1012 v z ( m / s ) x 10 PSTD/HydroPSATD/Hydro
Plasma Vacuum
FIG. 6. (a) Single-cycle pulse propagation with PSTD/Hydro solver (red dashed-dotted line) and with PSATD/Hydro solver (solid blue line).The laser pulse interacts in normal incidence on a plasma ramp starting at x ≤ − µ m as shown in light blue on the figure. The PSATD/Hydrocurve is superimposed on the analytical solution shown as a black dashed line. We observe artifact pre-pulses generated by numerical dispersionof PSTD in the plasma region. (b) Velocity component v z of electron fluid at the same time. [22] Liu, Q. H. The PSTD algorithm: A time-domain method re-quiring only two cells per wavelength. Microw. Opt. Technol.Lett. , 158–165 (1997).[23] Liska, R. & Wendroff, B. Composite Schemes for ConservationLaws. SIAM J. Numer. Anal. , 2250–2271 (1998).[24] Vay, J.-L., Haber, I. & Godfrey, B. B. A domaindecomposition method for pseudo-spectral electromagneticsimulations of plasmas. J. Comput. Phys. , 260–268 (2013). URL .[25] Strang, G. On the Construction and Comparison of DifferenceSchemes.
SIAM J. Numer. Anal. , 506–517 (1968).[26] Hakim, A., Loverich, J. & Shumlak, U. A high resolution wavepropagation scheme for ideal Two-Fluid plasma equations. J.Comput. Phys. , 418–442 (2006).[27] Bittencourt, J. A.
Fundamentals of Plasma Physics (Springer,2013).[28] Yee, K. Numerical solution of initial boundary value problemsinvolving Maxwell’s equations in isotropic media.
IEEE Trans.Antennas Propag. , 302–307 (1966).[29] Shapoval, O., Vay, J.-L. & Vincenti, H. Two-step perfectlymatched layer for arbitrary-order pseudo-spectral analyticaltime-domain methods. Comput. Phys. Commun. , 102– 110 (2019). URL .[30] Richtmyer, D.
A Survey of Difference Methods for Non-SteadyFluid Dynamics (Technical Note NCAR/TN-63-2+STR, 1962).[31] Shampine, L. F. Two-step Lax–Friedrichs method.
Appl. Math.Lett. , 1134–1136 (2005).[32] Leveque, R. J. Finite Volume Methods for Hyperbolic Problems (Cambridge University Press, 2002).[33] Jackson, J. D.
Classical electrodynamics (John Wiley & Sons,1998). Third edition.[34] Speziale, T. & Catto, P. J. Linear wave conversion in an unmag-netized, collisionless plasma.
Phys. Fluids , 990–997 (1977).[35] Hinkel-Lipsker, D. E., Fried, B. D. & Morales, G. J. Analyticexpression for mode conversion of Langmuir and electromag-netic waves. Phys. Rev. Lett. , 2680–2682 (1989).[36] Forslund, D. W., Kindel, J. M., Lee, K., Lindman, E. L. &Morse, R. L. Theory and simulation of resonant absorption in ahot plasma. Phys. Rev. A , 679–683 (1975).[37] Mishra, S. K., Andreev, A. & Kalashinikov, M. P. Reflectionof few cycle laser pulses from an inhomogeneous overdenseplasma. Opt. Express25