Mathieu Drouin
email: mathieu.drouin@cea.fr

Laurent Gremillet
email: laurent.gremillet@cea.fr

Jean-Claude Adam

Anne Héron

J.-C Adam

Particle-in-cell modelling of relativistic laser-plasma interaction with the adjustable-damping, direct implicit method

Keywords: particle-in-cell method, implicit scheme, laser-plasma interaction, relativistic plasma

published or not. The documents may come

Particle-in-cell modelling of relativistic laser-plasma interaction with the adjustable damping, direct implicit method

Introduction

Particle-in-cell (PIC) codes have become widely used plasma simulation tools owing to their ability to mimic real plasma behavior. Yet the standard PIC algorithm employs an explicit time-differencing, and hence suffers from strict stability constraints on the time step, which needs to resolve the highest-frequency modes of the system [START_REF] Birdsall | Plasma physics via computer simulation[END_REF]. Furthermore, the mesh size must be comparable to the Debye length λ D in order to prevent the finite-grid instability [START_REF] Birdsall | Plasma physics via computer simulation[END_REF]. As a consequence, explicit PIC codes may find it difficult to cope with the large spatial and temporal scales associated with a number of physical scenarios, thus requiring massively parallel computing facilities [START_REF] Bowers | Ultrahigh performance three-dimensionnal electromagnetic relativistic kinetic plasma simulation[END_REF]. Several alternatives have been developed over the past decades to relax these constraints so that the choice of the space and time steps can be dictated by physical accuracy rather than stability conditions. The simplest way to do so is to suppress high-frequency processes within the mathematical model itself. Codes based on the Darwin-field approximation [START_REF] Hewett | Low-frequency electromagnetic (Darwin) applications in plasma simulations[END_REF][START_REF] Taguchi | Study of hot electron beam transport in high density plasma using 3D hybrid-Darwin code[END_REF], gyrokinetic equations [START_REF] Candy | An Eulerian gyrokinetic-Maxwell solver[END_REF] or hybrid particle-fluid models [START_REF] Mason | Monte Carlo hybrid modeling of electron transport in laser produced plasmas[END_REF][START_REF] Lipatov | The Hybrid Multiscale Simulation Technology. An Introduction with Application to Astrophysical and Laboratory Plasmas[END_REF][START_REF] Davies | Short-pulse high-intensity lasergenerated fast electron transport into thick solid targets[END_REF][START_REF] Gremillet | Filamented transport of lasergenerated relativistic electrons penetrating a solid target[END_REF][START_REF] Liljo | One-dimensional electromagnetic relativistic PIC-hydrodynamic hybrid simulation code H-VLPL (Hybrid Virtual Laser Plasma Lab)[END_REF] rely precisely on such an approach. The shortcoming inherent in these codes is the somewhat uncertain domain of validity of their basic assumptions. A second, more involved numerically, possibility retains a fully kinetic and electromagnetic description by using an implicit scheme for the entire Vlasov-Maxwell set of equations. This is the approach dealt with in this work.

The main feature, and difficulty, of a fully implicit PIC scheme is the prediction of the future particles' charge and current densities as functions of the future electromagnetic fields. Two main techniques have been designed to this goal. The first one to be published, the so-called moment method, makes use of the fluid equations to predict future source terms [START_REF] Denavit | Time-filtering particle simulations with ω pe δt ≫ 1[END_REF][START_REF] Mason | Implicit moment simulation of plasmas[END_REF][START_REF] Brackbill | An implicit method for electromagnetic plasma simulation in two dimensions[END_REF][START_REF] Mason | An electromagnetic field algorithm for 2d implicit plasma simulation[END_REF][START_REF] Vu | CELEST1D: an implicit, fully kinetic-model for low-frequency, electromagnetic plasma simulation[END_REF][START_REF] Lapenta | Kinetic approach to microscopicmacroscopic coupling in space and laboratory plasmas[END_REF]. and has been recently extended to the relativistic regime [START_REF] Noguchi | Formulation of the relativistic moment implicit particle-in-cell method[END_REF]. The present article will focus on the alternate approach, referred to as the direct implicit method, which is based on a direct linearization of the Lorentz equations [START_REF] Cohen | Implicit time integration for plasma simulation[END_REF][START_REF] Langdon | Direct implicit large time-step particle simulation of plasmas[END_REF][START_REF] Hewett | Electromagnetic direct implicit plasma simulation[END_REF][START_REF] Cohen | Performance and optimization of direct implicit particle simulation[END_REF]. Most implementations of the direct implicit method start with the so-called D 1 discretization of the Lorentz equation, first presented in Ref. [START_REF] Friedman | A direct method for implicit particle-in-cell simulation[END_REF]. The relativistic formulation, originally derived in Ref. [START_REF] Langdon | Relativistic extension of the electromagnetic direct implicit PIC algorithm[END_REF], was implemented, albeit in a simplified form, in the LSP code [START_REF] Welch | Implementation of a noniterative implicit electromagnetic field solver for dense plasma simulation[END_REF][START_REF] Welch | Implementation of a non-iterative implicit electromagnetic field solver for dense plasma simulation[END_REF][START_REF] Campbell | Collimation of PetaWatt laser-generated relativistic electron beams propagating though solid matter[END_REF][START_REF] Evans | Modelling short pulse, high intensity laser plasma interactions[END_REF][START_REF] Wei | Study of relativistic electron beam production and transport in high intensity laser interaction with a wire target by integrated LSP modeling[END_REF].

The direct implicit method proceeds as follows. First, particles' momenta and positions are advanced to an intermediate time level using known fields, yielding predicted charge and current densities. Second, by linearizing the latter quantities around the predicted momenta and positions, we can express correction terms as functions of the future fields and thus derive an implicit wave equation. Once this equation is solved, the particles' quantities are updated. Here we will show that the direct method can be derived as a simplified Newton scheme.

Our main motivation is the simulation of the interaction of an ultra-intense laser pulse with solid-density plasma slabs. The energetic particle beams originating from this interaction stir great interest in many fields spanning inertial confinement fusion [START_REF] Tabak | Ignition and high gain with ultrapowerful lasers[END_REF][START_REF] Campbell | Collimation of PetaWatt laser-generated relativistic electron beams propagating though solid matter[END_REF][START_REF] Mason | Heating mechanisms in short-pulse laser-driven cone targets[END_REF][START_REF] Honrubia | Three-dimensional fast electron transport for ignition-scale inertial fusion capsules[END_REF][START_REF] Atzeni | Fast ignitor target studies for the HiPER project[END_REF][START_REF] Chrisman | Intensity scaling in hot electron energy coupling in cone-guided fast ignition[END_REF], high energy density physics [START_REF] Patel | Isochoric heating of solid-density matter with an ultrafast proton beam[END_REF][START_REF] Martinolli | Fast-electron transport and heating of solid targets in high-intensity laser interactions measured by K α fluorescence[END_REF][START_REF] Akli | Laser heating of solid matter by light-pressure-driven shocks at ultrarelativistic intensities[END_REF][START_REF] Nilson | Bulk heating of solid-density plasmas during high-intensity-laser plasma interactions[END_REF], nuclear physics [START_REF] Cowan | Photonuclear fission from high energy electrons from ultraintense laser-solid interactions[END_REF][START_REF] Ledingham | Applications for nuclear phenomena generated by ultra-intense lasers[END_REF] or medical physics [START_REF] Bulanov | Feasibility of using laser ion accelerators in proton therapy[END_REF]. For the high plasma densities considered, the electron plasma frequency ω p largely exceeds the laser frequency.

Using an explicit PIC code, the space and time steps should resolve the highfrequency electron plasma modes of the plasma bulk. However, these modes are of no interest for the problem since they do not affect the laser-plasma interaction nor other potentially important related processes as the subsequent, fast electron-driven ion expansion. By contrast, resorting to an implicit scheme would allow a significantly increased time step, that is, determined only by the need to resolve the incoming laser wave. In this respect, one should realize that the strong wave damping inherent with implicit methods may be harmful in the context of laser-plasma interaction, for which light waves have to travel over many wavelengths. This prompted us to develop an electromagnetic solver with adjustable damping, based on a generalization of the scheme initially proposed by Friedman [START_REF] Friedman | A second-order implicit particle mover with adjustable damping[END_REF] for the Lorentz equation. We will demonstrate that our adjustable damping scheme tolerates abrupt spatial jumps in the controlling parameter. Our code therefore allows for dissipation-free laser propagation into vacuum, along with strong damping of undesirable plasma waves into the densest part of the target. As explicit codes, implicit codes suffer from the artificial heating arising from a crude discretization of the Debye length, as is commonplace when handling large-scale, high-density plasmas. This detrimental effect is generally attributed to the so-called grid-instability [START_REF] Birdsall | Plasma physics via computer simulation[END_REF]. To keep it at an acceptable level, we will exploit the well-known mitigating influence of high-order weight factors [START_REF] Abe | High-order spline interpolations in the particle simulation[END_REF][START_REF] Sentoku | Numerical methods for particle simulations at extreme densities and temperatures: Weighted particles, relativistic collisions and reduced currents[END_REF] by using quadratic weight factors. We will also take advantage of the stabilizing effect of the large time steps allowed by the implicit scheme.

The paper is organized as follows. In Sec. 2, we recall the basic principles of the PIC technique, give the implicit time-discretized equations to solve, and derive within a simplified Newton formalism the relativistic direct implicit method. In Sec. 3, we outline the numerical resolution of the wave equation as implemented in our newly developed, 2Dx-3Dv code ELIXIRS (ELectromagnetic Implicit X-dimensionnal Iterative Relativistic Solver). The introduction of implicit injecting/outgoing boundary conditions for the electromagnetic field is also discussed. Sec. 4 is devoted to the linear properties of the direct implicit method through the resolution of the electromagnetic and electrostatic dispersion relations. The effects of finite space and time steps, adjustable damping and high-order weight factors will be accounted for. Finally, in Sec. 5, our code is benchmarked against explicit simulations for two kinds of physical problems: the expansion of a plasma slab in vacuum, and the interaction of an ultra-intense laser pulse with an overcritical plasma target. The sensitivity of the simulation results to the damping parameter and the number of macro-particules will be addressed.

The relativistic direct implicit method as a simplified Newton scheme

In contrast to Ref. [START_REF] Langdon | Relativistic extension of the electromagnetic direct implicit PIC algorithm[END_REF], we present here a derivation of the electromagnetic direct implicit method for the relativistic case within a Newton iterative scheme and a weak formulation of Maxwell's equations. Anticipating our need of a dissipation-free propagation of light waves inside the vacuum region of the simulation domain, we introduce a generalization of the adjustable damping scheme proposed and used in the electrostatic regime by Friedman [START_REF] Friedman | A second-order implicit particle mover with adjustable damping[END_REF].

Basic equations

Consider Maxwell's equations

∇ × E = - ∂B ∂t , (1)
∇ × B = µ 0 j + 1 c 2 ∂E ∂t , (2)
and the collisionless Vlasov equation for the distribution function f s (x, u, t) of the sth particle species

∂f s ∂t + u γ ∂f s ∂x + q s m s E + u γ × B • ∂f s ∂u = 0 . (3)
Here q s and m s are the charge and the rest mass of the sth particle species, respectively. u denotes the relativistic momentum normalized by m s . The relativistic factor then writes γ = (1 + u 2 /c 2 ) 1/2 . The particle method consists in describing the distribution function f s as an ensemble of macro-particles in the form

f s (x, u, t) = Ns p=1 S(x -X p (t))δ(u -U p (t)) , (4)
where S is the shape function [START_REF] Birdsall | Plasma physics via computer simulation[END_REF], N s the total number of particles of the sth species, and δ the Dirac distribution. The relativistic motion of each macroparticle obeys the following equations:

dX p (t) dt = V p (t) = U p (t) γ p (t) , (5)
dU p (t) dt = q s m s E [X p (t), t] + U p (t) γ p (t) × B [X p (t), t] . (6)
We now make use of the implicit scheme with adjustable damping proposed by Friedman [START_REF] Friedman | A second-order implicit particle mover with adjustable damping[END_REF] for an electrostatic problem, which generalizes the so-called D 1scheme of Langdon et al. [START_REF] Cohen | Implicit time integration for plasma simulation[END_REF][START_REF] Langdon | Direct implicit large time-step particle simulation of plasmas[END_REF][START_REF] Hewett | Electromagnetic direct implicit plasma simulation[END_REF][START_REF] Langdon | Relativistic extension of the electromagnetic direct implicit PIC algorithm[END_REF]. The equations of motion are discretised as

X n+1 = X n + ∆t U n+1/2 γ n+1/2 , (7)
U n+1/2 = U n-1/2 + ∆t 2 (a n+1 + Ān-1 ) + q s ∆t 2m s U n+1/2 + U n-1/2 γ n × B n (X n ) , (8)
Ān-1 = θ f 2 a n + 1 - θ f 2 ān-2 , (9)
ān-1 = 1 - θ f 2 a n + θ f 2 ān-2 , (10)
where the index n denotes the time step index and we have defined

a n = q s m s E n , (11)
γ n = 1 + 1 c 2 U n-1/2 + ∆t 4 a n+1 + Ān-1 2 1/2 , (12)
γ n+1/2 = 1 + U 2 n+1/2 c 2 1/2 . (13)
Friedman's scheme can be readily applied to Maxwell's equations, which yields

E n+1 = E n + c 2 ∆t∇ × B n+1/2 - ∆t ǫ 0 j n+1/2 , (14)
B n+1/2 = B n-1/2 - ∆t 2 ∇ × E n+1 + Ēn-1 , (15)
B n = B n-1/2 - ∆t 2 ∇ × E n , (16)
Ēn-1 = θ f 2 E n + 1 - θ f 2 Ēn-2 , (17)
Ēn-1 = 1 - θ f 2 E n + θ f 2 Ēn-2 . ( 18
)
where j denotes the current density.

As will be demonstrated in Sec. 4, this scheme allows, via the parameter θ f , a flexible control of the damping of the high-frequency (electrostatic and electromagnetic) waves of the system. This property is of major interest for applications such as laser-plasma interaction involving a traveling electromagnetic wave into vacuum, for which the numerical damping associated with the standard D 1 method may prove too severe. It is worth noting that, even though referred to uniquely as θ f , the damping parameters involved in the electromagnetic scheme and the particle pusher may assume distinct values. The next sections will be devoted to the solution of the set of Eqs. ( 7)-( 18) within a Newton iterative scheme. We will show that for a proper choice of the initial conditions, this scheme reduces to the direct implicit method developed in Refs. [START_REF] Hewett | Electromagnetic direct implicit plasma simulation[END_REF][START_REF] Langdon | Relativistic extension of the electromagnetic direct implicit PIC algorithm[END_REF].

Weak formulation of the electric field equation

By replacing Eq. ( 15) into Eq. ( 14), one obtains the following wave equation

E n+1 + c 2 ∆t 2 2 ∇ × ∇ × E n+1 + ∆t ǫ 0 j n+1/2 = Q ′ , (19)
with the (known) source term

Q ′ = E n + c 2 ∆t∇ × B n-1/2 - c 2 ∆t 2 2 ∇ × ∇ × Ēn-1 . (20)
For any test function ψ, we assume the following weak formulation of the current density

j n+1/2 (x)ψ(x)dx = s q s 2 f s,0 (x, u)V n+1/2 (x, u) [ψ (X n+1 (x, u)) + ψ (X n (x, u))] dxdu , (21)
where f s,0 = f s (x, u, 0) is the initial particle distribution function and

V n+1/2 = U n+1/2 /γ n+1/2 .
The problem then consists in finding (E n+1 , X n+1 , U n+1/2 ) which solve

E n+1 (x)ψ(x)dx + c 2 ∆t 2 2 ∇ × ∇ × E n+1 (x)ψ(x)dx + ∆t ǫ 0 j n+1/2 (x)ψ(x)dx = Q ′ (x)ψ(x)dx (22)
together with Eqs. ( 7)- [START_REF] Brackbill | An implicit method for electromagnetic plasma simulation in two dimensions[END_REF]. We employ the Newton method to solve this system: for each quantity of interest Y , we introduce the ansatz

Y (k+1) n+α = Y (k) n+α + δY (k) n+α k = 0, 1, . . . ( 23
)
where α = (1/2, 1) depending on whether Y is centered at full or half time steps. The subscript n + 1 will be hereafter omitted for clarity. Substituting the above ansatz into Eq. ( 20) yields

E (k) (x) + δE (k) (x) ψ(x)dx + c 2 ∆t 2 2 ∇ × ∇ × E (k) (x) + δE (k) (x) ψ(x)dx + ∆t ǫ 0 j (k+1) (x)ψ(x)dx = Q ′ (x)ψ(x)dx . ( 24
)
The term involving j (k+1) is calculated with positions X (k+1) and velocities V (k+1)

j (k+1) ψ(x)dx = s q s 2 f s,0 (x, u)V (k) ψ(X (k) ) + ψ(X n ) dxdu + s q s 2 f s,0 (x, u)δV (k) ψ(X (k) ) + ψ(X n ) dxdu + s q s 2 f s,0 (x, u)V (k) ∇ψ(X (k) ) • δX (k) dxdu . ( 25
)
To obtain the equation solved for the electric field, we need to express the terms X (k) , δX (k) , V (k) and δV (k) as functions of the electric field. Before proceeding, let us first define the following quantities

γ (k) = 1 + U (k)2 c 2 1/2 , (26)
Γ (k) = 1 + 1 c 2 U n-1/2 + ∆t 4 q s m s E (k) (X (k) ) + Ān-1 2 1/2 , (27) θ
(X n ) = q s ∆t 2m s Γ (k) B n (X n ) , (28)
R(X n ) = 2 1 + θ 2 (I + θ ⊗ θ -θ × I) -I , (29)
M(U (k) ) = 1 γ (k) I - U (k) ⊗ U (k) γ (k)2 c 2 , ( 30
) N E (k) (X (k) ), U (k) = q s ∆t 4m s c 2 U n-1/2 + U (k) Γ (k)3 × B n (X n ) ⊗ U n-1/2 + ∆t 4 q s m s E (k) (X (k) ) + Ān-1 , (31)
with I the identity matrix. Straightforward calculations then yield

X (k) = X n + ∆tU (k) γ (k) , (32)
δX (k) = ∆tMδU (k) , (33)
V (k) = U (k) γ (k) , (34)
δV (k) = MδU (k) , (35)
Using the above expressions and the Newton ansatz [START_REF] Langdon | Relativistic extension of the electromagnetic direct implicit PIC algorithm[END_REF], the Lorentz equation becomes

U (k) + δU (k) = U n-1/2 + q s ∆t 2m s E (k) (X (k) ) + ∇E (k) (X (k) )δX (k) + δE (k) (X (k) ) + ∆t 2 Ān-1 + q s ∆t 2m s U (k) + δU (k) + U n-1/2 Γ (k) × B n (X n ) - q s ∆t 2m s N(E (k) (X (k) ), U (k) )∇E (k) (X (k) )δX (k) - q s ∆t 2m s N(E (k) (X (k) ), U (k) )δE (k) (X (k) ) , (36)
where we have dropped second-order terms. Assuming the electric field gradient term is negligible, this equation further simplifies as

U (k) + δU (k) = U n-1/2 + ∆t 4 [I + R(X n )] Ān-1 + q s m s E (k) (X (k) ) + ∆tq s 4m s [I + R(X n )] I -N E (k) (X (k) ), U (k) δE (k) (X (k) ) . ( 37
)
The set of equations ( 22)-( 38) constitutes the weak formulation of the problem.

We will now show how to recover the direct implicit method as a simplified Newton algorithm.

The direct implicit method

The simplest scheme consists in considering only one iteration in the above system and choosing the following initial values

   X (0) = X n+1 U (0) = U n+1/2 E (0) = 0    δX (0) = δX δU (0) = δU δE (0) = E (1) = E n+1 , (38)
where we have introduced the predicted position and momentum X n+1 and V n+ 1 2 computed from the known fields Ān-1 and B n . We have

X n+1 = X n + ∆t U n+1/2 γ n+1/2 , ( 39
)
U n+1/2 = R(X n )U n-1/2 + ∆t 4 [I + R(X n )] Ān-1 . ( 40
)
with γ n+1/2 = γ (0) . The correction terms then write

δU = q s ∆t 4m s [I + R(X n )][I -N( U n+1/2 )]E n+1 ( X n+1 ) , (41)
δV = MδU , ( 42
) δX = ∆tMδU , (43)
where we have defined

µ( U n+1/2 ) = N(0, U n+1/2 ) = q s ∆t 4m s c 2 U n-1/2 + U n+1/2 γ 3 n × B n (X n ) ⊗ U n-1/2 + ∆t 4 Ān-1 , (44)
and γ n = Γ (0) . After substituting the above equations into (25), using X n = X n+1 -∆t V n+1/2 and replacing the resulting expression into [START_REF] Welch | Implementation of a noniterative implicit electromagnetic field solver for dense plasma simulation[END_REF], we obtain

E n+1 (x)ψ(x)dx + c 2 ∆t 2 2 ∇ × ∇ × E n+1 (x)ψ(x)dx + s q s ∆t 2ǫ 0 f s,0 (x, u) V n+1/2 (x, u) ψ( X n+1 (x, u)) + ψ(X n (x, u)) dxdu + s q s ∆t ǫ 0 f s,0 (x, u)δV(x, u)ψ( X n+1 (x, u))dxdu + s q s ∆t 2ǫ 0 f s,0 (x, u) δX ⊗ V n+1/2 -V n+1/2 ⊗ δX ∇ψ( X n+1 )dxdu = Q ′ (x)ψ(x)dx . (45)
From Eq. ( 21), we identify

s q s ∆t 2ǫ 0 f s,0 (x, u) V n+1/2 (x, u) ψ( X n+1 (x, u)) + ψ(X n (x, u)) dxdu = ∆t ǫ 0 j n+1/2 (x)ψ(x)dx . ( 46
)
To reduce the next integral, it is convenient to introduce the weak formulation of the predicted charge density

ρ s (x)ψ(x)dx = q s f s,0 (x, u)ψ X n+1 (x, u) dxdu .
Approximating R(X n ) ≈ R( X n+1 ), we obtain

q s ∆t ǫ 0 f s,0 δVψ( X n+1 )dxdu = q s ∆t 2 4m s ǫ 0 ρ(x)M(x)(I + R(x)) [I -N(x)] E n+1 (x)dx . ( 47
)
Defining the implicit susceptibility χ as

χ(x) = s q s ∆t 2 4m s ǫ 0 M(x)(I + R s,n (x)) [I -N(x)] ρ s (x) , (48)
we have s q s ∆t ǫ 0 f s,0 (x, u)δV(x, u)ψ( X n+1 (x, u))dxdu = ψ(x)χ(x)E n+1 (x)dx .

(49) We treat the remaining integral by introducing the modified current j

+ s j + s (x)ψ(x)dx = q s f s,0 (x, u) V n+1/2 (x, u) ψ X n+1 (x, u) dxdu .
We then have

q s ∆t 2ǫ 0 f s,0 (x, u) δX ⊗ V n+1/2 -V n+1/2 ⊗ δX ∇ψ( X n+1 )dxdu = - q s ∆t 3 8m s ǫ 0 ∇ × j + s (x) × M(x) [I + R(x)] [I -N(x)] E n+1 (x) ψ(x)dx = - q s ∆t 3 8m s ǫ 0 ∇ × j + s (x) γ n+1/2 (x) × [I + R(x)] [I -N(x)] E n+1 (x) ψ(x)dx (50)
where use has been made of the identity U × U ⊗ U = 0. We are then led to define the tensor ζ as

ζ(x) = ∆t 2 8ǫ 0 s q s m s j + s γ n+1/2 × [I + R(x)] [I -N(x)] . (51)
There follows

q s ∆t 2ǫ 0 f s,0 δX ⊗ V n+1/2 -V n+1/2 ⊗ δX ∇ψdxdu = -∆t ∇ × (ζE n+1 )dx .
(52) Equation ( 25) supplemented by Eqs. ( 46), ( 49) and ( 52) should be satisfied for any test function ψ. As a result, we have to solve the local field equation

E n+1 + c 2 ∆t 2 2 ∇ × ∇ × E n+1 + χE n+1 -∆t∇ × ζE n+1 = Q , ( 53
)
where the source term now reads

Q = E n - ∆t ǫ 0 j n+1/2 + c 2 ∆t∇ × B n-1/2 - c 2 ∆t 2 2 ∇ × ∇ × Ēn-1 . ( 54
)
We have thus recovered the relativistic implicit method based on the D 1 scheme which was presented in Ref. [START_REF] Langdon | Relativistic extension of the electromagnetic direct implicit PIC algorithm[END_REF], with the only difference that the source term now involves the time-averaged field Ēn-1 . It then appears that the direct implicit method can be derived as a one-iteration Newton method with the starting values X (0) = X n+1 , U (0) = U n+1/2 and E (0) = 0.

Numerical resolution

Resolution of the field equation

In this section, we sketch the numerical procedure used to solve Eq. ( 53) in the case of a 2Dx-3Dv phase space with periodic boundary conditions along the transverse y axis. We have first to evaluate the implicit susceptibilities. These terms are computed for each macroparticle, yielding χ(X p , U p ) and ζ(X p , U p ), before being projected onto the (x, y) grid through the usual formulas:

χ(x) = s p S(X p -x)χ(X p , U p ) , (55)
ζ(x) = s p S(X p -x)ζ(X p , U p ) . ( 56
)
We then apply the iterative method of Concus and Golub [START_REF] Concus | Use of fast direct methods for the efficient numerical solution of nonseparable elliptic equations[END_REF] to solve the elliptic system defined by Eq. ( 53), which reads in the present case

E (m+1) + c 2 ∆t 2 2 ∇ × ∇ × E (m+1) + χE (m+1) -∆t∇ × ζE (m+1) = Q (m) (57)
The right-hand side of Eq. ( 57) is given by

Q (m) = Q -(χ -χ 0 )E (m) + ∆t∇ × ζ -ζ 0 E (m) ( 58
)
where m is the iteration index and χ 0 and ζ 0 denote the y-averaged susceptibilities. The fast convergence of the scheme implies, in principle, slow variations of the field quantities in the y direction, but this has not proved particularly constraining for the physical situations we have considered.

As is usual in electromagnetic PIC codes, two interleaved meshes are used for the spatial differencing of the grid quantities. The fields are discretized as follows: ρ i,j , J z,i,j , E z,i,j , J x,i+1/2,j , E x,i+1/2,j , B y,i+1/2,j , J y,i,j+1/2 , E y,i,j+1/2 , B x,i,j+1/2 and B z,i+1/2,j+1/2 . The χ and ζ are stored at (i, j) except for χ 11 , ζ 11 , ζ 21 , ζ 31 , which are located at (i + 1/2, j), and χ 22 , ζ 12 , ζ 22 , ζ 32 , located at (i, j + 1/2). Once space-discretized, the above equations are Fourier transformed along the y direction. Considering N y grid cells, we obtain N y one-dimensional equations to solve. Considering N x grid cells in the x direction, each equation gives a 6N x system of equations. These systems have a band-diagonal structure and are solved by a standard LU technique, using routines bandec and banbks of the numerical recipes library [START_REF] Press | Numerical Recipes in Fortran 90: The Art of Scientific Computing[END_REF]. Details on spatial discretisations and Fourier transformations used to solve Eq. (57) are given in Appendix A.

Charge correction

Our method to accumulate charge and current densities [Eqs. [START_REF] Cohen | Performance and optimization of direct implicit particle simulation[END_REF] and [START_REF] Esirkepov | Exact charge conservation scheme for particle-in-cell simulation with an arbitrary form-factor[END_REF]] does not satisfy charge conservation, which results into the violation of Poisson's equation. This is a common flaw of early electromagnetic PIC codes [START_REF] Birdsall | Plasma physics via computer simulation[END_REF] which may be corrected by a more sophisticated projection scheme [START_REF] Villasenor | Rigorous charge conservation for local electromagnetic field solvers[END_REF][START_REF] Esirkepov | Exact charge conservation scheme for particle-in-cell simulation with an arbitrary form-factor[END_REF]. A wellknown alternative approach, which will be implemented here, is to correct the electrostatic part of the electric field E n+1 solution of Eq. ( 53) so that it fulfills Poisson's equation [START_REF] Birdsall | Plasma physics via computer simulation[END_REF]. Using normalized quantities, our best statement of Gauss's law is

∇ • E * n+1 = ρ n+1 , (59)
where E * n+1 represents the sought-for electric field. Using

ρ n+1 = ρ n+1 -∇ • χE * n+1 , this can be reformulated as ∇ • (1 + χ)E * n+1 = ρ n+1 . (60)
Now, taking the divergence of Eq. ( 53) yields

∇ • [(1 + χ)E n+1 ] = ∇ • Q (61) with generally ∇ • Q = ρ n+1 . We may first think of introducing a potential ψ such that Q * = Q -∇ψ fulfills ∇ • Q * = ρ n+1
, but this correction has been shown to cause spurious effects [START_REF] Hewett | Electromagnetic direct implicit plasma simulation[END_REF]. A proper correction makes use of the following form [START_REF] Hewett | Electromagnetic direct implicit plasma simulation[END_REF]]

Q * = Q -(I + χ)∇ψ , (62)
There follows

∇ • [(1 + χ)∇ψ] = ∇ • Q -ρ n+1 , (63)
which is equivalent to

∇ • [(1 + χ)∇ψ] = ∇ • [(1 + χ)E n+1 ] -ρ n+1 , (64)
where the only unknown is the scalar field ψ. Eventually, the corrected field E * n+1 ensuring Eq. ( 60) is given by E * n+1 = E n+1 -∇ψ. Details on the numerical resolution of Eq. ( 64) are given in Appendix B.

Electromagnetic boundary conditions

In this section we describe the implementation of injecting/outgoing boundary conditions on both sides of the simulation box. Incident and scattered electromagnetic waves are assumed linearly polarized and depending on the phase term k • x -ωt only. Waves polarized in the (x, y) plane then verify

E inc y = B inc z cos θ i , (65)
E scat y = -B scat z cos θ s , (66)
where θ i and θ s denote respectively the incident and scattered angles. The total field becomes

E tot y = E scat y + E inc y (67) = -B tot z cos θ s + E inc y cos θ i (cos θ i + cos θ s ) (68)
Discretizing with centered finite differences in space and time gives

1 4 E n+1 y,1,j+1/2 + E n+1 y,0,j+1/2 + E n y,1,j+1/2 + E n y,0,j+1/2 = -B n+1/2 z,1/2,j+1/2 cos θ s + E inc,n+1/2 y,1/2,j+1/2 (cos θ i + cos θ s ) cos θ i . ( 69
)
Using Maxwell-Faraday's equation, we can express E n+1 y,0,j+1/2 as a function of the field values at inner grid points and previous time steps. We have

E n+1 y,0,j+1/2 = AE n+1 y,1,j+1/2 2 ∆t ∆x cos θ s -1 - 2A∆t ∆y cos θ s E n+1 x,1/2,j+1 -E n+1 x,1/2,j -4A cos θ s B n-1/2 z,1/2,j+1/2 + 2A∆t ∆x cos θ s Ēn-1 y,1,j+1/2 -Ēn-1 y,0,j+1/2 - 2A∆t ∆y cos θ s Ēn-1 x,1/2,j+1 -Ēn-1 x,1/2,j + 4A cos θ i (cos θ i + cos θ s ) E inc,n+1/2 y,1/2,j+1/2 -A E n y,1,j+1/2 + E n y,0,j+1/2 , (70)
where the coefficient A is given by

A = 1 + 2 ∆t ∆x cos θ s -1 . (71)
A similar equation can be established for z-polarized waves, which reads

E n+1 z,0,j = BE n+1 z,1,j 2∆t ∆x cos θ s -1 -B(E n z,0,j + E n z,1,j ) + 4B cos θ s B n-1/2 y,1/2,j + 2B∆t ∆x cos θ s Ēn-1 z,1,j -Ēn-1 z,0,j + 4BE inc,n+1/2 z,1/2,j 1 + cos θ s cos θ i , (72)
where we have defined the coefficient B as

B = 1 + 2∆t ∆x cos θ s -1 . ( 73
)
Note that the above equations only apply in vacuum. This is realized in practice by imposing boundary conditions on particles a few grid cells away from the outer boundaries of the computational domain.

Numerical analysis of the adjustable-damping, direct implicit method

Dispersion relation of electromagnetic waves in vacuum

Our aim here is to quantify the error in phase velocity and the damping associated with electromagnetic waves as functions of the space and time steps.

In particular, we will demonstrate the possibility to control the wave damping by adjusting the parameter θ f .

Combining Maxwell-Ampère's ( 14) and Maxwell-Faraday's (15) equations and assuming propagation in vacuum yield the wave equation

E n+1 = 2E n -E n-1 - c 2 ∆t 2 2 ∇ × ∇ × E n+1 + Ēn-1 . (74)
The time-filtered term involves the adjustable damping parameter θ f [Eq. ( 17)] and can be expanded as

E n+1 + Ēn-1 = E n+1 + θ f 2 E n + 1 - θ f 2 2 E n-1 + 1 - θ f 2 2 θ f 2 E n-2 + 1 - θ f 2 2 θ f 2 2 E n-3 + . . . (75)
In a 2-D geometry, taking the electric field in the form E n = E 0 Φ(x, y)z n with z = exp(-iω∆t) and i = √ -1, Eq. ( 75) becomes

E n+1 + Ēn-1 = E 0 Φ(x, y) z -1 1 - θ f 2 2 + θ f 2 z + z 2 + 1 - θ f 2 2 θ f 2 z -2 1 + θ f 2 z -1 + θ f 2 2 z -2 + . . . z n . ( 76
)
where the adjustable damping parameter θ f ∈ [0, 1]. Simplifying the series in the right-hand side of Eq. ( 76) yields

E n+1 + Ēn-1 = E 0 Φ(x, y) z -1 1 - θ f 2 2 + θ f 2 z + z 2 + 1 - θ f 2 2 θ f 2 2z -1 2z -θ f z n . ( 77
)
The electromagnetic wave is assumed polarized in the (x, y) plane with a harmonic dependence Φ(x, y) = exp [i(k x x + k y y)]. Substituting Eq. (77) into Eq.( 74) and space-differencing the Laplacian leads, we get after some straightforward algrebra the following third degree polynomial equation

z 2 = 2z -1 - 1 - θ f 2 2 + θ f 2 z + z 2 + 1 - θ f 2 2 θ f 2z -θ f Ω 2 2 , (78)
where we have introduced

Ω 2 = 4 c 2 ∆t 2 ∆x 2 sin 2 k x ∆x 2 + c 2 ∆t 2 ∆y 2 sin 2 k y ∆y 2 . ( 79
)
Equation ( 78) simplifies as

z 3 (2 + Ω 2 ) -z 2 (4 + θ f ) + z 2 + Ω 2 (1 -θ f ) + 2θ f -θ f = 0 . ( 80
)
Let us first examine the special case θ f = 0. The roots of interest are solutions of

z 2 (2 + Ω 2 ) -4z + (2 + Ω 2 ) = 0 (81)
The discriminant ∆ = 4 -(2 + Ω 2 ) 2 being always negative, we get the roots

z ± = (2 ± i √ -∆)/(2 + Ω 2 ), which statisfy |z + | = |z -| = 1.
We have therefore demonstrated the absence of damping when θ f = 0. Figure 1 plots the normalized phase velocity

v φ = ℜω kc (where k = k 2 x + k 2 y )
for different values of c∆t/∆x = c∆t/∆y. The phase velocity error grows for increasing ∆x and ∆t/∆x. A value c∆t/∆x > 1, that is, violating the stability constraint of the standard explicit scheme, therefore implies a moderate spatial step k x ∆x 0.38 (c∆t/∆x = 1.27) so as to avoid excessive (> 5%) phase velocity error, which, in presence of relativistic particles, may cause unphysical Cerenkov radiation [START_REF] Greenwood | On the elimination of numerical Cerenkov radiation in PIC simulations[END_REF].

Let us now address the case of nonzero θ f . Figures 2 and3 plot the normalized phase velocity v φ /c (left) and damping rate ℑω∆t (right) of the least damped root of Eq. (80) as functions of (k x ∆x, k y ∆y) for θ f = 1. Cuts of these two quantities in the plane k y = 0 are represented in Figures 4 and5 respectively. Again the phase velocity error grows for increasing ∆x and ∆t/∆x. A value c∆t/∆x > 1, therefore implies a reduced spatial step k x ∆x 0.28 (c∆t/∆x = 1.27) so as to keep phase velocity error below 5%. In this case the damping rate, which also increases with ∆x and ∆t/∆x, proves much too strong for applications relying on the propagation of an electromagnetic wave over several wavelengths. For example, assuming k x ∆x = 0.2 and c∆t/∆x = 1, a typical travel time of 200∆t requires |ℑω|∆t < 2.5 × 10 -4 for a tolerable wave dissipation (< 5%). As seen in Fig. 5(right), this condition cannot be fulfilled when θ f = 1, which further demonstrates the need for an adjustable-damping scheme for a proper modeling of laser-plasma interaction.

Dispersion relation of electrostatic plasma waves

We will now focus on the numerical relation dispersion of the electron plasma fluctuations in the case of a uniform, nonrelativistic Maxwellian plasma with a fixed neutralizing background. For this purpose, we shall adopt the formalism of Langdon [START_REF] Langdon | Analysis of the time integration in plasma simulation[END_REF] that accounts for both finite space and time steps, as well as allows for an arbitrary time-differencing scheme of the Lorentz equation. An infinite number of macroparticles is assumed, yielding a continuous velocity distribution function (taken in the Maxwellian form). In this framework, as detailed in Appendix C, the present adjustable-damping, direct implicit algorithm can be easily managed. The relation dispersion yielding the complex frequency ω as a function of the wave number k then reads

1 + (∆x/λ D ) 2 (k∆x) 2 sin(k∆x/2) k∆x/2 2 +∞ p=-∞ sin (k p ∆x/2) k p ∆x/2 2m+2 sin(k p ∆x) k p ∆x +∞ q=-∞ [1 + ξ q Z(ξ q )] + (ω p ∆t) 2 /2 (k∆x) 2 sin(k∆x/2) k∆x/2 2 +∞ p=-∞ (k p ∆x) 2 sin (k p ∆x/2) k p ∆x/2 2m+2 sin(k p ∆x) k p ∆x S(θ f ) = 0 , ( 82
)
where m is the order of the shape factor [START_REF] Birdsall | Plasma physics via computer simulation[END_REF]. k p = k -2πp/∆x and ω q = ω -2πq/∆t are the aliased wave number and frequency, respectively. Z denotes the plasma dispersion function [START_REF] Fried | The Plasma Dispersion Function[END_REF] whose argument is ξ q = ω q / √ 2k p v t (where v t is the electron thermal velocity). Moreover, we have defined the function S as

S(θ f ) = +∞ s=0 e i(ω/ωp)s(ωp∆t) (2/θ f ) s e -1 2 (λ D /∆x) 2 s 2 (k∆x) 2 (ωp∆t) 2 , ( 83
)
with the value S(0) = 1. We have numerically solved Eq. ( 82) using the nonlinear solver STRSCNE developed in Ref. [START_REF] Bellavia | STRSCNE: A scaled trust-region solver for constrained nonlinear equations[END_REF] and the algorithm of Ref. [START_REF] Weideman | Computation of the complex error function[END_REF] to compute the Z function. We will restrict the following analysis to systems characterized by a crude resolution of the Debye length (∆x/λ D > 1), as is commonplace in simulations of large-scale, high-density plasmas. Figure 6 displays the k-dependence of the complex frequency of the fastest growing (or least damped) mode solution of Eq. (82) for θ f = 1, ω p ∆t = 2 and various values of ∆x/λ D . For ∆x/λ D = 32 (i.e., v t ∆t/∆x = 0.06), most of the k-spectrum is damped except for a bounded unstable region located near k∆x ∼ 2.6 with a maximum growth rate ℑω/ω p ∼ 0.011. This corresponds to the well-known finite-grid instability [START_REF] Birdsall | Plasma physics via computer simulation[END_REF] commonly afflicting PIC simulations with ∆x/λ D ≫ 1, and responsible for nonphysical field energy growth and plasma heating. This instability originates from the interplay of the aliased wave numbers in Eq. (82). Note also the nonphysical k-dependence of the real frequency obtained at large ω p ∆t : ℜω is significantly below ω p at k = 0 and further drops with increasing k∆x. As seen in Fig. 6, decreasing ∆x/λ D eventually leads to a complete stabilization of the system along with a displacement of the dominant mode towards low k values. For ∆x/λ D = 4 (i.e., v t ∆t/∆x = 0.5), the least damped mode is thus located at k∆x = 0.76 with ℑω/ω p ∼ -0.1. This evolution points to a transition between spatial step-dominated and time-step-dominated regimes.

The dependence of the characteristics of the dominant mode on the ratio ∆x/λ D ≫ 1 and the weight factor order is summarized in Table 1 for θ f = 1 and ω p ∆t = 2. The benefit of a high-order interpolation scheme is clearly evidenced: the system turns out to be entirely stabilized up to ∆x/λ D = 32 with a quadratic weight factor, and ∆x/λ D = 64 with a cubic weight factor. In addition, the wavenumber of the increasingly damped dominant mode is shifted downward.

A connection between the present calculations and previously published simulation results [START_REF] Brackbill | An implicit method for electromagnetic plasma simulation in two dimensions[END_REF][START_REF] Cohen | Performance and optimization of direct implicit particle simulation[END_REF] is provided by Tables 2 and3, which display the dependence of the dominant mode on the ratio v t ∆t/∆x = ω p ∆t/(∆x/λ D ), as well as on the damping parameter (the time step being fixed to ω p ∆t = 2). An extensive set of implicit electrostatic PIC simulations using the D 1 scheme (i.e., θ f = 1) and linear interpolation has indeed revealed that satisfactory energy conservation can be achieved in the range [START_REF] Brackbill | An implicit method for electromagnetic plasma simulation in two dimensions[END_REF][START_REF] Cohen | Performance and optimization of direct implicit particle simulation[END_REF] 0.

1 v t ∆t ∆x 1 (84)
Even though the present stability analysis alone is not expected to account for the complex issue of numerical self-heating [START_REF] Birdsall | Plasma physics via computer simulation[END_REF][START_REF] Ueda | A study of the numerical heating in electrostatic particle simulations[END_REF], the results of Table 2 are found in reasonable agreement with the lower bound of the above heuristic range, as they indicate a complete stabilization of the system for v t ∆t/∆x 0.1 in case of a linear weigth factor and θ f = 1. For lower θ f values, stabilization is reached for increased v t ∆t/∆x. Moreover, Table 3 shows that the use of a quadratic weight factor permits to suppress the finite-grid instability at reduced v t ∆t/∆x ( 0.06 for θ f = 1). Similarly to Fig. 6,a regime to the low-k temporal regime is evidenced when raising v t ∆t/∆x. As expected, a high-order (m > 1) weight factor, which enables to filter out high spatial frequencies, proves beneficial only in the high-k, grid-instability regime (for v t ∆t/∆x 0.25). Note that we have not considered values v t ∆t/∆x > 1 since, in the present case, this would imply ∆x/λ D < 2, a parameter range of little practical interest for the aforementioned applications.

Further insight into the stability properties of the adjustable-damping scheme is given by fixing the ratio v t ∆t/∆x = 0.09 and varying accordingly the space and time steps. Equivalently, within the laser-plasma context which we propose to address, this can be achieved by fixing the parameters ω 0 ∆x/c and ω 0 ∆t (where ω 0 is the incident laser frequency) and varying the plasma density. The resulting data is displayed in Table 4 in the ranges 1.26 ≤ ω p t ≤ 8.94 and 14.3 ≤ ∆x/λ D ≤ 101.1. One can see that a linear shape factor proves rather inappropriate for most of the parameter range considered. By contrast, complete stabilization is achieved for n ≥ 2 weight factors. It is worth noting that, in terms of laser-plasma parameters, the rightmost column of Table 4 corresponds to a 2000n c , 1 keV plasma (where n c is the critical density at the laser frequency ω 0 ) discretized with ω 0 ∆t = 0.2 and ω 0 ∆x/c = 0.1. In addition to accessing such extreme plasma conditions, employing a cubic weight factor may give the opportunity to reduce the damping parameter θ f .

θ f 0 0.1 0.

Numerical applications

Wave propagation in vacuum

Here, we illustrate the capability of the adjustable damping, implicit scheme implemented in the code ELIXIRS to manage the propagation of electromagnetic waves in vacuum. Let us consider a plane wave, with normalized vector potential a 0 = 3 and frequency ω 0 , entering the left-hand side of a 1024∆x × 4∆y box, with ∆x = 0.2c/ω 0 , ∆y = 0.8c/ω 0 and ∆t = 0.2ω -1 0 . The wave is injected and absorbed using the procedure detailed in 3.3. Figure 7(left) shows the expected monotonous damping of the incident wave induced when a spatially uniform damping parameter θ f = 1 is applied. After propagating across the simulation box, the wave amplitude is measured to be 46% of the initial value, which is close to the theoretical value (49%). The opposite, dissipation-free case corresponding to θ f = 0 is displayed in Fig. 7(right). Finally, with the problem of laser plasma interaction in mind, we address the case of a spatially varying θ f profile in the form   

θ f = 0, 0 < ω 0 x/c < 51.2 θ f = 1, 51.2 < ω 0 x/c < 153.6 θ f = 0, 153.6 < ω 0 x/c < 204.8 (85)
Figure 7(center) shows that the discontinuity in θ f does not cause significant spurious effects. This sought-for property is of major interest for modeling laser-plasma interaction as it allows the laser wave to travel unperturbed in vacuum over several wavelengths before reaching the overcritical target, whose numerical stability calls for finite numerical damping. For the sake of completeness, we have checked that the weak (∼ 0.1% in the present case) re-flection arising at the discontinuity surface is consistent with Fresnel's formula R = (N(1) -N(0)) 2 /(N(1) + N(0)) 2 , where N(θ f ) = c/v φ (θ f ) is the numerical refraction index derived in Sec. 4.1.

Plasma expansion into vacuum: benchmarking against explicit simulations

As a first test of the implicit Vlasov-Maxwell solver, we simulate the dynamics of a plasma slab freely expanding into vacuum. The results of the implicit code ELIXIR are confronted to refined, explicit simulations performed with the code CALDER [START_REF] Lefebvre | Electron and photon production from relativistic laser-plasma interactions[END_REF]. We consider a 0.6c/ω p plasma slab composed of hot (10 keV) electrons and cold ions. In the implicit case, the simulation box is 103∆x × 4∆y large, with ∆x = 2c/ω p and ∆y = 0.4c/ω p (yielding the ratios ∆x/λ D = 14 and v t ∆t/∆x = 0.14), whereas the explicit simulation handles a 1024∆x × 8∆y box, with ∆x = ∆y = 0.2c/ω p . A linear weight factor is used in all cases. Figures 8, 9 and 10 plot the time evolution of the ion density profile, the ion phase space and the time evolution of the plasma kinetic energies, as simulated by the implicit and explicit codes. The implicit damping parameter is chosen to be θ f = 1, whereas the total number of macroparticles N p is 6 × 10 4 and 6 × 10 5 in the implicit and explicit cases, respectively. Overall, albeit roughly resolved and strongly damped (as expected from Table 1), the implicit scheme manages to satisfactorily capture the finely resolved, explicit results. Yet, the wave damping gives rise to artificial electron cooling, which results into a weakened ion acceleration as seen in Figs. 9 and 10. More quantitatively, the total energy drops by ∼ 3%, yielding a maximum ion energy of ∼ 160 keV, as compared to ∼ 220 keV in the explicit case. For the sake of completeness, we have carried out additional calculations so as to assess the influence of the damping parameter and the number of macroparticules. For each simulation, we have measured the energy variation and the peak ion energy. The data thus obtained is summarized in Tables 5 and6. The implicit scheme behaves reasonably well up to θ f = 0.15 with an energy variation < 10%, comparable or better than its explicit counterpart for an equal number of macroparticles. Increasing the latter from 6 × 10 4 to 6 × 10 5 approximately halves the energy variation but hardly changes the peak ion energy. The transition from numerical electron cooling and heating occurs between θ f = 1 and θ f = 0.5. Finally, the undamped (θ f =0) case is subject to a much stronger, if still limited, electron heating, which translates into a twofold overestimate of the peak ion energy.

A parametric study of plasma self-heating and cooling

We have carried out a series of simulations of the free evolution of an electronion plasma to gauge the potential discrepancy between the idealized linear analysis of Sec. 4.2 and the actual predictor-corrector numerical scheme. Evidently, the objective is to gain further insight into the energy conservation properties of the latter and the predictive capability of the former. These calculations draw upon and extend the work of Ref. [START_REF] Cohen | Performance and optimization of direct implicit particle simulation[END_REF] to the electromagnetic regime. The system consists of a bounded electron-ion plasma with T e = T i = 1 keV and m i /m e = 900, extending over half a 300∆x × 4∆y simulation box. We have scanned the (∆x/λ D , ω p ∆t) parameter space in the range [5, 60] × [START_REF] Birdsall | Plasma physics via computer simulation[END_REF][START_REF] Candy | An Eulerian gyrokinetic-Maxwell solver[END_REF]. In practice, after introducing ω 0 , the frequency of a fictitious electromagnetic wave, and n c , the corresponding critical density, we have set ∆x = 0.2c/ω 0 and varied the ratio n e /n c and the time step so that ∆x/λ D ∈ {5, 10, 20, 30, 60} and ω p ∆t ∈ {1, 2, 5}. The damping parameter is θ f = 1. The total simulation time is kept fixed at 1000ω -1 0 . For each simulation, we have calculated the relative variation of the total kinetic energy per time step (∆K/K 0 )/N (where ∆K is the kinetic variation, K 0 the initial kinetic energy and N the number of time steps). To be complete, we have also performed electrostatic calculations, whereby the electric field is directly computed through the Poisson equation (64). 3.2 × 10 -4 9.9 × 10 -4 1.7 × 10 -3 2.8 × 10 -3 2 -1.1 × 10 -4 1.3 × 10 -4 6.9 × 10 -4 1.2 × 10 -3 2.6 × 10 -3 5 -5.8 × 10 -5 -2.4 × 10 -4 2.9 × 10 -5 3 × 10 -4 9.5 × 10 -4 The results are summarized in Tables 789. The associated plots of the kinetic energies are shown in Figs. 111213: each column corresponds to a specific value of ∆x/λ D and each line to a specific value of ω p ∆t. Note that we have excluded in these plots the case ∆x/λ D = 60 as it always gives rise to significant numerical heating. We have checked that the plasma kinetic energy makes up for most of the system energy. Overall, the electrostatic results prove close to the electromagnetic ones. Satisfactory energy conservation ( 10 -4 ) is obtained for v t ∆t/∆x 0.2 and v t ∆t/∆x 0.1 in the linear and quadratic interpolation cases, respectively. These lower bound values are in fairly good agreement, albeit slightly higher, with the linear results of Sec. 4.2. Larger v t ∆t/∆x ratios eventually lead to plasma cooling,

High intensity laser interaction with an overdense plasma slab 5.4.1. Quasi-one-dimensional simulation

Let us now address the problem of the interaction of a relativistic-intensity laser pulse with an overcritical plasma, which is the prime motivation behind this work.

As a first illustration, we consider the case of a quasi-1D laser-plasma system. The irradiated target consists of a 60c/ω 0 -long, 1 keV, 200n c plasma slab preceded by a 18c/ω 0 -long density ramp rising linearly from 0 to 200n c . The incident electromagnetic plane wave has a 120ω -1 0 constant-intensity profile with a 22ω -1 0 rise time and a normalized amplitude a 0 = eE 0 /m e cω 0 = 3. The implicit simulation employs a 2048∆x × 4∆y grid, with ∆x = ∆y = 0.1c/ω 0 and ∆t = 0.14ω -1 0 , yielding, in terms of plasma parameters, ∆x/λ D = 32 and ω p ∆t = 2 (v t ∆t/∆x = 0.06). The damping parameter in the electromagnetic solver, as well as in the particle pusher, is set to zero in the vacuum region and the moderately dense plasma region up to n e = 60n c , and to unity in the denser plasma region. Guided by the results of Sec. 5.3, we make use of a quadratic weight factor to reduce the numerical heating. The number of macroparticles per cell N p is varied from 100 to 1300. These calculations are compared with explicit simulations using the same parameters except for a decreased time step ∆t = 0.05ω -1 0 so as to fulfill the Courant stability condition.

Explicit Implicit (θ

f = 0) Implicit (θ f = 1 if n e > 60n c ) N p = 1300 +14.4% +6% -3% N p = 400 +15.3% +10.5% -1% N p = 100 +22% +25.5% +12.7%
Table 10: Quasi-1D laser-plasma interaction: energy variation in the explicit simulations with ∆t = 0.05ω -1 0 and the implicit simulations with ∆t = 0.14ω -1 0 and varying θ f . See text for other simulation parameters.

Table 10 compares the values of the total energy variation (calculated after complete reflection of the laser pulse) as obtained in the explicit and implicit cases. Results from implicit simulations with zero damping are also displayed. Overall, except for N p = 100, for which case the three schemes behave similarly, the implicit simulations are found to achieve better energy conservation than their explicit counterparts. The benefit of a strongly damped scheme in the densest region of the plasma is mostly evidenced for N p = 1300 and 400. The not-so-good performances of the explicit calculations prompted us to carry out an additional, more refined explicit simulation that can serve more properly as a reference calculation. This simulation made use of a 4096∆x × 8∆y grid with ∆x = ∆y = 0.05c/ω 0 and ∆t = 0.03ω -1 0 , as well as of a third-order weight factor with N p = 650. It yielded a total energy variation of 4%.

The electron (x, p x ) phase space (integrated in the y-direction) is displayed in Fig. 14 for both explicit and implicit schemes. Consistently with the wellknown ponderomotive heating mechanism arising at relativistic laser intensities, fast electrons are accelerated into the target as bunches separated by half the laser wavelength [START_REF] Adam | Dispersion and transport of energetic particles created during the interaction of intense laser pulses with overdense plasmas[END_REF]. The explicit simulation predicts maximum electron momenta about 20% higher than that predicted by the implicit simulation. Also, as a result of the damping of longitudinal beam-plasma modes, the implicit simulation exhibits a longer-lived separation between the thermal electrons and the fast electrons as the latter propagate through the target. In an actual solid-density configuration, though, the beam-plasma wave mixing observed in the explicit case should be suppressed by collisions as demonstrated in Ref. [START_REF] Kemp | Collisional relaxation of superthermal electrons generated in relativitsic laser pulses in dense plasmas[END_REF]. Yet, these discrepancies do not translate into major differences in the electron energy dis- tribution as shown at three successive times in Fig. 15. In particular, the slope of the high-energy tail of the spectra is satisfactorily reproduced. The reduced electron heating gives rise in turn to a ∼ 15% slower, space-charge-driven ion acceleration into vacuum as depicted by the ion (x, p x ) phase spaces of Fig. 16.

Two-dimensional simulations

We now consider a fully two-dimensional laser-plasma system. The electronion plasma slab has a peak density of 200n c , a temperature of 1 keV and a thickness of 6c/ω 0 . A 12c/ω 0 -long linear density ramp is added in front of the target. The simulation box consists of a 1024 × 512 grid with ∆x = ∆y = 0.1c/ω 0 (∆x/λ D = 32). The incoming laser pulse has unchanged parameters except for a 12c/ω 0 FWHM Gaussian transverse profile. Open and periodic boundary conditions are applied for the electromagnetic fields along the xand y-axis, respectively. Due to memory constraints, we use a rather small number of macroparticles N p = 40. So as to stabilize the system, in addition to using a quadratic weight factor, the time step is significantly increased as compared to the previous simulations: ∆t = 0.3ω -1 0 , which corresponds to ω p ∆t = 4.2 and v t ∆t/∆x = 0.13. Particles are subject to periodic boundary conditions in the y-direction, and reinjected with their initial temperature in the x-direction. The damping parameter in the electromagnetic solver, as well as in the particle pusher, is set to zero in the vacuum region and the moderately dense plasma region up to n e = 30n c . Two maximum values of the spatially varying damping parameter have been tried in the denser plasma region: θ f = 0.1 and 0.5. The explicit simulation of reference makes use of a third-order weight factor with the parameters ∆x = ∆y = 0.08c/ω 0 , ∆t = 0.05ω -1 0 and N p = 160. This parallel calculation takes 4.5h on 64 1.6 GHz Itanium 2 processors. By contrast, the (sequential) implicit simulations take 27h on a 2.66 GHz Intel Xeon X5355 processor.

The time evolution of the particle kinetic energies is displayed in Fig. 18. All simulations predict about the same peak electron energy. Yet, the damped implicit calculations yield a faster decreasing electron energy. The total energy variation, evaluated over the time interval 215 < ω 0 t < 715 (that is, after complete reflection of the laser pulse and before the fastest ions hit the box boundaries) is -12% and -15% for the θ f = 0.1 and θ f = 0.5 implicit cases, respectively, as compared to +5% in the explicit case. Despite their crude time resolution and limited number of macroparticles, the implicit calculations manage to reproduce quite accurately the salient features of the fast electron and ion generation. This is evidenced by the electron and ion (x, p x ) phase spaces of Figs. 19 and20, as well as by the electron energy spectra of Fig. 22. As in the previous Section, if to a lesser extent due to the weaker numerical damping employed here, the implicit simulations somewhat underestimate the maximum electron energies. A 2-D picture of the fast electron generation is provided by the map of the electron kinetic energy density shown in Fig. 23. A reasonable agreement is observed between the three cases, each calculation showing the characteristic 2ω 0 -bunched propagation of the fast electrons and their breakout into vacuum.

Conclusion

This paper has been devoted to the application of the relativistic direct implicit method to the problem of laser-plasma interaction. In contrast to closely related works [START_REF] Campbell | Collimation of PetaWatt laser-generated relativistic electron beams propagating though solid matter[END_REF][START_REF] Evans | Modelling short pulse, high intensity laser plasma interactions[END_REF][START_REF] Wei | Study of relativistic electron beam production and transport in high intensity laser interaction with a wire target by integrated LSP modeling[END_REF], our scheme, implemented inside the 2Dx-3Dv code ELIXIRS, allows for high-order weight functions and adjustable damping of the high-frequency waves. The latter capability, which extends to electromagnetic waves a method originally designed by Friedman [START_REF] Friedman | A second-order implicit particle mover with adjustable damping[END_REF] for electrostatic waves, permits to manage within a unified algorithm the dissipation-free, Courant conditionfree propagation of the incident laser pulse through vacuum, while suppressing the need to resolve the high-frequency collective modes inside the dense plasma region. After having presented an original derivation of the adjustable-damping, direct implicit method as a simplified, one-iteration Newton scheme, we have carried out a thorough analysis of its numerical properties regarding both electromagnetic and electrostatic waves. The latter study, accounting for the effects of finite ∆t and ∆x, the weight factor order and the damping parameter is found to provide useful hints when compared to the simulation results of the free evolution of a plasma slab. Several numerical tests have been presented and successfuly benchmarked against finely resolved explicit simulations. In particular, we have demonstrated the ability of the code to capture the main features of the laserplasma interaction despite cruder space-time resolution. Yet, our code being still sequential, its increased stability domain remains insufficient to access the large space-and time-scales managed nowadays by massively parallel explicit codes. The parallelization of our code is therefore required and will be the subject of a future work.

Acknowledgments

We gratefully acknowledged the work of U. Voss on the application of the direct implicit method to the problem of laser-plasma interaction. This study, which provided us with important guidelines, was carried out in 1998 at the CMAP/ École Polytechnique and supported by the EU TMR grant FMBICT972082.

A. Numerical implementation of the field equation

We detail here the numerical procedure to solve Eq. ( 57) within a 2D geometry. The Concus and Golub iterative method [START_REF] Concus | Use of fast direct methods for the efficient numerical solution of nonseparable elliptic equations[END_REF] is applied to the three components of Eq. ( 57). The x-component writes

E n+1 x,i+1/2,j + c 2 ∆t 2 2∆x∆y E n+1 y,i+1,j+1/2 -E n+1 y,i+1,j-1/2 -E n+1 y,i,j+1/2 + E n+1 y,i,j-1/2 - c 2 ∆t 2 2∆y 2 E n+1 x,i+1/2,j+1 -2E n+1 x,i+1/2,j + E n+1 x,i+1/2,j-1 + χ 11,0 i+1/2 E n+1 x,i+1/2,j + 1 4 χ 12,0 i E n+1 y,i,j+1/2 + χ 12,0 i E n+1 y,i,j-1/2 +χ 12,0 i+1 E n+1 y,i+1,j-1/2 + χ 12,0 i+1 E n+1 y,i+1,j+1/2 + 1 2 χ 13,0 i E n+1 z,i,j + 1 2 χ 13,0 i+1 E n+1 z,i+1,j - ∆t 2∆y ζ 31,0 i+1/2 E n+1 x,i+1/2,j+1 -ζ 31,0 i+1/2 E n+1 x,i+1/2,j-1 - ∆t 2∆y ζ 32,0 i E n+1 y,i,j+1/2 + ζ 32,0 i+1 E n+1 y,i+1,j+1/2 -ζ 32,0 i E n+1 y,i,j-1/2 -ζ 32,0 i+1 E n+1 y,i+1,j-1/2 - ∆t 4∆y ζ 33,0 i+1 E n+1 z,i+1,j+1 + ζ 33,0 i E n+1 z,i,j+1 -ζ 33,0 i+1 E n+1 z,i+1,j-1 -ζ 33,0 i E n+1 z,i,j-1 = Q x,i+1/2,j . (86)
The y-component writes

E n+1 y,i,j+1/2 - c 2 ∆t 2 2∆x 2 E n+1 y,i+1,j+1/2 -2E n+1 y,i,j+1/2 + E n+1 y,i-1,j+1/2 + c 2 ∆t 2 2∆x∆y E n+1 x,i+1/2,j+1 -E n+1 x,i-1/2,j+1 -E n+1 x,i+1/2,j + E n+1 x,i-1/2,j + χ 21,0 i 4 E n+1 x,i-1/2,j + E n+1 x,i+1/2,j + E n+1 x,i-1/2,j+1 + E n+1 x,i+1/2,j+1 + χ 22,0 i E n+1 y,i,j+1/2 + χ 23,0 i 2 (E n+1 z,i,j + E n+1 z,i,j+1 ) + ∆t 2∆x ζ 31,0 i+1/2 (E n+1 x,i+1/2,j + E n+1 x,i+1/2,j+1 ) -ζ 31,0 i-1/2 (E n+1 x,i-1/2,j + E n+1 x,i-1/2,j+1 ) + ∆t 2∆x ζ 32,0 i+1 E n+1 y,i+1,j+1/2 -ζ 32,0 i-1 E n+1 y,i-1,j+1/2 + ∆t 4∆x ζ 33,0 i+1 (E n+1 z,i+1,j + E n+1 z,i+1,j+1 ) -ζ 33,0 i-1 (E n+1 z,i-1,j + E n+1 z,i-1,j+1 ) = Q y,i,j+1/2 . ( 87
)
The z-component writes

E n+1 z,i,j - c 2 ∆t 2 2∆x 2 E n+1 z,i+1,j -2E n+1 z,i,j + E n+1 z,i-1,j - c 2 ∆t 2 2∆y 2 E n+1 z,i,j+1 -2E n+1 z,i,j + E n+1 z,i,j-1 + χ 31,0 i 2 E n+1 x,i-1/2,j + E n+1 x,i+1/2,j + χ 32,0 i 2 E n+1 y,i,j-1/2 + E n+1 y,i,j+1/2 + χ 33,0 i E n+1 z,i,j - ∆t ∆x ζ 21,0 i+1/2 E n+1 x,i+1/2,j -ζ 21,0 i-1/2 E n+1 x,i-1/2,j - ∆t 4∆x ζ 22,0 i+1 E n+1 y,i+1,j-1/2 + E n+1 y,i+1,j+1/2 -ζ 22,0 i-1 E n+1 y,i-1,j+1/2 + E n+1 y,i-1,j-1/2 - ∆t 2∆x ζ 23,0 i+1 E n+1 z,i+1,j -ζ 23,0 i-1 E n+1 z,i-1,j + ∆t 4∆y ζ 11,0 i E n+1 x,i+1/2,j+1 + E n+1 x,i-1/2,j+1 -E n+1 x,i+1/2,j-1 -E n+1
x,i-1/2,j-1

+ ∆t ∆y ζ 12,0 i E n+1 y,i,j+1/2 -E n+1 y,i,j-1/2 + ∆t 2∆y ζ 13,0 i E n+1 z,i,j+1 -E n+1 z,i,j-1 = Q z,i,j . (88)
The right-hand sides of Eqs. ( 86)-( 88) are given by -

Q(m) x,i+1/2,j = Q x,i+1/2,j -(χ 11 i+1/2,j -χ 11,0 i+1/2,j )E (m) x,i+1/2,
-ζ 32 i,j-1/2 -ζ 32,0 i E (m) y,i,j-1/2 -ζ 32 i+1,j-1/2 -ζ 32,0 i+1 E (m) y,i+1,j-1/2 + ∆t 4∆y ζ 33 i+1,j+1 -ζ 33,0 i+1 E (m) z,i+1,j+1 + ζ 33 i,j+1 -ζ 33,0 i E (m) z,i,j+1 -ζ 33 i+1,j-1 -ζ 33,0 i+1 E (m) z,i+1,j-1 -ζ 33 i,j-1 -ζ 33,0 i E (m) z,i,j-1 , (89)
ζ 31 i-1/2,j -ζ 31,0 i-1/2 E (m)
x,i-1/2,j -

ζ 31 i-1/2,j+1 -ζ 31,0 i-1/2 E (m) x,i-1/2,j+1 - ∆t 2∆x ζ 32 i+1,j+1/2 -ζ 32,0 i+1 E (m) y,i+1,j+1/2 -ζ 32 i-1,j+1/2 -ζ 32,0 i-1 E (m) y,i-1,j+1/2 - ∆t 4∆x ζ 33 i+1,j -ζ 33,0 i+1 E (m) z,i+1,j + ζ 33 i+1,j+1 -ζ 33,0 i+1 E (m) z,i+1,j+1 -ζ 33 i-1,j -ζ 33,0 i-1 E (m) z,i-1,j -ζ 33 i-1,j+1 -ζ 33,0 i-1 E (m) z,i-1,j+1 , Q(m) z,i,j =Q z,i,j - 1 2 (χ 31 i,j -χ 31,0 i ) E (m) x,i-1/2,j + E (m) x,i+1/2,j - 1 2 χ 32 i,j -χ 32,0 i E (m) y,i,j-1/2 + E (m) y,i,j+1/2 -χ 33 i,j -χ 33,0 i E (m) z,i,j + ∆t ∆x ζ 21 i+1/2,j -ζ 21,0 i+1/2 E (m)
x,i+1/2,j -

ζ 21 i-1/2,j -ζ 21,0 i-1/2 E (m) x,i-1/2,j + ∆t 4∆x ζ 22 i+1,j-1/2 -ζ 22,0 i+1 E (m) y,i+1,j-1/2 + ζ 22 i+1,j+1/2 -ζ 22,0 i+1 E (m) y,i+1,j+1/2 -ζ 22 i-1,j-1/2 -ζ 22,0 i-1 E (m) y,i-1,j-1/2 -ζ 22 i-1,j+1/2 -ζ 22,0 i-1 E (m) y,i-1,j+1/2 + ∆t 2∆x ζ 23 i+1,j -ζ 23,0 i+1 E (m) z,i+1,j -ζ 23 i-1,j -ζ 23,0 i-1 E (m) z,i-1,j - ∆t 4∆y ζ 11 i+1/2,j+1 -ζ 11,0 i+1/2 E (m) x,i+1/2,j+1 + ζ 11 i-1/2,j+1 -ζ 11,0 i-1/2 E (m) x,i-1/2,j+1 -ζ 11 i+1/2,j-1 -ζ 11,0 i+1/2 E (m)
x,i+1/2,j-1 -

ζ 11 i-1/2,j-1 -ζ 11,0 i-1/2 E (m) x,i-1/2,j-1 - ∆t ∆y ζ 12 i,j+1/2 -ζ 12,0 i E (m) y,i,j+1/2 -ζ 12 i,j-1/2 -ζ 12,0 i E (m) y,i,j-1/2 - ∆t 2∆y ζ 13 i,j+1 -ζ 13,0 i E (m) z,i,j+1 -ζ 13 i,j-1 -ζ 13,0 i E (m) z,i,j-1 .
Assuming periodicity of the electric field along the y direction, we Fourier transform Eqs. ( 86)-(88) in this direction. We introduce E R k and E I k the real and imaginary parts of the Fourier transformed electric field. For notational simplicity, the index k will be omitted in the following. The real part of the Fourier transform of Eq. ( 86

+ E R z i+1 χ 13,0 i+1 2 + E I z i+1 ∆t 2∆y ζ 33,0 i+1 sin( k∆y) = Q R x i+1/2 . (90)
The imaginary part of the Fourier transform of Eq. ( 86) reads

+ E R z i+1 - ∆t 2∆y ζ 33,0 i+1 sin( k∆y) + E I z i+1 χ 13,0 i+1 2 = Q I x i+1/2 . ( 91
) E R y i-1 - c 2 ∆t 2 2∆x 2 - ∆t 2∆x ζ 32,0 i-1 + E R z i-1 - ∆t 4∆x ζ 33,0 i-1 cos( k∆y) + 1 + E I z i-1 - ∆t 4∆x ζ 33,0 i-1 sin( k∆y) + E R x i-1/2 c 2 ∆t 2 2∆x∆y cos( k∆y) -1 + χ 21,0 i 4 cos( k∆y) + 1 - ∆t 2∆x ζ 31,0 i-1/2 cos( k∆y) + 1 + E I x i-1/2 c 2 ∆t 2 2∆x∆y + χ 21,0 i 4 - ∆t 2∆x ζ 31,0 i-1/2 sin( k∆y) + E R y i 1 + c 2 ∆t 2 ∆x 2 + χ 22,0 i + E R z i χ 23,0 i 2 cos( k∆y) + 1 + E I z i χ 23,0 i 2 sin( k∆y) + E R x i+1/2 - c 2 ∆t 2 2∆x∆y cos( k∆y) -1 + χ 21,0 i 4 cos( k∆y) + 1 + ∆t 2∆x ζ 31,0 i+1/2 cos( k∆y) + 1 + E I x i+1/2 - c 2 ∆t 2 2∆x∆y + χ 21,0 i 4 + ∆t 2∆x ζ 31,0 i+1/2 sin( k∆y) + E R y i+1 - c 2 ∆t 2 2∆x 2 + ∆t 2∆x ζ 32,0 i+1 + E R z i+1 ∆t 4∆x ζ 33,0 i+1 cos( k∆y) + 1 + E I z i+1 ∆t 4∆x ζ 33,0 i+1 sin( k∆y) = Q R y i . (92)
The imaginary part of the Fourier transform of Eq. (87) reads

E I y i-1 - c 2 ∆t 2 2∆x 2 - ∆t 2∆x ζ 32,0 i-1 + E R z i-1 ∆t 4∆x ζ 33,0 i-1 sin( k∆y) + E I z i-1 - ∆t 4∆x ζ 33,0 i-1 cos( k∆y) + 1 + E R x i-1/2 -
+ E I y i+1 - c 2 ∆t 2 2∆x 2 + ∆t 2∆x ζ 32,0 i+1 + E R z i+1 - ∆t 4∆x ζ 33,0 i+1 sin( k∆y) + E I z i+1 ∆t 4∆x ζ 33,0 i+1 cos( k∆y) + 1 = Q I y i . (93)
The real part of the Fourier transform of Eq. ( 88) reads

E R y i-1 ∆t 4∆x ζ 22,0 i-1 cos( k∆y) + 1 + E I y i-1 - ∆t 4∆x ζ 22,0 i-1 sin( k∆y) + E R z i-1 - c 2 ∆t 2 2∆x 2 + ∆t 2∆x ζ 23,0 i-1 + E R x i-1/2 χ 31,0 i 2 + ∆t ∆x ζ 21,0 i-1/2 + E I x i-1/
+ E R z i 1 + c 2 ∆t 2 ∆x 2 + c 2 ∆t 2 ∆y 2 1 -cos( k∆y) + χ 33,0 i + E I z i - ∆t ∆y ζ 13,0 i sin( k∆y) + E R x i+1/2 χ 31,0 i 2 - ∆t ∆x ζ 21,0 i+1/2 + E I x i+1/2 - ∆t 2∆y ζ 11,0 i sin( k∆y) + E R y i+1 - ∆t 4∆x ζ 22,0 i+1 cos( k∆y) + 1 + E I y i+1 ∆t 4∆x ζ 22,0 i+1 sin( k∆y) + E R z i+1 - c 2 ∆t 2 2∆x 2 - ∆t 2∆x ζ 23,0 i+1 = Q R z i . ( 94
)
The imaginary part of the Fourier transform of Eq. ( 88) reads

E R y i-1 ∆t 4∆x ζ 22,0 i-1 sin( k∆y) + E I y i-1 ∆t 4∆x ζ 22,0 i-1 cos( k∆y) + 1 + E I z i-1 - c 2 ∆t 2 2∆x 2 + ∆t 2∆x ζ 23,0 i-1 + E R x i-1/2
+ E I z i 1 + c 2 ∆t 2 ∆x 2 + c 2 ∆t 2 ∆y 2 1 -cos( k∆y) + χ 33,0 i + E R x i+1/2 ∆t 2∆y ζ 11,0 i sin( k∆y) + E I x i+1/2 χ 31,0 i 2 - ∆t ∆x ζ 21,0 i+1/2 + E R y i+1 - ∆t 4∆x ζ 22,0 i+1 sin( k∆y) + E I y i+1 - ∆t 4∆x ζ 22,0 i+1 cos( k∆y) + 1 + E I z i+1 - c 2 ∆t 2 2∆x 2 - ∆t 2∆x ζ 23,0 i+1 = Q I z i . (95)
Considering N x grid points along x-direction Eqs. ( 90)-(95) can be formulated as a band-diagonal system of equations, which we solve using a LU technique [START_REF] Press | Numerical Recipes in Fortran 90: The Art of Scientific Computing[END_REF] for each of the N y modes of the discrete Fourier transform. Then we compute the field solution in real space by inverse Fourier transformation.

B. Numerical implementation of the charge correction step

We detail here the numerical procedure to solve Eq. ( 64) within a 2D geometry. As for the wave equation, we make use of the Concus and Golub iterative method [START_REF] Concus | Use of fast direct methods for the efficient numerical solution of nonseparable elliptic equations[END_REF], which writes in the present case

-∇ • (1 + χ 0 )∇ψ (m+1) = ρ -∇ • [(1 + χ)E n+1 ] + ∇ • (χ -χ 0 )∇ψ (m) (96)
where χ 0 = χ kl,0 1≤k,l≤3 denotes the y-averaged χ susceptibility tensor with χ kl,0 =< χ kl > y . E n+1 is the solution of the wave equation ( 53) and m denotes the iteration index. Omitting the latter, we discretize the above equation in the form

- 1 ∆x 1 + χ 11,0 i+1/2,j 1 ∆x (ψ i+1,j -ψ i,j ) -1 + χ 11,0 i-1/2,j 1 ∆x (ψ i,j -ψ i-1,j ) - 1 2∆x χ 12,0 i+1,j 1 2∆y (ψ i+1,j+1 -ψ i+1,j-1 ) -χ 12,0 i-1,j 1 2∆y (ψ i-1,j+1 -ψ i-1,j-1 ) - 1 2∆y χ 21,0 i,j+1 1 2∆x
(ψ i+1,j+1 -ψ i-1,j+1 ) -χ 21,0 i,j-1

1 2∆x
(ψ i+1,j-1 -ψ i-1,j-1 )

-
+ ρ i,j (99)
The above equations are Fourier transformed along the y direction. Considering N y grid cells we have to solve N y one-dimensional equations. Assuming N x grid cells in the x direction, each equation turns out into a 2N x system of equations. These systems have a band-diagonal structure and are solved with a LU technique [START_REF] Press | Numerical Recipes in Fortran 90: The Art of Scientific Computing[END_REF].

C. Derivation of the dispersion relation of electron plasma waves with finite ∆x and ∆t

We restrict our analysis to a one-dimensional, nonrelativistic electrostatic plasma with immobile ions. In the following, we adopt the methodology and notations of Ref. [START_REF] Birdsall | Plasma physics via computer simulation[END_REF]. For a single macro-particle, the adjustable-damping scheme ( 7)- [START_REF] Liljo | One-dimensional electromagnetic relativistic PIC-hydrodynamic hybrid simulation code H-VLPL (Hybrid Virtual Laser Plasma Lab)[END_REF]

where n stands for the time step index. We now assume a harmonic form for the interpolated electric force F (1) = F (k)e i(kx-ωt) . As a direct consequence of the PIC interpolation scheme, we have the relation [START_REF] Birdsall | Plasma physics via computer simulation[END_REF] F (k) = qE(k)S(-k)

where E(k) and S(k) are the discrete Fourier transforms of the electric field and the m-order weight function, respectively. The latter reads

S(k) = sin (k∆x/2) k∆x/2 m+1 . (102)
The first-order acceleration term can then be expressed as

a n = F (k) m exp i(kx (0) n -ωt n ) = F (k) m exp ik(x 0 + v (0) t n ) -iωt n = F (k) m exp ikx 0 exp [i(kv -ω)n∆t] . (103)
Defining A(k) = F (k) m e ikx 0 and z = e i(kv-ω)∆t , Eq. ( 103) reads

x n+1 -2x n + x n-1 = ∆t 2 2 A(k) z n+1 + 1 2 z n + 1 2 2 z n-1 + 1 2 3 z n-2 + . . . x n+1 -2x n + x n-1 = ∆t 2 2 A(k)z n z -1 1 - θ f 2 2 + θ f 2 z + z 2 + 1 - θ f 2 2 θ f 2 z -2 1 + θ f 2 z -1 + θ f 2 2 z -2 + . . . . (104)
This equation can be further simplified as

x n+1 -2x n + x n-1 = ∆t 2 2 A(k)2z n [(1 -θ f ) + z 2 ] 2z -θ f . (105)
We linearize x n = x

Where the polynomial P reads

P(k) = 2z n [(1 -θ f ) + z 2 ] 2z -θ f (107) (108)
We deduce that x

n (x 0 , v 0 , t n ) varies as e i(kv-ω)n∆t = z n . Hence we find the solution

x (1) n = ∆t 2 m F (k)e i(kx-ωt) z (z -1) 2 + z 2z -θ f (109)
To evaluate the charge density, we introduce the dipole density P (x, t) = n 0 q dvf 0 (v)x (1) n (x, v, t)

= -n 0 q m F (k)e i(kx-ωt) dvf 0 (v) 1

Figure 1 :Figure 2 :

 12 Figure 1: Phase velocity of the least damped root of Eq. (80) as a function of (k x ∆x, k y ∆y), for different values of c∆t/∆x = c∆t/∆y ∈ {0.05, 0.66, 1.28, 1.9, 2.5} (from top to bottom) and θ f = 0. A narrower (k x ∆x, k y ∆y) range is represented on the right.

Figure 3 :

 3 Figure 3: Same as Fig. 2 but with a narrower (k x ∆x, k y ∆y) range.

Figure 4 :

 4 Figure 4: Phase velocity (left) and damping rate ℑω∆t (right) of the least damped root of Eq. (80) as a function of (k x ∆x), for different values of c∆t/∆x = c∆t/∆y ∈ {0.05, 0.66, 1.28, 1.9, 2.5} (from top to bottom on the left and bottom to top on the right) and θ f = 1. Phase velocity without damping (θ f = 0) is represented by dotted-dashed line.

Figure 5 :

 5 Figure 5: Same as Fig. 4 but with a narrower (k x ∆x) range.

Figure 6 :

 6 Figure 6: Real frequency (blue) and growth rate (red) vs k∆x of the dominant mode solving Eq. (82) with ω p ∆t = 2, θ f = 1 and a linear weight factor (n = 1): ∆x/λ D = 32 (left), 20 (center) and 4 (right).

Figure 7 :

 7 Figure 7: Propagation of a plane wave with θ f = 1 (top, left), θ f = 0 (top, right), and a spatially varying θ f profile according to Eq. (85) (bottom).

Figure 8 :

 8 Figure 8: Time evolution of the ion density profile: explicit (left) and implicit (right) simulations with ∆x = 0.2c/ω p , ∆t = 0.1ω -1 p , N p = 6 × 10 5 and ∆x = 2c/ω p , ∆t = 2ω -1 p , N p = 6 × 10 4 , respectively. The implicit damping parameter is θ f = 1.

Figure 9 :

 9 Figure 9: Ion phase space at t = 2600ω -1 p : explicit (left) and implicit (right) simulations with ∆x = 0.2c/ω p , ∆t = 0.1ω -1 p , N p = 6 × 10 5 and ∆x = 2c/ω p , ∆t = 2ω -1 p , N p = 6 × 10 4 , respectively. The implicit damping parameter is θ f = 1.

Figure 10 :

 10 Figure 10: Time evolution of the electron (red) and ion (green) kinetic energies: explicit (left) and implicit (right) simulations with ∆x = 0.2c/ω p , ∆t = 0.1ω -1 p , N p = 6 × 10 5 and ∆x = 2c/ω p , ∆t = 2ω -1 p , N p = 6 × 10 4 , respectively. The implicit damping parameter is θ f = 1.

2 × 2 - 9 . 2 × 0 - 1 . 6 ×

 2292016 10 -4 1.1 × 10 -3 2.1 × 10 -3 4.7 × 10 -3 10 -5 1.5 × 10 -4 7.9 × 10 -4 1.5 × 10 -3 4.5 × 10 -3 5 10 -4 1.2 × 10 -4 4.7 × 10 -4 1.7 × 10 -3

Figure 11 :

 11 Figure 11: Time evolution of the total (blue), ion (red) and electron (green) energies: electrostatic case with linear weight factor. ∆x/λ D = (5, 10, 20, 30) from left to right and ω p ∆t = (1, 2, 5) from top to bottom.

Figure 12 :

 12 Figure 12: Time evolution of the total (blue), ion (red) and electron (green) energies: electromagnetic case with linear weight factor. ∆x/λ D = (5, 10, 20, 30) from left to right and ω p ∆t = (1, 2, 5) from top to bottom.

Figure 13 :

 13 Figure 13: Time evolution of the total (blue), ion (red) and electron (green) energies: electromagnetic case with quadratic weight factor.

Figure 14 :

 14 Figure 14: Electron (x, p x ) phase space at t = 198ω -1 0 : explicit simulation (left) and implicit simulation with θ f = 1 (right). In both cases, N p = 1300. See text for other simulation parameters.

Figure 15 :

 15 Figure 15: Electron energy distribution at different times: explicit simulation (red) and implicit simulation (blue). Energy is normalized by m e c 2 .

Figure 16 :

 16 Figure 16: Ion (x, p x ) phase space at t = 792ω -1 0 : explicit simulation (left) and implicit simulation with θ f = 1 (right). In both cases, N p = 1300. See text for other simulation parameters.

Figure 17 :

 17 Figure 17: Time evolution of the electron (red) and ion (green) kinetic energies: explicit simulation (left) and implicit simulation with θ f = 1 (right). In both cases, N p = 1300. See text for other simulation parameters.

Figure 18 :

 18 Figure 18: Time evolution of the electron (red) and ion (green) kinetic energies: explicit simulation (left), implicit simulations with θ f = 0.1 (center) and θ f = 0.5 (right).

Figure 19 :

 19 Figure 19: Electron (x, p x ) phase space at t = 96ω -1 0 : explicit simulation (left) and implicit simulations with θ f = 0.1 (center) and θ f = 0.5 (right).

Figure 20 :

 20 Figure 20: Ion (x, p x ) phase space at t = 523ω -1 0 : explicit simulation (left) and implicit simulations with θ f = 0.1 (center) and θ f = 0.5 (right).

Figure 21 :

 21 Figure 21: Electron energy distribution at different times: explicit simulation (red) and implicit simulation with θ f = 0.1 (blue). Energy is normalized by m e c 2 .

Figure 22 :

 22 Figure 22: Electron energy distribution at different times: explicit simulation (red) and implicit simulation with θ f = 0.5 (blue). Energy is normalized by m e c 2 .

Figure 23 :

 23 Figure 23: Electron kinetic energy density (normalized by m e c 2 n c ) at t = 67ω -1 0 and t = 86ω -1 0 : explicit simulation (top) and implicit simulations with θ f = 0.1 (center) and θ f = 0.5 (bottom).

2 ∆t sin(ω -kv) ∆t 2 2 + n 0 q∆t 2 2mF 1 vt √ 2π exp -v √ 2vt 2 ,θ f ) s e -v 2 t 2 (sk∆t) 2 (F 2 t 2 (sk∆t) 2 ( 112 ) 1 ω 2 t 2 ( 2 p κ 2 2 p κ 2 (k) ∆t 2 2 p 1 2 v 2 t

 222122222221121222222212 (k)e i(kx-ωt) dvf 0 (v) ∞ s=0 e i(ω-kv)s∆t (2/θ f ) s(110)The first and second terms of the right-hand side correspond to the explicit leapfrog scheme and the implicit correction, respectively. Assuming a Maxwellian distribution f 0 (v) = the latter can be writtendvf 0 (v) ∞ s=0 e i(ω-kv)s∆t (2/θ f ) s = ∞ s=0 e iωs∆t (2/θ f ) s dvf 0 (v)e -ikvs∆t = ∞ s=0 e iωs∆t (2/θ f ) s F (f 0 )(ks∆t) 111)where F denotes the Fourier transform. Thus the polarisation becomesP (x, t) = n 0 q m F (k)e i(kx-ωt) (k)e i(kx-ωt) ∞ s=0 e iωs∆t (2/θ f ) s e -vWe can develop cotan as a series in the formcotan (ω -kv) -kv -qω g (113)The continuous charge density is given by ρ p = -∇ • P, which writes in Fourier space ρ p (k) = -ikP (k). The discrete charge density is then given byρ(k) = p S(k p )ρ p (k p ) = -i p k p S(k p )P (k p ) = -i p |S(k p )| 2 n 0 q 2 m E(k p ) +∞ q=-∞ dv ∂f 0 (v)/∂v ω -k p v -qω g -i p k p |S(k p )| 2 n 0 q 2 ∆t 2 2m E(k p ) ∞ s=0 e iωs∆t (2/θ) s e -v sk∆t) 2 . (114)Using centered space-differencing, discrete Fourier transform of the relationE = -∂φ/∂x gives E(k) = -iK(k)φ(k) = -iK(k)φ(k), as modified by the direct implicit method reads∇ • (∇φ n+1 ) = -ρ n+1 ǫ 0(117)Centered space-differencing followed by a Fourier transformation givesκ 2 (k)φ(k) = ρ(114)-(119), we obtain the dispersion relation for an infinite electrostatic one dimensional plasma taking into account both spatial and temporal discretizations ǫ(ω, k) = 1 + ω (k) p |S(k p )| 2 K(k p ) +∞ q=-∞ dv ∂f 0 (v)/∂v ω -k p v -qω g + ω k p |S(k p )| 2 K(k p ) +∞ s=0 e iωs∆t (2/θ) s e -(sk∆t) 2 = 0 , (120)

Table 1 :

 1 Imaginary frequency ℑω/ω p (wavenumber k∆x) of the dominant mode as a function of the ratio ∆x/λ D and the weight factor order for ω p ∆t = 2 and θ f = 1.

	∆x/λ D	14.3	22.6	32	64
	linear	-0.024 3.3 × 10 -3 (2.11) (2.42)	0.011 (2.58)	0.01 (2.85)
	quadratic -0.04 (1.96)	-0.015 (2.30)	-3.7 × 10 -3 2.8 × 10 -3 (2.48) (2.70)
	cubic	-0.039	-0.018

-8.6 × 10 -3 -2 × 10 -4 (1.84) (2.14) (2.36) (2.67)

  clear transition from the high-k spatial

	θ f	0	0.1	0.5	1
	v t ∆t/∆x
	0.05	0.0166	0.016	0.0150	0.012
		(2.64)	(2.64)	(2.67)	(2.67)
	0.0625	0.0192	0.0187	0.0161	0.011
		(2.51)	(2.51)	(2.54)	(2.58)
	0.1	0.0204 (2.18)	0.0185 (2.18)	0.01 (2.27)	-1.8 × 10 -3 (2.33)
	0.25	8 × 10 -4 -7.4 × 10 -3 (1.05) (1.11)	-0.04 (1.28)	-0.08 (1.46)
	0.5	0	-0.01	-0.0508	-0.105
		(0.39)	(0.54)	(0.63)	(0.76)
	1	0	-0.0102	-0.0532	-0.112
		(0.14)	(0.27)	(0.33)	(0.39)

Table 2 :

 2 Imaginary frequency ℑω/ω p (wave number k∆x) of the dominant mode as a function of the ratio v t ∆t/∆x and the damping parameter θ f for ω p ∆t = 2 and a linear weight factor.

Table 3 :

 3 Imaginary frequency ℑω/ω p (wave number k∆x) of the dominant mode as a function of the ratio v t ∆t/∆x and the damping parameter θ f for ω p ∆t = 2 and a quadratic (n = 2) weight factor.

	ω p ∆t	1.26	2	2.83	3.46	4	5.66	6.32	8.94
	∆x/λ D	14.3	22.6	32	39.1	45.2	64	71.5	101
	linear	-0.0036 0.0034 0.0048 0.0047 0.0044 0.0036 0.0033 0.0024
		(2.09) (2.41) (2.59) (2.67)	(2.74)	(2.85)	(2.87)	(2.96)
	quadratic -0.021 -0.015 -0.01 -0.0078 -0.0066 -0.0044 -0.0039 -0.0026
		(1.95)	(2.3)	(2.5)	(2.62)	(2.68)	(2.82)	(2.85)	(2.92)
	cubic	-0.022 -0.019 -0.015 -0.013 -0.011 -0.0079 -0.0071 -0.0051
		(1.83) (2.16) (2.36) (2.48)	(2.56)	(2.7)	(2.76)	(2.85)

Table 4 :

 4 Imaginary frequency ℑω/ω p (wave number k∆x) of the dominant mode as a function of the space and time steps and the weight factor order, for a fixed ratio v t ∆t/∆x = 0.09 and θ f = 1.

Table 5 :

 5 Total energy variation and ion peak kinetic energy (keV) at 2600ω -1 p with N p = 6×10 4 .

		∆E/E 0 Ion peak energy (keV)
	Explicit	+9.3 %	232
	Implicit (θ f = 1)	-2.8%	162
	Implicit (θ f = 0.5)	+3.1%	208
	Implicit (θ f = 0.15)	+9%	273
	Implicit (θ f = 0)	+19.7%	451
		∆E/E 0 Ion peak energy (keV)
	Explicit	+1 %	221
	Implicit (θ f = 1)	-1.4%	162
	Implicit (θ f = 0.5)	+1.5%	198
	Implicit (θ f = 0.15) +4.5%	256
	Implicit (θ f = 0)	+12.4%	418

Table 6 :

 6

Total energy variation and ion peak kinetic energy (keV) at 2600ω -1 p with N p = 6×10 5 .

Table 7 :

 7 Relative variation of the total kinetic energy per time step (∆K/K 0 )/N : electrostatic case and linear weight factor.

	∆x/λ D	5	10	20	30	60
	ω p ∆t
	1	2.8 × 10 -5

Table 8 :

 8 Relative variation of the total kinetic energy per time step (∆K/K 0 )/N : electromagnetic case and linear weight factor. × 10 -4 4.6 × 10 -4 1.1 × 10 -3 2 -1.1 × 10 -4 -3.5 × 10 -5 1.4 × 10 -4 3.2 × 10 -4 8.3 × 10 -4

	∆x/λ D	5	10	20	30	60
	ω p ∆t
	1 2.3 5 -3 × 10 -5 4 × 10 -5 -1.3 × 10 -4 -2.2 × 10 -4 -10 -4	0	2.4 × 10 -4

Table 9 :

 9 Relative variation of the total kinetic energy per time step (∆K/K 0 )/N : electromagnetic case and quadratic weight factor.

  =∂ x (χ 11 -χ 11,0 )∂ x ψ + (χ 12 -χ 12,0 )∂ y ψ +∂ y (χ 21 -χ 21,0 )∂ x ψ + (χ 22 -χ 22,0 )∂ y ψ + ρ -∂ x (1 + χ 11 )E x -∂ x χ 12 E y -∂ x χ 13 E z -∂ y χ 21 E x -∂ y (1 + χ 22 )E y -∂ y χ 23 E z(98)A centered spatial discretization of Eq. (98) is given by

	1 ∆y where we have defined the source term 1 + χ 22,0 i,j+1/2 1 ∆y (ψ i,j+1 -ψ i,j ) -1 + χ 22,0 i,j-1/2 = S i,j , 1 ∆x (χ 11 i+1/2,j -χ 11,0 i+1/2 ) 1 ∆x (ψ i+1,j -ψ i,j ) -(χ 11 i-1/2,j -χ 11,0 1 ∆y (ψ i,j -ψ i,j-1 ) (97) i-1/2 ) 1 ∆x (ψ i,j -ψ i-1,j ) + 1 2∆x (χ 12 i+1,j -χ 12,0 i+1 ) 1 2∆y (ψ i+1,j+1 -ψ i+1,j-1 ) -(χ 12 i-1,j -χ 12,0 i-1 ) 1 2∆y (ψ i-1,j+1 -ψ i-1,j-1 ) + 1 2∆y (χ 21 i,j+1 -χ 21,0 i ) 1 2∆x (ψ i+1,j+1 -ψ i-1,j+1 ) -(χ 21 i,j-1 -χ 21,0 i ) 1 2∆x (ψ i+1,j-1 -ψ i-1,j-1 ) S S i,j = + + 1 ∆y (χ 22 i,j+1/2 -χ 22,0 i ) 1 ∆y (ψ i,j+1 -ψ i,j ) -(χ 22 i,j-1/2 -χ 22,0 i ) 1 ∆y (ψ i,j -ψ i,j-1 )
	--	1 ∆x 1 2∆x	(1 + χ 11 i+1/2,j )E x,i+1/2,j -(1 + χ 11 i-1/2,j )E x,i-1/2,j χ 12 i+1,j 2 E y,i+1,j+1/2 + E y,i+1,j-1/2 -χ 12 i-1,j 2 E y,i-1,j+1/2 + E y,i-1,j-1/2
	--	1 2∆x 1 2∆y	χ 13 i+1,j E z,i+1,j -χ 13 i-1,j E z,i-1,j χ 21 i,j+1 2 E x,i+1/2,j+1 + E x,i-1/2,j+1 -	χ 21 i,j-1 2	E x,i+1/2,j-1 + E x,i-1/2,j-1
	-	1 ∆y	1 + χ 22 i,j+1/2 E y,i,j+1/2 -1 + χ 22 i,j-1/2 E y,i,j-1/2
	-	1 2∆y	χ 23 i,j+1 E z,i,j+1 -χ 23 i,j-1 E z,i,j-1

  can be formulated asx n+1 -2x n + x n-1 =

					∆t 2 2	a n+1 +	a n 2	+	a n-1 2 2 +	a n-2 2 3 + . . .
	=	∆t 2 2	a n+1 +	θ f 2	a n + 1 -	θ f 2	2	a n-1 +	θ f 2	a n-2 +	θ f 2

2

a n-3 + . . .

where k g = 2π/∆x, ω g = 2π/∆t, ω q = ω -qω g and k p = k -pk g .

Exploiting the Maxwellian form of f 0 , we have

where ξ q = ωq √ 2kpvt and Z denotes the Fried and Conte plasma dispersion function Z [START_REF] Fried | The Plasma Dispersion Function[END_REF], defined by

Finally, substituting Eq. (121) into Eq. 120 yields Eq. ( 82).