Filipe Da Silva

Martin Campos Pinto

Bruno Després
email: despres@ann.jussieu.fr

Stephane Heuraux
email: stephane.heuraux@univ-lorraine.fr

Martin Campos Pinto

Stéphane Heuraux

Martin Campos
email: campos@ann.jussieu.fr

Stable coupling of the Yee scheme with a linear current model

Keywords: Maxwell's equation, fusion plasma, stable scheme, FDTD, linear current model

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The understanding of macroscopic phenomenons resulting from turbulent processes in magnetic plasmas requires huge numerical simulations. Such simulations should capture the statistical properties so that one can realize relevant averaging in time and space. This is particularly true to explain the macroscopic wave behavior which takes place on time scales that diagnostics cannot measure. An example which motivates our work is the recovery of the smooth envelope of a probing beam: to be able to compare with the theory, we need to simulate the non stationary Maxwell equations with a linear current model, see equations ( 1)-(2) below, in a domain with very fine non homogeneities (since the underlying plasma is turbulent) and over a simulation time equivalent to a half-million wave periods. This is required to study in details the role of the wavenumber spectrum on the beam widening [START_REF] Silva | A numerical study of forward-and back-scattering signatures on Doppler reflectometry signals[END_REF][START_REF] Sysoeva | Microwave beam broadening in turbulent magnetized plasmas: comparison theory-2D full-wave computations[END_REF], and is illustrated with a snapshot of a numerical calculation of a wave in a plasma in Fig. 1. The Gaussian shape is recovered at coarse scale using some time averaging, but at the same time very fine scale details are visible. Since the energy of the wave is conserved during its propagation (we will prove this property for the model used in this work), it seems unavoidable to ask for similar energy preservation properties for the numerical methods in order to capture both the fine and coarse scales in the computations.

The previous example corresponding to beam broadening induced by wave propagation in turbulent plasma is illustrative of diagnostics and wave heating systems used in magnetized fusion plasmas. Hot topics for the ITER design concern the wavenumber resolution of the so-called Doppler reflectometry as well as the beam widening induced by the turbulence on the electron cyclotron heating to determine whether neo-classical tearing modes inducing big islands can be controlled or not [START_REF] Peysson | RF current drive and plasma fluctuations[END_REF]. The same questions arise for the lower-hybrid heating system to evaluate correctly what is the wavenumber spectrum launched into the plasma in view of the efficiency prediction of the non-inductive current driven system and the definition of the energy deposition zone. In the near future, wave polarization changes will be a subject of great interest, as it is becoming with O-X mode heating scenario for the stellarator or tokamak [START_REF] Kohn | Full-wave modeling of the O-X mode conversion in the Pegasus Toroidal Experiment[END_REF] and also for the new diagnostic concepts based on wave polarization changes induced by linear mode conversion or turbulence in inhomogeneous magnetized plasmas of space, astrophysics or fusion.

To our knowledge the state of the art of the computations done for such problems use simplified wave models taking into account only one physical mechanism (refraction in [START_REF] Peysson | RF current drive and plasma fluctuations[END_REF], diffraction effects in [START_REF] Tsironis | Electron-cyclotron wave scattering by edge density fluctuations in ITER[END_REF] and recently both in [START_REF] Joo | Three-Dimensional FDTD Simulation of Electromagnetic Wave Transformation in a Dynamic Inhomogeneous Magnetized Plasma[END_REF][START_REF] Jandieri | Spatial power spectrum of multiple scattered ordinary and extraordinary waves in magnetized plasma with electron density fluctuations[END_REF]) but not all the wave structure as a Maxwell's equation solver with a linear current model can predict. This strongly questions the numerical stability of the simulations, indeed our numerical tests show a fundamental stability issue on extraordinary mode (also called TE or X-mode) which are more demanding in terms of stability than O-mode computations. Because an ideal scheme should be fast enough to reach the requirements stressed above for the comparison between theory and full-wave simulations, we are strongly interested in FDTD (Finite Difference Time Domain) schemes, indeed they are still the cheapests in terms of CPU. Many FDTD methods were developed to improve the performances and the stability in the past few years, such as EJ collocated FDTD [START_REF] Yu | An EJ collocated 3-D FDTD model of electromagnetic wave propagation in magnetized cold plasma[END_REF], Runge-Kutta exponential time differencing formulation (RKETD) FDTD [START_REF] Liu | Runge-Kutta Exponential Time Differencing FDTD Method for Anisotropic Magnetized Plasma[END_REF], matrix exponential (ME) FDTD [START_REF] Heh | FDTD modeling for dispersive media using matrix exponential method[END_REF], and exponential time differencing (ETD) FDTD [START_REF] Huang | FDTD implementation for magnetoplasma medium using exponential time differencing[END_REF]. Unfortunately they cannot be easily used to model the general dispersive media. On the other hand, in order to overcome the Courant-Friedrich-Levy (CFL) constraint of the conventional FDTD method, the one-step leapfrog ADI-FDTD method has been developed [START_REF] Cooke | A leapfrog formulation of the 3D ADI-FDTD algorithm[END_REF]. It originates from the conventional ADI-FDTD for which we refer to [START_REF] Zheng | Toward the development of a three dimensional unconditionally stable finitedifferent time-domain method[END_REF]. A re-cent variant with higher computational efficiency is the implicit method [START_REF] Tan | Fundamental schemes for efficient unconditionally stable implicit finite-difference time-domain methods[END_REF]. We finally mention the fact that conservation of energy for a fourth-order FDTD scheme is addressed in [START_REF] Liang | The spatial fourth-order energy-conserved S-FDTD scheme for Maxwell's equations[END_REF], and that massive direct simulations of wave propagation in a different physical context can be found in the recent references [START_REF] Desroziers | Simulation of laser propagation in a plasma with a frequency wave equation[END_REF][START_REF] Del Pino | 3D finite volume simulation of acoustic waves in the earth atmosphere[END_REF][START_REF] Smithe | Finite-difference time-domain simulation of fusion plasmas at radiofrequency time scales[END_REF].

The general model problem considered in this work is the non stationary Maxwell system

     -ε 0 ∂ t E + curl H = J µ 0 ∂ t H + curl E = 0 (1)
coupled with a linear equation for the electronic current density J = eN e u e , m e ∂ t u e = e (E + u e ∧ B 0 ) .

(

) 2
The unknowns are the electromagnetic field (E, H), with the usual notation that H = B µ 0 , and the electronic current is J. Equation ( 2) is the linearized Newton law for electrons where the Lorentz force on the right hand side is obtained by linearization [START_REF] Branbilla | Kinetic Theory of Plasma Waves-Homogeneous Plasmas[END_REF][START_REF] Chen | Amplification and Absorption of Electromagnetic Waves in Overdense Plasmas[END_REF][START_REF] Dumont | Effects of non-Maxwellian species on ion cyclotron waves propagation and absorption in magnetically confined plasmas[END_REF][START_REF] Stix | The theory of plasma waves, Mac Graw Hill Advanced physics monograph series[END_REF] of the full Newton law in a strongly magnetized medium with bulk magnetic field B 0 0, and vanishing electric field and current E 0 = J 0 = 0. The coefficients are the mass m e and charge e < 0 of the electron, the permittivity ε 0 and permeability µ 0 of vacuum. The electronic density N e = N e (t, x) is a prescribed function of the space and time variables. The background magnetic field could as well be dependent of the space and time variables B 0 = B 0 (t, x). The plasma frequency and cyclotron frequency are defined by ω p = e 2 N e mε 0 and

ω c = |eB 0 | m e . (3)
At any time, the total energy of the system is naturally the sum of the electromagnetic energy and the kinetic energy of electrons

E = Ω ε 0 |E| 2 2 + |B| 2 2µ 0 + m e |J| 2 2|e|N e dx. ( 4
)
This quantity will play a key role in our analysis. The specific numerical difficulty addressed here is the stability of such couplings in the regime where the electronic density N e has extremely strong gradients in space

|∇N e | >> 1.
(

) 5
The reason for such a strong gradient is that the underlying plasma is turbulent and strongly non homogeneous. See for example the density profiles used in the numerical section. On the other hand the gradient of B 0 can be considered small, even zero in first approximation. One may also consider that N e is constant in time since this very simplified situation is representative of some reflectometry problems of fusion plasmas in X-mode polarization, see the simulations on a Cartesian grid and coupled with the Yee scheme [START_REF] Yee | Numerical Solution of Initial Boundary Value Problems Involving Maxwell's Equations in Isotropic Media[END_REF] reported in [START_REF] Bohner | Simulation of Microwave Propagation in a Fusion Plasma[END_REF][START_REF] Silva | Developments on reflectometry simulations for fusion plasmas : applications to ITER position reflectometry[END_REF] and references therein. See also [START_REF] Xu | FDTD Formulations for Scattering From 3-D Anisotropic Magnetized Plasma Objects[END_REF]. This hypothesis also simplifies the analysis since it can be checked that the total energy (4) is constant in time in this regime, up to the boundary contributions. It has been observed (see e.g. Fig. 6 below) that numerical instabilities sometimes develop in the regime ∂ t N e = 0 with long simulations using well-established coupling strategies between a Yee scheme discretizing 3

(1) and the linear current equation [START_REF] Branbilla | Kinetic Theory of Plasma Waves-Homogeneous Plasmas[END_REF]. For simplicity reasons and since it does not restrict the physics neither the numerics, the bulk of this article is devoted to the case where ∂ t N e = 0, whereas the time dependent case is addressed in Appendix A. In the context of coupling the linear current equation ( 2) with Maxwell's equations (1), these numerical instabilities which are clearly related to strong gradients of the electronic density can also be characterized with the existence of different cut-offs and resonances which show up because the electronic density N e varies with space, the source terms at the boundary are time harmonic with pulsation ω and the time of simulation is large enough to reach the harmonic regime. Standard theory of waves in plasma [START_REF] Stix | The theory of plasma waves, Mac Graw Hill Advanced physics monograph series[END_REF][START_REF] Swanson | Plasma waves[END_REF] pinpoints the role of the cyclotron resonance where ω = ω c (x) and the hybridresonance where ω 2 = ω p (x) 2 +ω c (x) 2 . Since for a given ω there may always exist a solution x to these equations, it can yield a singular solution in time harmonic domain [START_REF] Després | Hybrid resonance of Maxwell's equations in slab geometry[END_REF] and so be the germ of numerical instabilities for time domain simulations. This phenomenon will be evidenced in our numerical section for the hybrid resonance. The cut-off ω = ω p (x) is not a singularity but a mild transition for propagative zone to elliptic (non propagative) and that is why it does not seem to yield the same kind of numerical difficulties as resonance do.

To be more specific let us assume now that b = B 0 |B 0 | = (0, 0, 1) t is constant. It is not difficult to see that the initial problem ( 1)-( 2) can be decomposed in two different polarizations [START_REF] Stix | The theory of plasma waves, Mac Graw Hill Advanced physics monograph series[END_REF][START_REF] Swanson | Plasma waves[END_REF] . The first polarization is transverse magnetic so the electric field together with the electric current is always parallel to b. In this case the vector product b ∧ J vanishes and the equation can be simplified a lot. This situation that does not yield resonances but only cut-off is called O-mode in plasma physics. The other polarization is transverse electric so that b ∧ J is a priori non zero. This polarization is called X-mode in the language of plasma physics. It is fundamentally the one we are interested in since the cyclotron resonance ω = ω c (x) and the hybrid-resonance 2 can show up as the result of the resonant interaction of b ∧ J 0 with the rest of the system. The equations write

ω 2 = ω p (x) 2 + ω c (x)
                   -ε 0 ∂ t E x + ∂ y H z = J x , J x = eN e u x , ε 0 ∂ t E y -∂ x H z = J y , J y = eN e u y , µ 0 ∂ t H z + ∂ x E y -∂ y E x = 0, m e ∂ t u x = eE x + eu y B 0 z , m e ∂ t u y = eE y -eu x B 0 z . (6)
This system, as the previous one (1)-( 2), is linear and will be used for the numerical tests, even if the analysis is done for the full 3D problem.

The main contributions of this work are in two directions: first we develop a theoretical framework that provides stable couplings in the energy norm between the Yee scheme (in any dimension) and a linear current model. Second, we design two simple explicit coupling schemes that are proved to be stable in the energy norm, a property that is generally not properly satisfied by existing FDTD schemes for such problems. The stability is achieved for all electronic densities, including cases with cut-offs or resonances in the computational domain. We will also prove that the CFL condition of one of our method is exactly the same as the one of the Yee scheme for Maxwell's equations in vacuum. Numerical tests performed for [START_REF] Silva | A numerical study of forward-and back-scattering signatures on Doppler reflectometry signals[END_REF] with physically based coefficients will illustrate the properties of the new scheme.

The plan is as follows. In section 2 we recall the basics of the Yee scheme coupled with a linear current model. We shall outline the formal procedure that allows to derive explicit schemes and apply a standard stability analysis to different leap-frog time discretizations. Because in Yee schemes the different components of the vector fields are defined on staggered grids, care must be taken when defining the discrete curl and vector product operators. For this reason we introduce abstract notations for these operators and we identify the conditions under which both the stability analysis and the formal derivation of an explicit scheme can be reproduced at the discrete level. In Section 3 we describe specific operators that satisfy these conditions, a property that is not fulfilled in usual FDTD schemes, see Remark 3.5. Equipped with these operators we then propose two explicit schemes that are proved to be stable under standard CFL conditions, and we give detailed formulas specifying the abstract operators. In Section 4 we finally show numerical experiments that demonstrate the improved stability of this approach.

Yee scheme with linear current model

The standard method to compute a numerical approximation of ( 1)-( 2) on a Cartesian grid consists of coupling the Yee scheme with a linear discretization of the current, see [START_REF] Xu | FDTD Formulations for Scattering From 3-D Anisotropic Magnetized Plasma Objects[END_REF][START_REF] Bohner | Simulation of Microwave Propagation in a Fusion Plasma[END_REF]. To have better insight in the structure of the coupled problem (1)-( 2) we eliminate the velocity in terms of the current and discuss a discretization procedure. Here the data N e is constant in time. Extension to the time-dependent case is postponed to Appendix A and Section 4.2.4.

Thus, the model problem is rewritten in a form which is common in plasma physics,

             -ε 0 ∂ t E + curl H = J, µ 0 ∂ t H + curl E = 0, ∂ t J = ε 0 ω 2 p (x)E + ω c b(x) ∧ J. (7)
As above, in inhomogeneous steady-state plasmas ω p (x) = |e| 2 N e (x)

mε 0 and ω c (x) = |eB 0 (x)| m e
are the plasma and cyclotron frequencies, and b(x) = B 0 (x) |B 0 (x)| is the direction of the background magnetic field.

Let n denote the time step (t n = n∆t). In the standard approach of [START_REF] Xu | FDTD Formulations for Scattering From 3-D Anisotropic Magnetized Plasma Objects[END_REF][START_REF] Bohner | Simulation of Microwave Propagation in a Fusion Plasma[END_REF], the electrical field is discretized on the integer time steps n, n + 1 ... , while the other fields are on the half time steps n -1 2 , n + 1 2 ... The algorithmic structure of the standard method reads then

                           ε 0 E n+1 -E n ∆t = curl H n+ 1 2 -J n+ 1 2 , µ 0 H n+ 1 2 -H n-1 2 ∆t = -curl E n , J n+ 1 2 -J n-1 2 ∆t = ε 0 ω 2 p E n + ω c b ∧ J n+ 1 2 + J n-1 2 2 . (8)
The equation for J h can be obtained as the limit for small friction of the one from [START_REF] Xu | FDTD Formulations for Scattering From 3-D Anisotropic Magnetized Plasma Objects[END_REF], and is equal to the one from [START_REF] Bohner | Simulation of Microwave Propagation in a Fusion Plasma[END_REF].

Formal derivation of an explicit scheme

Explicit versions of the above scheme can be derived by solving for J n+ 1 2 the third equation, which has the structure

(I -θb∧)J n+ 1 2 = W n (9) with θ = 1 2 ∆tω c and W n = ∆tε 0 ω 2 p E n + J n-1 2 + 1 2 ∆tω c b ∧ J n-1 2 .
Formally this is done by observing that since |b| = 1, the identity (b∧) 3 = -(b∧) [START_REF] Després | Hybrid resonance of Maxwell's equations in slab geometry[END_REF] holds in the sense of operators acting on vector fields, i.e., (b∧) 3 V = -(b∧)V for all V. Using (10) one easily verifies

(I -θb∧) -1 = I + θ 1 + θ 2 (b∧) + θ 2 1 + θ 2 (b∧) 2 ,
hence a formal solution for Equation ( 9) is given by

J n+ 1 2 = I + θ 1 + θ 2 (b∧) + θ 2 1 + θ 2 (b∧) 2 W n ( 11
)
and this is the argument that is used in standard approaches to derive explicit solvers for [START_REF] Silva | Development of a 2D full-wave JE-FDTD Maxwell X-mode code for reflectometry simulation[END_REF], see also Remark 3.5 below.

In the Yee scheme it is necessary to specify how the curl and vector product operators act on discrete vector fields defined on staggered grid, and from the above discussion we already see that a discrete analogue of property [START_REF] Després | Hybrid resonance of Maxwell's equations in slab geometry[END_REF] is likely to play a key role in the derivation of explicit schemes.

Discretizing on staggered grids

We consider a space discretization of (8) on a staggered mesh with spatial step h = ∆x = ∆y = ∆z, which means that E, H and J are not all on the same grid. To continue the construction, we assume that J is discretized on the edges of a Cartesian 3D cubic mesh. The center of cubic cells are denoted as x i, j,k := (ih, jh, kh), for all triplets (i, j, k) ∈ D ⊂ Z 3 where D implicitly refers to the domain of computation.

With usual notations the discretization of the electric field is

E h = E x | i, j+ 1 2 ,k+ 1 2 , E y | i+ 1 2 , j,k+ 1 2 , E z | i+ 1 2 , j+ 1 2 ,k (i, j,k)∈D
, and similarly the discrete current density has the form

J h = J x | i, j+ 1 2 ,k+ 1 2 , J y | i+ 1 2 , j,k+ 1 2 , J z | i+ 1 2 , j+ 1 2 ,k (i, j,k)∈D .
As for the magnetic field, it is discretized on the middle of faces,

H h = H x | i+ 1 2 , j,k , H y | i, j+ 1 2 ,k , H z | i, j,k+ 1 2 (i, j,k)∈D .
It is convenient to introduce some notations for the spaces of vector fields defined on these staggered grids. Thus, the space of edgewise vector fields will be denoted [START_REF] Dumont | Effects of non-Maxwellian species on ion cyclotron waves propagation and absorption in magnetically confined plasmas[END_REF] so that E h , J h ∈ V e . The space of facewise quantities is

V e = V h = V x | i, j+ 1 2 ,k+ 1 2 , V y | i+ 1 2 , j,k+ 1 2 , V z | i+ 1 2 , j+ 1 2 ,k (i, j,k)∈D
V f = V h = V x | i+ 1 2 , j,k , V y | i, j+ 1 2 ,k , V z | i, j,k+ 1
2 (i, j,k)∈D [START_REF] Heh | FDTD modeling for dispersive media using matrix exponential method[END_REF] so that H h ∈ V f . For later purposes we also introduce a notation for the space of cellwise vector fields,

V c = V h = V x | i, j,k , V y | i, j,k , V z | i, j,k (i, j,k)∈D . ( 14
)
The space discretization of (8) is easily described with these notations.

• The discretization of curl H is realized as in the classical Yee scheme [START_REF] Yee | Numerical Solution of Initial Boundary Value Problems Involving Maxwell's Equations in Isotropic Media[END_REF] by a matrix R acting from facewise quantities into edgewise quantities

R : V f -→ V e , (15)
see also Equation (54) below.

• Again as in the classical Yee scheme, the transpose operator R t : V e -→ V f will be used to compute curl E.

• Assuming that the normalized background field x → b(x) is known everywhere, the space discretization of (8) requires a discrete operator

b∧ h : V e → V e
to represent the vector product by a unitary vector field b within the space of edgewise vector fields. The exact definition and the compatibility properties of this operator will be detailed in Section 3.1.

• Finally we will need a notation for the multiplication of an edgewise vector by a scalar continuous function x → u(x). We let

S (u) : V e → V e
denote the corresponding discrete operator, where S is a reminder of its scalar nature.

Being essentially a diagonal operator, S (u) can be defined without loss of generality as

           (S (u)V h ) x | i, j+ 1 2 ,k+ 1 2 = u e | i, j+ 1 2 ,k+ 1 2 V x | i, j+ 1 2 ,k+ 1 2 (S (u)V h ) y | i+ 1 2 , j,k+ 1 2 = u e | i+ 1 2 , j,k+ 1 2 V y | i+ 1 2 , j,k+ 1 2 (S (u)V h ) z | i+ 1 2 , j+ 1 2 ,k = u e | i+ 1 2 , j+ 1 2 ,k V z | i+ 1 2 , j+ 1 2 ,k , (16)
where u e should be seen as a discrete representation of u on the edges of the grid. Of course a natural choice would consist of simply taking

           u e | i, j+ 1 2 ,k+ 1 2 = u(x i, j+ 1 2 ,k+ 1 2 ) u e | i+ 1 2 , j,k+ 1 2 = u(x i+ 1 2 , j,k+ 1 2 ) u e | i+ 1 2 , j+ 1 2 ,k = u(x i+ 1 2 , j+ 1 2 ,k ) (17)
and in this case there would be no need for a new notation. However, it will be seen in Section 3.2 that this choice does not allow to easily derive explicit schemes that are stable, for the reason that the resulting S (u) does not commute with a proper definition of the discrete operator b∧ h . With these notations, the space discretization of (8) writes

                               ε 0 E n+1 h -E n h ∆t = RH n+ 1 2 h -J n+ 1 2 h µ 0 H n+ 1 2 h -H n-1 2 h ∆t = -R t E n h J n+ 1 2 h -J n-1 2 h ∆t = ε 0 S (ω 2 p )E n h + S (ω c )b ∧ h J n+ 1 2 h + J n-1 2 h 2 . ( 18
)
To analyze further the stability of the scheme it remains to study specific operators. Note that if the scalar fields (e.g., ω p ) are constant in space, one has S (ω 2 p ) = ω 2 p , but in general this is not the case. Let us finally emphasize the fact that for the moment ( 18) is only an abstract discrete scheme: in particular it has not even been discussed whether ( 18) is an explicit or an implicit method. As we will show later, an explicit version is available for a proper design of the discrete operators that mimics the algebra ( 9)- [START_REF] Desroziers | Simulation of laser propagation in a plasma with a frequency wave equation[END_REF], see e.g. Section 3.4 below for detailed formulas.

Stability analysis

Our goal is here to analyze the stability of the abstract scheme [START_REF] Liang | The spatial fourth-order energy-conserved S-FDTD scheme for Maxwell's equations[END_REF], and obtain some conditions that the various operators should satisfy to obtain stability in the energy norm. To this end we symmetrize the continuous system and define

E = E c , J = J ω p cε 0 and B = B = µ 0 H.
Still in the case where the plasma is stationary, that is

∂ t ω p = 0, one gets                ∂ t E -c curl B = -ω p J ∂ t B + c curl E = 0 ∂ t J = ω p E + ω c b ∧ J. (19)
It can be written under the symmetric form of a Friedrichs system

∂ t U + A∂ x U + B∂ y U + C∂ z U = MU ( 20
)
with unknown U = E, B, J t . The blockwise structure of the matrices is

A = c           0 à 0 Ãt 0 0 0 0 0           B = c           0 B 0 Bt 0 0 0 0 0           C = c           0 C 0 Ct 0 0 0 0 0           with à =           0 0 0 0 0 I 0 -I 0           B =           0 0 -I 0 0 0 I 0 0           C =           0 I 0 -I 0 0 0 0 0           ,
and finally

M =           0 0 -ω p I 0 0 0 ω p I 0 M           with M = ω c           0 -b z I b y I b z I 0 -b x I -b y I b x I 0           .
In these notations the size of the identity matrix is defined in accordance with the size of the global matrix. The total physical energy in a domain Ω is

E(t) = Ω 1 2µ 0 |U(t, x)| 2 dx = Ω ε 0 |E(t, x)| 2 2 + |B(t, x)| 2 2µ 0 + m e |J(t, x)| 2 2|e|N e (x) dx. ( 21
)
Taking the scalar product of ( 20) with U, and performing basic integration by parts in the computational domain Ω, one gets the fundamental law that

d dt E(t) = 0 (22)
which expresses the conservation of the physical energy, up to boundary conditions of course. We now extend this analysis to the scheme [START_REF] Liang | The spatial fourth-order energy-conserved S-FDTD scheme for Maxwell's equations[END_REF]. Here the discrete symmetrized variables

E n h = E n h c , J n+ 1 2 h = S 1 ω p J n+ 1 2 h cε 0 and B n+ 1 2 h = B n+ 1 2 h = µ 0 H n+ 1 2 h ( 23
) are solution of                                  E n+1 h -E n h ∆t = cR B n+ 1 2 h -S (ω p ) J n+ 1 2 h B n+ 1 2 h -B n-1 2 h ∆t = -cR t E n h J n+ 1 2 h -J n-1 2 h ∆t = S (ω p ) E n h + S (ω c )b ∧ h J n+ 1 2 h + J n-1 2 h 2 . ( 24
)
We also consider the natural hypothesis that the discrete cross product satisfies

V h , b ∧ h V h h = 0, V h ∈ V e ( 25
)
where

V h , W h h = h 3 (i, j,k)∈D (V x W x )| i, j+ 1 2 ,k+ 1 2 + (V y W y )| i+ 1 2 , j,k+ 1 2 + (V z W z )| i+ 1 2 , j+ 1 2 ,k
is the natural L 2 scalar product in V e , and similar definition with a shift of indices holds concerning the norms in V f and in V c . This of course will depend on the precise definition of the operator b∧ h , but we note that it is obviously true at the continuous level. We define the discrete pseudo-energy by

Ẽn h = E n h 2 h + B n-1 2 h 2 h + J n-1 2 h 2 h -∆t E n h , cR B n-1 2 h -S (ω p ) J n-1 2 h h . (26)
Here

V h h = √ V h , V h h is the natural L 2 norm in V e ,

and later we will denote by

|R | h = max V h h ≤1 RV h h the associated matrix norm.
Proposition 2.1. Assume the orthogonality [START_REF] Tan | Fundamental schemes for efficient unconditionally stable implicit finite-difference time-domain methods[END_REF]. One has Ẽn+1 h = Ẽn h . Proof. The proof, which is a matter of easy algebra, nevertheless helps to understand the role of the key property [START_REF] Tan | Fundamental schemes for efficient unconditionally stable implicit finite-difference time-domain methods[END_REF]. We take the scalar product of the first (resp. second, resp. third) equation of ( 24) by

E n+1 h + E n h (resp. B n+ 1 2 h + B n-1 2 h , resp. J n+ 1 2 h + J n-1 2 h
). Then the difference

∆ := 1 ∆t E n+1 h 2 h + B n+ 1 2 h 2 h + J n+ 1 2 h 2 h -E n h 2 h -B n-1 2 h 2 h -J n-1 2 h 2 h reads ∆ = cR B n+ 1 2 h -S (ω p ) J n+ 1 2 h , E n+1 h + E n h h + -cR t E n h , B n+ 1 2 h + B n-1 2 h h + S (ω p ) E n h + ω c 2 b ∧ h ( J n+ 1 2 h + J n-1 2 h ), J n+ 1 2 h + J n-1 2 h
h . Many terms on the right hand side cancel, in particular the last one due to property [START_REF] Tan | Fundamental schemes for efficient unconditionally stable implicit finite-difference time-domain methods[END_REF]. One gets then

∆ = R B n+ 1 2 h -S (ω p ) J n+ 1 2 h , E n+1 h h -cR t E n h , B n-1 2 h h + S (ω p ) E n h , J n-1 2 h h = cR B n+ 1 2 h -S (ω p ) J n+ 1 2 h , E n+1 h h -cR B n-1 2 h -S (ω p ) J n-1 2 h , E n h h
which ends the proof of the claim.

Remark 2.2. In the case J ≡ 0,

ω p = ω c ≡ 0,
|R | h ≤ 2 √ 3/h in 3D. ( 27
)
In this case the CFL condition becomes c∆t < h/ √ 3 which is less stringent. The seminal Yee contribution [START_REF] Yee | Numerical Solution of Initial Boundary Value Problems Involving Maxwell's Equations in Isotropic Media[END_REF] points out c∆t < h √ 3 which is only a necessary condition.

To continue the analysis of the general scheme we set ω p L ∞ = sup x |ω p (x)|. We observe that a Cauchy-Schwarz inequality yields

Ẽn h ≥ E n h 2 h + B n-1 2 h 2 h + J n-1 2 h 2 h -c∆t |R | h E n h h B n-1 2 h h -∆t ω p L ∞ E n h h J n-1 2 h h , that is Ẽn h ≥ P α,β E n h h , B n-1 2 h h , J n-1 2 h h , α = c∆t |R | h , β = ∆t ω p L ∞ .
If the polynomial P α,β (u, v, w) = u 2 + v 2 + w 2 -αuvβuw is positive for all non zero values of (u, v, w), then Ẽn h is (the square of) a norm. This is if and only if. Assuming for example that u 0, one can define x = v u and y = v w so that P α,β (u, v, w) = x 2 + y 2 -αxβy + 1. A further reduction is based on the decomposition (x, y) = γ a + δ b where two vectors are orthonormal:

a = 1 √ α 2 +β 2 (α, β) and b = 1 √ α 2 +β 2 (β, -α).
Plugging this decomposition in the polynomial, one gets P α,β (u, v, w) = γ 2 -γ α 2 + β 2 + 1 + δ 2 . So P is positive for all non zero values of (u, v, w) if and only the discriminant of the first term in parenthesis is negative, that is

α 2 + β 2 < 4.
Thus we obtain that if the CFL condition

∆t 2 c 2 |R | 2 h + ω p 2 L ∞ < 1 ( 28
)
is satisfied, then Ẽn h is actually (the square of) a norm, which establishes the numerical stability of the scheme [START_REF] Liang | The spatial fourth-order energy-conserved S-FDTD scheme for Maxwell's equations[END_REF].

Using the estimate [START_REF] Tsironis | Electron-cyclotron wave scattering by edge density fluctuations in ITER[END_REF] for the norm of the operator R, one gets an estimate on the CFL condition ∆t 2

12c 2 h 2 + ω p 2 L ∞ 1 2 < 1.
This formula shows that the actual CFL condition is influenced by local extrema of the plasma density N e : the higher the plasma frequency ω p L ∞ , the smaller the time step. On the other hand the cyclotron frequency does not show up in this CFL condition.

Optimization of the CFL condition

Since the pseudo-energy (26) contains a term coupling the electric current with the electric field, one gets at the end of the analysis a term proportional to the maximal norm of the plasma frequency ω p L ∞ in the CFL condition [START_REF] Xu | FDTD Formulations for Scattering From 3-D Anisotropic Magnetized Plasma Objects[END_REF]. In view of this CFL condition, it is natural to ask if a slightly different method could get rid of that term ω p L ∞ , so that the CFL condition is equal to the one in vacuum,

∆t 2 c R h < 1. ( 29
)
The answer to this question is positive as is now shown. We consider the scheme

                             E n+1 h -E n h ∆t = cR B n+ 1 2 h -S (ω p ) J n+1 h + J n h 2 B n+ 1 2 h -B n-1 2 h ∆t = -cR t E n h J n+1 h -J n h ∆t = S (ω p ) E n+1 h + E n h 2 + S (ω c )b ∧ h J n+1 h + J n h 2 (30)
and the pseudo-energy

E n h = E n h 2 h + B n-1 2 h 2 h + J n h 2 h -∆t E n h , cR B n-1 2 h h . (31)
Proposition 2.3. The pseudo-energy (31) is preserved time steps after time steps by [START_REF] Yu | An EJ collocated 3-D FDTD model of electromagnetic wave propagation in magnetized cold plasma[END_REF], hence the latter is endowed with the CFL condition in vacuum.

Proof. The proof is immediate by taking the scalar product of the first equation of this system (resp. second, resp. third) by [START_REF] Cooke | A leapfrog formulation of the 3D ADI-FDTD algorithm[END_REF]. This CFL condition is independent of all plasma parameters, which may have a practical interest if ω p L ∞ is very large.

E n+1 h + E n h (resp. B n+ 1 2 h + B n-1 2 h , resp. J n+1 h + J n h ). Remark 2.
In the original variables the scheme [START_REF] Yu | An EJ collocated 3-D FDTD model of electromagnetic wave propagation in magnetized cold plasma[END_REF] reads

                             ε 0 E n+1 h -E n h ∆t = RH n+ 1 2 h - J n+1 h + J n h 2 µ 0 H n+ 3 2 h -H n+ 1 2 h ∆t = -R t E n+1 h J n+1 h -J n h ∆t = ε 0 S (ω 2 p ) E n+1 h + E n h 2 + S (ω c )b ∧ h J n+1 h + J n h 2 . ( 32
)
Remark 2.5. The extension of (32) to the time dependent case ∂ t N e 0 is addressed in Appendix A where a stable scheme is provided, and in Section 4.2.4 where preliminary numerical experiments are discussed.

Staggered grid operators and explicit schemes

In this section we will design a vector product operator acting on edgewise vector fields b∧ h : V e → V e that allows to derive stable and explicit schemes when used in conjunction with a proper operator

S (u) : V e → V e
for the multiplication of an edgewise vector field by a scalar field u.

Stable vector products on staggered grids

In the previous sections we have identified two key properties that the discrete vector product operator b∧ h should satisfy. The first one is linked with the energy stability of the scheme: it is [START_REF] Tan | Fundamental schemes for efficient unconditionally stable implicit finite-difference time-domain methods[END_REF]. The second one is linked with the existence of a local procedure to solve Equation ( 9), and hence with the existence of an explicit version of the scheme: it is the analogue of (10) and reads

(b∧ h ) 3 V = -(b∧ h )V for all V h ∈ V e . ( 33
)
To construct discrete operators on the staggered grids it will be convenient to denote by τ α x , α ∈ {-1, +1}, the half cell translation operator in the x direction,

(τ α x v)| i ′ , j ′ ,k ′ := v| i ′ + α 2 , j ′ ,k ′ for (i ′ , j ′ , k ′ ) ∈ 1 2 Z. (34)
Similarly we denote by τ β y and τ γ z the translations in the other directions. We then observe that any triplet α, β, γ ∈ {-1, +1} defines a translation operator

τ = τ α,β,γ : V e → V c (35) that maps any V h ∈ V e into V c , (τV)| i, j,k := (τ β y τ γ z V x , τ α x τ γ z V y , τ α x τ β y V z )| i, j,k for i, j, k ∈ Z (36)
see Figure 2. Of course these notations must be technically adapted near the boundaries of the domain of computation.

Averaged (centered, second order) vector products

Since the components of V h ∈ V e are defined on the edges of the grid, it is natural to design b∧ h through an averaging procedure. Specifically, one could define

                                 (b ∧ h V h ) x | i, j+ 1 2 ,k+ 1 2 := 1 4 α,β=±1 b y τ α x τ β z V z - α,β=±1 b z τ α x τ β y V y i, j+ 1 2 ,k+ 1 2 (b ∧ h V h ) y | i+ 1 2 , j,k+ 1 2 := 1 4 α,β=±1 b z τ α y τ β x V x - α,β=±1 b x τ α y τ β z V z i+ 1 2 , j,k+ 1 2 (b ∧ h V h ) z | i+ 1 2 , j+ 1 2 ,k := 1 4 α,β=±1 b x τ α y τ β z V y - α,β=±1 b y τ α x τ β z V x i+ 1 2 , j+ 1 2 ,k (37)
and indeed such formulas are often encountered in the practical implementation of Yee schemes coupled with a linear current model, see Remark 3.5

Proposition 3.1. The operator (37) is formerly second order. It satisfies [START_REF] Tan | Fundamental schemes for efficient unconditionally stable implicit finite-difference time-domain methods[END_REF] if b is a constant vector, and if D = Z 3 , but not (33).

Proof. The formulas (37) are centered so second order. The fact that (33) is not true is a standard consequence of the averaging technique. Indeed such an averaging enlarges the stencil, so at the end the stencil of (b∧ h ) 3 is strictly greater than the one of b∧ h , which implies the non satisfaction of (33). Finally property ( 25) is elementary to check on an infinite Cartesian grid.

The fact that the averaged vector product does not satisfy relation (33) has an important consequence on the stability of the numerical scheme, see Remark 3.5 below.

Clustered (non-centered, first order) vector products

In order to satisfy (33), the main idea is to design an operator b∧ h that is block diagonal, so that iterated applications on a given V h with a single non-zero entry should not result in vector fields with increasing number of non-zero entries. This is easily obtained by using local clusters of nodes and can be formalized as follows. Although the vector field b is in general known everywhere, we first introduce a discrete representation of it on the cell centers,

b c = (b c x | i, j,k , b c y | i, j,k , b c z | i, j,k ) (i, j,k)∈D ∈ V c
that may or may not coincide with the values b(x i, j,k ), but in every case satisfies

|(b c | i, j,k )| = (b c x | i, j,k ) 2 + (b c y | i, j,k ) 2 + (b c z | i, j,k ) 2 1 2 = 1 for all (i, j, k) ∈ D.
We next define the operator b∧ h as τ -1 b c ∧ τ, i.e.,

τ(b ∧ h V h ) := b c ∧ (τV h ) for all V h ∈ V e ( 38
)
where we observe that the standard product "∧" can be used since all three components of the vector fields in V c are defined on the nodes (i, j, k) ∈ D. Specifically, expanding the above definition gives

               (b ∧ h V h ) x | i, j+ β 2 ,k+ γ 2 := b c y | i, j,k V z | i+ α 2 , j+ β 2 ,k -b c z | i, j,k V y | i+ α 2 , j,k+ γ 2 , (b ∧ h V h ) y | i+ α 2 , j,k+ γ 2 := b c z | i, j,k V x | i, j+ β 2 ,k+ γ 2 -b c x | i, j,k V z | i+ α 2 , j+ β 2 ,k , (b ∧ h V h ) z | i+ α 2 , j+ β 2 ,k := b c x | i, j,k V y | i+ α 2 , j,k+ γ 2 -b c y | i, j,k V x | i, j+ β 2 ,k+ γ 2 . ( 39
)
Clearly the resulting operator b∧ h is block diagonal (the size of each block being 3 × 3 as in the standard vector product), and hence the locality of the operator is guaranteed. The degrees of freedom which are coupled by this strategy are depicted with black circles on Figure 2. Proof. The proof of these various properties is an immediate consequence of the definition of the operator. First one has that Here the vector field τV h is defined by (36) on each cell center x i, j,k (in white) from the components of V h located on the neighbouring edges (in black), namely

b ∧ h V h , V h h = τ(b ∧ h V h ), τV h h = b c ∧ τV h , τV h h = 0
V x | i, j+ 1 2 ,k+ 1 2 , V y | i-1 2 , j,k+ 1 2 and V z | i-1 2 , j+ 1
2 ,k . From (38) we verify that the discrete operator b∧ h is block-diagonal, with 3 × 3 blocks corresponding to these clusters of components.

τ((b∧ h ) 3 V h ) = (b c ∧) 3 τV h = -b c ∧ τV h = -τ(b ∧ h V h )
which ends the formal proof of (33). The fact that the operator is first order is immediate.

In particular, we observe that the formally less accurate definition (38) has stronger algebraic properties than the second order approximation (37).

Commutative multiplication by a scalar field

If x → u(x) is a continuous scalar function, the most natural way to multiply some V h ∈ V e by u consists of applying the operator S (u) defined by ( 16) together with the choice [START_REF] Kohn | Full-wave modeling of the O-X mode conversion in the Pegasus Toroidal Experiment[END_REF]. In general, however, the resulting S (u) does not commute with b∧ h , which makes the resolution procedure ( 9)- [START_REF] Desroziers | Simulation of laser propagation in a plasma with a frequency wave equation[END_REF] ineffective. Therefore, it seems natural to use a commutativity preserving operator, i.e., satisfying

S (u)(b ∧ h V h ) = b ∧ h (S (u)V h ) for all V h ∈ V e . ( 40
)
To this end we easily verify the following result. where (α, β, γ) is the triplet defining the translation operator τ in (35) and the associated vector product b∧ h in (38).

Explicit resolution

In this section we summarize the main results of this paper in a concise form using the staggered grid operators defined above. In particular we provide explicit versions of the resulting stable schemes, and refer to Section 3.4 below for detailed formulas involving the clustered (noncentered) operators introduced in Sections 3.1 and 3.2. Proposition 3.4. If the operators b∧ h and S (u) satisfy (25), ( 33) and (40), then the schemes [START_REF] Liang | The spatial fourth-order energy-conserved S-FDTD scheme for Maxwell's equations[END_REF] and (32) are stable under the respective CFL conditions [START_REF] Xu | FDTD Formulations for Scattering From 3-D Anisotropic Magnetized Plasma Objects[END_REF] and [START_REF] Yee | Numerical Solution of Initial Boundary Value Problems Involving Maxwell's Equations in Isotropic Media[END_REF]. Moreover, the scheme [START_REF] Liang | The spatial fourth-order energy-conserved S-FDTD scheme for Maxwell's equations[END_REF] has an explicit form given by

                   H n+ 1 2 h = H n-1 2 h -∆t µ 0 R t E n h J n+ 1 2 h = I + S θ 1+θ 2 (b∧ h ) + S θ 2 1+θ 2 (b∧ h ) 2 W n h E n+1 h = E n h + ∆t ε 0 (RH n+ 1 2 h -J n+ 1 2 h ) (42) with θ = ∆tω c 2 and W n h = ∆tε 0 S (ω 2 p )E n h + (I + ∆t 2 S (ω c )b∧ h )J n-1 2
h . Finally the scheme (32) has an explicit form given by

                 H n+ 1 2 h = H n-1 2 h -∆t µ 0 R t E n h J n+1 h = I + S θ 1+θ 2 (b∧ h ) + S θ 2 1+θ 2 (b∧ h ) 2 W n h E n+1 h = E n h + ∆t ε 0 RH n+ 1 2 h -1 2 (J n+1 h + J n h ) (43)
now with θ = ∆tω c 2η where η = 1 + ∆tω p 2 2 , and

W n h = S 1 η ∆tε 0 S (ω 2 p )(E n h + ∆t 2ε 0 RH n+ 1 2 h ) + (I -∆t 2 4 S (ω 2 p ) + ∆t 2 S (ω c )b∧ h )J n h .
Proof. The stability properties have been established in Sections 2.3 and 2.4 under the assumption [START_REF] Tan | Fundamental schemes for efficient unconditionally stable implicit finite-difference time-domain methods[END_REF]. To prove the explicit formulas we observe that (33) and (40) allow to reproduce the formal resolution procedure ( 9)-( 11) at the discrete level. Indeed, using these relations it is easily verified that

I + S θ 1+θ 2 (b∧ h ) + S θ 2 1+θ 2 (b∧ h ) 2 (I -S (θ)b∧ h )V h = V h ( 44
)
holds for all V h ∈ V e . Thus, rewriting the last equation from ( 18) as

(I -S (θ)b∧ h )J n+ 1 2 h = W n h ( 45
)
leads to the announced form (42). As for the new scheme (32), we observe that its implicit part reads

       E h n+1 = -∆t 2ε 0 J h n+1 + U n h , J h n+1 = ∆tε 0 2 S (ω 2 p )E h n+1 + ∆t 2 S (ω c )b ∧ h J h n+1 + V n h , (46)
with source terms

         U n h = E h n + ∆t ε 0 RH n+ 1 2 h -∆t 2ε 0 J n h V n h = J n h + ∆tε 0 2 S (ω 2 p )E n h + ∆t 2 S (ω c )b ∧ h J n h . (47)
Eliminating E n+1 h in the second equation of (46) and applying S ( 1 η ), we find

I -∆t 2 S ( ω c η )b ∧ h J n+1 h = S 1 η ∆tε 0 2 S (ω 2 p )U n h + V n h = W n h ( 48
)
so that the desired form ( 43) is again obtained by applying the identity (44).

Remark 3.5. The fact that the averaged vector product (37) does not satisfy relation (33) -and in fact, is not block-diagonal -has two practical consequences: First, it means that if one replaces the formal vector products by (37) in the semi-discrete scheme [START_REF] Silva | Development of a 2D full-wave JE-FDTD Maxwell X-mode code for reflectometry simulation[END_REF], then one obtains a scheme that is indeed stable but such that the discrete operator involved in (9) becomes non block-diagonal.

In particular, there is a priori no chance of deriving an explicit version of it. Second, it means that if one performs first the formal derivation (9)-( 11) and only then discretizes (11) with (37) as is done in [START_REF] Xu | FDTD Formulations for Scattering From 3-D Anisotropic Magnetized Plasma Objects[END_REF], then one obtains a scheme that is indeed explicit but that is not equivalent with the original formulation [START_REF] Liang | The spatial fourth-order energy-conserved S-FDTD scheme for Maxwell's equations[END_REF]. In particular, the stability analysis developped in Section 2.3 does not apply to the resulting scheme. We finally note that the difference between the averaged vector product and the clustered one is proportional to the local gradient of the density, so that the difference between the two schemes becomes significant when turbulences are present.

Detailed explicit formulas

We conclude this section by writing down the detailed formulas that result from the above procedure. To this end we assume that we have chosen an arbitrary triplet (α, β, γ) ∈ {-1, +1} 3 , and to fix the ideas we consider staggered grid operators b∧ h and S (u) defined respectively by (38) with b c | i, j,k := b| i, j,k = b(x i, j,k )

and by [START_REF] Jandieri | Spatial power spectrum of multiple scattered ordinary and extraordinary waves in magnetized plasma with electron density fluctuations[END_REF] with

u e | i, j+ β 2 ,k+ γ 2 = u e | i+ α 2 , j,k+ γ 2 = u e | i+ α 2 , j+ β 2 ,k := u| i, j,k = u(x i, j,k ).
For conciseness we now denote ĩ = i + α 2 ,  = j + β 2 and k = k + γ 2 . With the above choices the staggered grid operators correspond to

             (b ∧ h V h ) x | i, , k := b y | i, j,k V z | ĩ, j,k -b z | i, j,k V y | ĩ, j, k (b ∧ h V h ) y | ĩ, j, k := b z | i, j,k V x | i, , k -b x | i, j,k V z | ĩ, j,k (b ∧ h V h ) z | ĩ, j,k := b x | i, j,k V y | ĩ, j, k -b y | i, j,k V x | i, , k. ( 49
)
and

         (S (u)V h ) x | i, , k = u| i, j,k V x | i, , k (S (u)V h ) y | ĩ, j, k = u| i, j,k V y | ĩ, j, k (S (u)V h ) z | ĩ, j,k = u| i, j,k V z | ĩ, j,k (50)
and from Propositions 3.2 and 3.3 we know that the hypotheses of Proposition 3.4 are then satisfied. However we note that other choices could be made for b c and u e . Thus, in the explicit scheme (43) one first computes an auxiliary vector field W n h given by

                                             W n x | i, , k = ∆tε 0 ω 2 p | i, j,k E n x | i, , k + J n-1 2 x | i, , k + ∆t 2 ω c | i, j,k b y | i, j,k J n-1 2 z | ĩ, j,k -b z | i, j,k J n-1 2 y | ĩ, j, k W n y | ĩ, j, k = ∆tε 0 ω 2 p | i, j,k E n y | ĩ, j, k + J n-1 2 y | i, , k + ∆t 2 ω c | i, j,k b z | i, j,k J n-1 2 x | i, , k -b x | i, j,k J n-1 2 z | ĩ, j,k W n z | ĩ, j,k = ∆tε 0 ω 2 p | i, j,k E n z | ĩ, j,k + J n-1 2 z | i, , k + ∆t 2 ω c | i, j,k b x | i, j,k J n-1 2 y | ĩ, j, k -b y | i, j,k J n-1 2 x | i, , k (51)
and expanding the formula J

n+ 1 2 h = I + S θ 1+θ 2 (b∧ h ) + S θ 2 1+θ 2 (b∧ h ) 2 W n h for the current, one obtains                                                              J n+ 1 2 x | i, , k = 1 -θ 2 1+θ 2 | i, j,k (b 2 y -b 2 z )| i, j,k W n x | i, , k + θ 1+θ 2 | i, j,k (-b z + θb x b y )| i, j,k W n y | ĩ, j, k + θ 1+θ 2 | i, j,k (b y + θb x b z )| i, j,k W n z | ĩ, j,k J n+ 1 2 y | ĩ, j, k = 1 -θ 2 1+θ 2 | i, j,k (b 2 z -b 2 x )| i, j,k W n y | ĩ, j, k + θ 1+θ 2 | i, j,k (-b x + θb y b z )| i, j,k W n z | ĩ, j,k + θ 1+θ 2 | i, j,k (b z + θb y b x )| i, j,k W n x | i, , k J n+ 1 2 z | ĩ, j,k = 1 -θ 2 1+θ 2 | i, j,k (b 2 x -b 2 y )| i, j,k W n z | ĩ, j,k + θ 1+θ 2 | i, j,k (-b y + θb z b x )| i, j,k W n x | i, , k + θ 1+θ 2 | i, j,k (b x + θb z b y )| i, j,k W n y | ĩ, j, k. (52)
In the scheme (43), one computes an auxiliary vector field W n h given by a formula similar to (51), specifically

                                                                           W n x | i, , k = 1 η| i, j,k ∆tε 0 ω 2 p | i, j,k E n x | i, , k + ∆t 2ε 0 (RH n+ 1 2 h ) x | i, , k + (1 -∆t 2 4 ω 2 p | i, j,k )J n x | i, , k + ∆t 2 ω c | i, j,k b y | i, j,k J n z | ĩ, j,k -b z | i, j,k J n y | ĩ, j, k W n y | ĩ, j, k = 1 η| i, j,k ∆tε 0 ω 2 p | i, j,k E n y | ĩ, j, k + ∆t 2ε 0 (RH n+ 1 2 h ) y | ĩ, j, k + (1 -∆t 2 4 ω 2 p | i, j,k )J n y | i, , k + ∆t 2 ω c | i, j,k b z | i, j,k J n x | i, , k -b x | i, j,k J n z | ĩ, j,k W n z | ĩ, j,k = 1 η| i, j,k ∆tε 0 ω 2 p | i, j,k E n z | ĩ, j,k + ∆t 2ε 0 (RH n+ 1 2 h ) z | ĩ, j,k + (1 -∆t 2 4 ω 2 p | i, j,k )J n z | i, , k + ∆t 2 ω c | i, j,k b x | i, j,k J n y | ĩ, j, k -b y | i, j,k J n x | i, , k . (53)
The current field J n+1 h is then given by the same formula as before, i.e., (52) applies again to compute the entries J n+1

x | i, , k, J n+1 y | ĩ, j, k and J n+1 z | ĩ, j,k . For the sake of completeness we also recall the classical formulas corresponding to the matrix (15) that approximates the curl operator between the spaces V f and V e (the adjoint curl operator being approximated by the transpose matrix).

                                                       (RH n+ 1 2 h ) x | i, j+ 1 2 ,k+ 1 2 = 1 h (H n+ 1 2 z | i, j+1,k+ 1 2 -H n+ 1 2 z | i, j,k+ 1 2 ) - 1 h (H n+ 1 2 y | i, j+ 1 2 ,k+1 -H n+ 1 2 y | i, j+ 1 2 ,k ) (RH n+ 1 2 h ) y | i+ 1 2 , j,k+ 1 2 = 1 h (H n+ 1 2 x | i+ 1 2 , j,k+1 -H n+ 1 2 x | i+ 1 2 , j,k ) - 1 h (H n+ 1 2 z | i+1, j,k+ 1 2 -H n+ 1 2 z | i, j,k+ 1 2 ) (RH n+ 1 2 h ) z | i+ 1 2 , j+ 1 2 ,k = 1 h (H n+ 1 2 y | i+1, j+ 1 2 ,k -H n+ 1 2 y | i, j+ 1 2 ,k ) - 1 h (H n+ 1 2 x | i+ 1 2 , j+1,k -H n+ 1 2 x | i+ 1 2 , j,k ).
(54)

Numerical results

Our numerical tests are relevant for reflectometry studies in fusion plasmas or ionospheric probing. The domain of computation is 2D and we consider X-mode configuration, meaning that polarization is transverse electric as in [START_REF] Silva | A numerical study of forward-and back-scattering signatures on Doppler reflectometry signals[END_REF]. The 2D implementation maintains all the difficulties studied in this work and is simpler to implement and validate. We present three test problems on time-independent plasmas which are the main subject of this work, together with one additional time-dependent plasma case, a subject treated in Appendix A. The tests have been conducted within the REFMULX code, a 2D FTDT X-mode code developed at the IPFN (Instituto de Plasmas e Fusão Nuclear) [START_REF] Silva | Development of a 2D full-wave JE-FDTD Maxwell X-mode code for reflectometry simulation[END_REF]. Simulations were done on a grid of 1500 ×1000 points with launched frequency f = 35 GHz (that is in terms of vacuum wavelengths 65.6λ 0 × 43.8λ 0 , given a spatial resolution of h = ∆x = ∆y = λ 0 /20, where λ 0 is the the vacuum wavelength at 40 GHz). Figure 3 shows that the coefficients of the problem (in terms of f c = ω c 2π and f p = ω p 2π ) are of the same order.

Implementation details

The former implementation of REFMULX, version 1.0 (v1.0) [START_REF] Silva | Development of a 2D full-wave JE-FDTD Maxwell X-mode code for reflectometry simulation[END_REF], used a kernel referred to as the Xu-Yuan Kernel (XYK) whose original description appears in [START_REF] Xu | FDTD Formulations for Scattering From 3-D Anisotropic Magnetized Plasma Objects[END_REF]. We use it in some of the tests to show its limitations in the context of this work.

The new version of REFMULX, version 2.0 (v2.0), includes two kernels: (i) the explicit scheme (42) equivalent to [START_REF] Liang | The spatial fourth-order energy-conserved S-FDTD scheme for Maxwell's equations[END_REF], which corresponds to a rewriting of the XYK scheme where the symmetric vector products (37) have been substituted by the non-centered formulas (38)-(39) (with a given arbitrary translation pattern fixed once and for all) ; and (ii) the new conservative scheme (43) equivalent to (32), proposed in this work to improve the CFL constraint.

The CFL condition has been checked to satisfy the numerical formula [START_REF] Xu | FDTD Formulations for Scattering From 3-D Anisotropic Magnetized Plasma Objects[END_REF] which is enough to get the numerical stability in the energy norm, but can be also considered as a physical evaluation of the speed of numerical waves.

Tests

The tests use physically based coefficients.

f U = (4π) -1 (ω p (x) 2 + ω 2 c ) 1/2
+ ω c (above the cyclotron frequency f c = ωc 2π ≈ 26.6 GHz) are presented for a magnetic field of B 0 = 0.95T. In zones where n e = 0 (meaning no plasma), we arbitrarily set f U = 0 thus following the physical intuition. This is visible in the bottom figure where f U vanishes in the same zones than n e . The insets show the sharp discontinuities at the left plasma edge, and the smooth transitions used in some of the test cases.

The plasma setup has a background density denoted as N e 0 (x) which increases linearly along x (radial direction), with its edge at x = 500 grid points (21.9λ 35 GHz ). Turbulence δN e (x, y) is added to the background plasma N e (x, y) = N e 0 (x) + δN e (x, y). In Figure 3 a typical plasma density profile is depicted in the top figure where the electronic density is plotted against the horizontal axis (measured in cell points knowing that h = ∆x = ∆y = λ 0 /20, where λ 0 is the vacuum wavelength at 40 GHz). The background density profile N e 0 (x) (without turbulence, i.e. 0%) is used in the simulations described in 4.2.1-4.2.2. Simulation tests on 4.2.3-4.2.4 are done with turbulent plasma with levels δN e (x 35 )/N e (x 35 ) ranging from 5%-50% in RMS: here x 35 is the position of the cut-off for a probing wave launched at 35 GHz. In Figure 3 density profiles at different turbulence levels 10%, 30% and 50% are shown. In the inset a zoom of the edge density (only for the 50% curve, for clarity) is shown. Its meaning will be explained in 4.2.3. The bottom of Figure 3 shows the upper cut-off frequency f U for a constant external magnetic field of B 0 = 0.95T: the upper cut-off frequency is above the electronic cyclotron frequency f c = ω c 2π ≈ 26.6 GHz. Be careful that for physical reasons, we decided arbitrarily that a vanishing electronic density n e corresponds to f U = 0.

The first test presents a comparison of ( 42) and ( 43) to verify that the pseudo-energy is conserved.

In the second set of tests very long simulations are performed to show the influence of roundoff errors on the numerical stability/instability when second order symmetric formulas like (37) are used. Here there are no density fluctuations, so the problem is without turbulence.

After these first numerical tests, the following ones are more closely connected to scenarios of interest for the fusion plasma community. Different types of turbulence have been tested to determine if the numerical stability is affected by the turbulence properties. As Gaussian wavenumber spectrum is representative of the turbulence behavior, this kind of turbulence is used in all the following tests. The different test cases appear from the simplest to the more general one.

In the third family of tests we consider steady state turbulent cases and we assess the improved stability offered by the clustered (non-centered) vector-product (38). Results are shown for different levels of turbulence, and we specify some of the parameters used to stabilize the computations over simulation times representing at least one million time steps.

The fourth tests address an extension of the proposed scheme in order to assess the stability for time-dependent turbulence-wave interaction simulations, which are closer to experiments. There the turbulence has a temporal dependency N e (x, y, t) = N e 0 (x) + δN e (x, y, t), being able to evolve with a velocity along y (poloidal) and x (radial) directions.

Test 1: pseudo-energy conservation

Here we compare the results obtained with the new explicit schemes (42) and (43), both implemented with the clustered (non-centered) vector product (38). This test intends to verify the theoretical conservation properties in the numerical results for an extremely steep initial signal. The initial data is a Dirac mass, and we represent in figure 4 the 2D version of the pseudo-energy [START_REF] Zheng | Toward the development of a three dimensional unconditionally stable finitedifferent time-domain method[END_REF] versus the time steps. The Dirac mass is located 100 grid cells inside the plasma, therefore being at i s = 600 where i s = 500 is the beginning of the electronic density ramp, see Figure 3 (curve 0%). The external magnetic field used in this set of runs was B 0 = 0.95T .

At the position (i s , j s ) = (600, 500), an impulse is excited at iteration n = 0, with a value

H z (is, js, 0) = 1 µ 0 1 √ ∆x∆y
to have a energy density of 1 ∆x∆y resulting in an initial pseudo-energy E 0 h = 1, see [START_REF] Zheng | Toward the development of a three dimensional unconditionally stable finitedifferent time-domain method[END_REF]. The number of iterations considered is n = 700, not long enough to allow a response from the boundary zone. Results appear in Figure 4, where the exact preservation of the pseudo-energy by the scheme (43) confirms the theory, see Proposition 2.3. We also note that according to Proposition 2.1, the pseudo-energy should be exactly preserved by the scheme (42).

Test 2: stability of wave propagation (without turbulence)

In this series of tests we assess the long time stability of the code in the absence of turbulence. Thus the same linear plasma density N e 0 is considered as in the previous test cases (cf. Figure 3, curve 0%), and a wave is launched with frequency f = 35 GHz from the left part of the domain, as in Figure 8 . Here the standard XYK scheme (involving averaged vector products and implemented in REFMULX v1.0) shows to be stable up to 100 thousand iterations but for longer runs is unstable. By contrast, on REFMULX v2.0 both kernels (42) and (43) show an augmented stability up to 3.5 million iterations. Increasing further the duration of the runs we have verified that they became unstable when using FLOAT. Written in DOUBLE the new kernels reached 8 million iterations in our trials. The natural interpretation is that the control of roundoff errors is important for such difficult problems. In Figure 5 we present the longest run tried in these tests, using a logarithmic scale to better visualize the exponential growth of the unstable run.

Test 3: wave propagation in steady state turbulent plasma

The aim of this part is to detail some points on the stability domain of the computation of X-mode propagation in turbulent plasma. We start with the standard scheme XYK [START_REF] Xu | FDTD Formulations for Scattering From 3-D Anisotropic Magnetized Plasma Objects[END_REF], used on v1.0 (with averaged vector products), and further ahead use both the schemes (42) and (43) implemented in v2.0.

The original XYK, which was already seen to be unstable without turbulence for simulation times of about 100 thousand iterations, shows very low stability properties when adding high levels of turbulence. In fact, even before 100 thousand iterations the runs are seldom stable. The characterization of these instabilities with turbulence has some interest to better understand the phenomenon. In Figure 6 we present the results obtained with XYK (REFMULX v1.0), using a logarithmic scale to better visualize the exponential growth of the instabilities. In the top part Here the use of FLOAT results in a run that is unstable after ≈ 3.5 million time steps, and a logarithmic scale is used to better visualize the exponential growth.

21

Using DOUBLE the run remains stable during the 8 million time steps of the complete simulation. Similar results were obtained with the new scheme (43) and are not shown here.

of Figure 6 we show the runs obtained using turbulent plasma density profiles where on the left plasma edge a sharp discontinuity is present, as seen in Figure 3 for three levels of turbulence. In the lower part of Figure 6 we show the runs obtained by replacing those sharp discontinuities by a smooth transition, as seen in the inset of Figure 3 (only for the 50% case, for clarity). Although very localized, this slight modification of N e has a considerable impact on the overall results: for instance, with a smooth plasma edge the 10% case becomes stable.

In Figure 7 we next represent the time steps where the instabilities are triggered (upper part) together with their growth rates (lower part) for several levels of turbulence. We thus verify that in the presence of a sharp discontinuity at the plasma edge, the instabilities are triggered sooner and they grow faster than in the case where the transition is smooth. And in both cases we observe that the same occurs (earlier instabilities and higher growth rates) as the turbulence amplitude increases. Note that these results are to be looked at as a behavior guideline, not as an absolute value: indeed, using a different turbulence matrix profile δN e with same spectrum we may very well obtain an unstable run in the 10% case with a smooth plasma edge.

In Figures 8 and9 we next show numerical results obtained with the new explicit scheme (43) involving non-centered vector products, as implemented in REFMULX v2.0. Here we essentially observe that all the problems that were unstable before become stable, and that the spurious exponential growths illustrated in Figure 6 are prevented. This confirms in practice the theoretical stability result proved in Proposition 3.4, in particular this shows that in this case the presence of boundary conditions does not spoil the numerical stability of the scheme. In Figure 8 we represent a typical snapshot that is representative of all these stable runs. There the field H z is plotted at iteration 80, 000 for a turbulence level of 10%, together with the distorted contours of the upper cutoff f U = 35 GHz (defined in Figure 3). In Figure 9 we next show the runs obtained with the scheme (43) with different levels of turbulence and a sharp plasma edge (see Figure 3). When compared with the runs obtained with the centered vector product, the result is outstanding in terms of stability. Note that simulations are now stable at least up to 1 million iterations. Results with the other new scheme (42) present the same level of stability, and are not shown here.

As we have identified that the spatio-temporal evolution (stability) is highly connected to the 6 and on bottom their growth rates. Here the framed cases correspond to a smooth transition at the plasma edge, see Figure 3.

behavior to sharp density gradients at the plasma edge, the use of a smoothed density fluctuations profile can play a role on stability improvement. Another important point when implementing the kernels is to optimize the numerical operations by minimizing the number of operations, especially divisions, which not only increases performance but also improves stability for such challenging problems. This improvement is enhanced when one switches from single to double precision, extending the validity domain of both the new kernels (42) and (43) to more than one million iterations. This number is enough to treat cases relevant for fusion reactors or ITER. However if one wants to study the probing beam widening induced by density fluctuations [START_REF] Sysoeva | Microwave beam broadening in turbulent magnetized plasmas: comparison theory-2D full-wave computations[END_REF], the requirements are even higher. The turbulence profile δN e needs to move so that the wave intensity can be averaged over a number of time steps, and we are closer to the limits of HPC machines. This motivates the following section, devoted to time-dependent cases.

Test 4: time-dependent plasma

Although not the main subject of this article, the case of a time-dependent plasma density is addressed here because of its importance on some problems in plasma physics. Thus, our objective is not to provide a full treatment of this issue, which we divert to a later study, but to alert to some of the problems it poses, namely as stability is concerned. In the time-dependent case indeed an additional term ∂ t N e /N e is introduced, and our numerical scheme needs to be modified as described in Appendix A. As a consequence, we observe that numerical problems can arise at the plasma edge where the plasma density N e can be very low and the variations ∂ t N e very high. Again, this can induce numerical instabilities. Moreover this problem is more difficult to analyze than before since the total energy is no more conserved, due to the source term proportional to ∂ t N e , see (A.5) below.

To reduce this potential instability effect, we will use some heuristics deduced from the previous tests. Indeed one can use the smoothed edge plasma (see Figure 3) to reduce the local density Figure 8: Snapshot of the positive part of the magnetic field, that is max (H z , 0), taken at iteration 80, 000 with a turbulence level of 10%, using the new scheme (43). The distorted countour lines corresponds to the upper cutoff f U = 35 GHz as defined in the caption of Figure 3. gradient and to regularize the plasma density profile. We introduce two directions in which the turbulence matrix δN e moves. One is parallel to the plasma edge and it is called poloidal (using tokamak terminology as a reference) and the other is co-linear to the plasma density gradient and is the radial direction. Different studies have been done as a function of the turbulence level, the direction, and velocity of the motion. The movement of the density is done sliding the profile matrix δN e along the required direction (poloidally or radially), shifting the lines or columns within the discretization grid. The velocity of a fusion plasma is much slower than the time frame of the electromagnetic wave (V ≪ c), with top values in the order of km/s. To speed up simulations we allow the plasma to move faster than the real velocities, correcting the results to accommodate for the higher simulation speed [START_REF] Silva | A numerical study of forward-and back-scattering signatures on Doppler reflectometry signals[END_REF]. Nevertheless, there are limits to that scaling and the plasma velocity must be kept at a fraction of the speed of light to maintain physical meaning [START_REF] Silva | Studies on O-mode reflectometry spectra simulations with velocity shear layer[END_REF], [START_REF] Silva | A numerical study of forward-and back-scattering signatures on Doppler reflectometry signals[END_REF], usually around 2% of c, always less than 10%. Therefore the number of iterations for a δ x,y displacement of the matrix δN e to occur is ∆n = 2(V/c) -1 . While ensuring the correct velocity, this shift of the matrix every ∆n iterations involves a sudden variation in the value of N e at each point. A better method is to check at each iteration what fraction of ∆n has already elapsed and set N e as a weight of the present density and a (1weight) of the next value, securing a much smoother motion with the same velocity and showing much cleaner spectra. The movement of the turbulence matrix δN e reduces the stability and the correction ∂ t N e /N e limits it even further. In all tests the time dependence reduces severely the level of turbulence were the code holds its stability as well as the duration of the stable runs. Introducing a poloidal movement with V/c = 1.4% on -0. Here every case is stable, although with more amplitude for a high level of turbulence. Note that the number of time steps is much greater than in Figure 6 to illustrate the improved stability of the method. Runs obtained with the scheme (42) also implemented in REFMULX v2.0 show the same level of stability and are not presented here.

the simulation set-up used in the previous tests, neglecting the correction ∂ t N e /N e , a stable run of 140 thousand iterations is obtained for a turbulence level of 10%. Including the correction ∂ t N e /N e , the maximum turbulence level drops to 5%. Preliminary tests with radial displacement of the plasma seem to be more stable than the poloidal case when neglecting ∂ t N e /N e but the degradation in stability brought by the correction is much higher than in the poloidal case. Partial conclusions to have a stable code in the time-dependent case are: reduce the turbulence level at the edge, be sure that turbulent matrix δN e has no discontinuity inside, optimize the numerical operation avoiding division, use double precision.

Conclusion

We have detailed the influence of strong gradients in the plasma parameters on the stability in the energy norm of a Yee scheme coupled with a linear current model. The main asset is the identification of an instability phenomenon triggered by strong gradients of N e . This is related to the use of a centered evaluation of the vector product. Different solutions have been proposed to get a scheme guaranteed stable in energy norm, all of them being related to the use of a clustered vector product acting on edgewise vector fields that satisfies two compatibility relations. The first one,

V h , b ∧ h V h h = 0 for V h ∈ V e ,
is in our context necessary to get a scheme stable in the energy norm. The second one reads

(b∧ h ) 3 V h = -(b∧ h )V h for V h ∈ V e
and expresses the fact that, as an operator acting on edgewise vector fields, b∧ h is block-diagonal.

As such, it allows to get an explicit scheme with the same CFL condition as in vacuum. Both conditions are very natural since there are true of course for the exact operator, given that |b| = 1.

Our numerical tests confirm the theoretical analysis and they assess the remarkable stability property of the compact first order vector product. Further research will be devoted to the study of moving plasmas for which ∂ t N e 0 and to full 3D implementations involving all wave components. These are of main interest to treat cases closer to experiments where the full description should include turbulence and plasma motions. This stage is a path to the introduction of synthetic diagnostics in turbulence codes to predict signals coming from either reflectometer signals or scattering diagnostics. Overall, we emphasize that the applications of the numerical schemes introduced in this article are wider than those presented. They include the beam widening of electron cyclotron heating as well as other heating system using waves, space physics or telecommunications.

Appendix A. Time dependent plasma densities

If the electronic density N e varies with time, then the current density satisfies Turning to the symetrized variables E = E c , J = J ω p cε 0 and B = B = µ 0 H we observe that the time-varying plasma density yields ∂ t J = ∂ t J ω p cε 0 -∂ t ω p ω p J. We thus obtain

             ∂ t E -c curl B = -ω p J ∂ t B + c curl E = 0 ∂ t J = ω p E + ω c b ∧ J + ∂ t ω p ω p J.
(A.2)

Discretizing this system as in (30) yields then . Clearly, this is not an intuitive discretization of (A.1).

                                       E n+1 h -E n h ∆t = cR B n+ 1 2 h -S (ω n+ 1 2 p ) J n+ 1 2 h B n+ 1 2 h -B n-
                                                     ε 0 E n+1 h -E n h ∆t = RH n+ 1 2 h - 1 2          S ω n+ 1 2 p ω n p J n h + S ω n+ 1 2 p ω n+1 p J n+1 h          µ 0 H n+ 3 2 h -H n+ 1 2 h ∆t = -R t E n+1 h J n+1 h -J n h ∆t = ε 0 S (

Appendix B. Resolution with time averaged operators

We here detail how the formulas in Section 3.3 are modified in the case where some averaging in time is introduced through the use of a time-dependent cluster pattern

(α n , β n , γ n ) ∈ {-1, +1} 3 , n ∈ N. (B.1)
In such a case we observe that in order to satisfy Properties ( 25), ( 33) and (40) the staggered grid operators should also change in time. Let us then denote by b∧ n h and S n (u) those defined with the triplet (α n , β n , γ n ).

Repeating the stability analysis from Section 2.4 we then verify that the discrete scheme

                             E n+1 h -E n h ∆t = cR B n+ 1 2 h -S n (ω p ) J n+1 h + J n h 2 B n+ 1 2 h -B n-1 2 h ∆t = -cR t E n h J n+1 h -J n h ∆t = S n (ω p ) E n+1 h + E n h 2 + S n (ω c )b ∧ n h J n+1 h + J n h 2 (B.2)
(which can be made explicit with the procedure outlined in Section 3.3) allows to compute

( E n h , B n-1 2 h , J n h ) → ( E n+1 h , B n+ 1 2
h , J n+1 h ) in such a way that the pseudo-energy (31)

E n h = E n h 2 h + B n-1 2 h 2 h + J n h 2 h -∆t E n h , cR B n-1 2 h h (B.3)
is preserved. In particular, we see that the use of a time-dependent cluster pattern (B.1) does not change anything on the symetrized formulation. However it affects the scheme formulated in the original variables, due to the fact that being pattern-dependent the change of variables now becomes also time-dependent, namely

J n h = cε 0 S n (ω p ) J n h . (B.4)
The first effect is that the discrete pseudo-energy, seen as a functional of the field J n h , changes with time. The second effect is a slight modification of the numerical scheme (32). Indeed applying (B.4) we find that (32) reformulated in the original variables reads

                                     ε 0 E n+1 h -E n h ∆t = RH n+ 1 2 h - Jn+1 h + J n h 2 , µ 0 H n+ 3 2 h -H n+ 1 2 h ∆t = -R t E n+1 h , Jn+1 h -J n h ∆t = ε 0 S n (ω 2 p ) E n+1 h + E n h 2 + S n (ω c )b ∧ h Jn+1 h + J n h 2 , J n+1 h = S n+1 (ω p )S n (ω -1 p ) Jn+1 h .
(Although it has not been used here, we may note that S n (u) and S n+1 (v) commute as diagonal operators.) Clearly the explicit version of this scheme is obtained by a straightforward application of Proposition 3.4: it reads

                           H n+ 1 2 h = H n-

Figure 1 :

 1 Figure 1: Contours at a given time of the non negative part of the electric field, that is max (E z , 0), of an ordinary wave (known as TM or O-mode) with a Gaussian shape in unperturbed plasma with a linear density profile (left) and with turbulence of δn/n c = 5% (RMS value, with n c the density at the cut-off layer for ν = 40 GHz) where the sub-beams are present (right).

Proposition 3 . 2 .

 32 The operator b∧ h defined in (38) is a first order approximation of the exact operator. It satisfies both[START_REF] Tan | Fundamental schemes for efficient unconditionally stable implicit finite-difference time-domain methods[END_REF] and (33) for general b and D ⊂ Z 3 .

Figure 2 :

 2 Figure2: Graphical depiction of a cluster of staggered degrees of freedom defined by a translation τ = τ α,β,γ which pattern is (α, β, γ) = (-1, +1, +1). Here the vector field τV h is defined by (36) on each cell center x i, j,k (in white) from the components of V h located on the neighbouring edges (in black), namely V x | i, j+ 1

  holds for any vector field b and domain D, since every term involved in the (i, j, k)-sum underlying b c ∧ τV h , τV h h is zero. Next from |b c | = 1 one verifies that (b c ∧) 3 = -(b c ∧). Applying several times (38) this gives

Proposition 3 . 3 .2 ,k+ γ 2 = u e | i+ α 2 , j,k+ γ 2 =

 33222 The operator S (u) defined by[START_REF] Jandieri | Spatial power spectrum of multiple scattered ordinary and extraordinary waves in magnetized plasma with electron density fluctuations[END_REF] satisfies (40) if and only if the edgewise representation u e of the scalar field u satisfies u e | i, j+ β u e | i+ α 2 , j+ β 2 ,k for (i, j, k) ∈ D, (41) 14

Figure 3 :

 3 Figure3: On top, the plasma background density N e 0 (x) (0%) used for the non turbulent tests is shown together with turbulent densities N e (x, y) = N e 0 (x) + δN e (x, y) for three levels of turbulence δN e (x 35 )/N e (x 35 ) of 10%, 30% and 50% RMS (here x 35 is the position of the cut-off for a probing wave launched at 35 GHz). On bottom, the corresponding upper cutoff frequencies f U = (4π) -1 (ω p (x)2 + ω 2 c ) 1/2 + ω c (above the cyclotron frequency f c = ωc 2π ≈ 26.6 GHz) are presented for a magnetic field of B 0 = 0.95T. In zones where n e = 0 (meaning no plasma), we arbitrarily set f U = 0 thus following the physical intuition. This is visible in the bottom figure where f U vanishes in the same zones than n e . The insets show the sharp discontinuities at the left plasma edge, and the smooth transitions used in some of the test cases.

Figure 4 :

 4 Figure 4: Pseudo-energy (31) of the numerical solutions computed with the schemes (42) (in green) and (43) (in red), starting from an initial Dirac impulse far from the boundaries. The exact preservation of the pseudo-energy by (43) is in accordance with the theory, see Proposition 2.3.

Figure 5 :

 5 Figure5: Amplitude of the solution H z computed with the new scheme (42) implemented in REFMULX v2.0 versus iteration number for the non-turbulent test case described in Section 4.2.2. Here the use of FLOAT results in a run that is unstable after ≈ 3.5 million time steps, and a logarithmic scale is used to better visualize the exponential growth. Using DOUBLE the run remains stable during the 8 million time steps of the complete simulation. Similar results were obtained with the new scheme (43) and are not shown here.

Figure 6 :

 6 Figure 6: Return signals 20 log 10 H z L ∞ for wave simulations in a steady state turbulent plasma using the XYK scheme implemented in REFMULX v1.0 with a sharp edge (top) and a smooth transition (bottom) in the plasma density profile, see Figure 3 and details in the text.

2 ]

 2 δn e (x 35GHz )/n e (x 35GHz )[% RMS] w/o frame w/ frame

Figure 7 :

 7 Figure7: On top the trigger instants of the instabilities presented in Figure6and on bottom their growth rates. Here the framed cases correspond to a smooth transition at the plasma edge, see Figure3.

Figure 9 :

 9 Figure9: Return signals are shown in linear scale for different levels of turbulence using the new scheme (43) implemented in REFMULX v2.0 which involves the clustered (non-centered) vector product (38). Here every case is stable, although with more amplitude for a high level of turbulence. Note that the number of time steps is much greater than in Figure6to illustrate the improved stability of the method. Runs obtained with the scheme (42) also implemented in REFMULX v2.0 show the same level of stability and are not presented here.

∂= 2 ∂

 2 t J = eN e ∂ t u e + J ∂ t N e N eand using ∂ t N e N e t ω p ω p we see that[START_REF] Silva | Developments on reflectometry simulations for fusion plasmas : applications to ITER position reflectometry[END_REF] becomes                 -ε 0 ∂ t E + curl H = J, µ 0 ∂ t H + curl E = 0, ∂ t J = ε 0 ω 2 p E + ω c b ∧ J + 2J ∂ t ω p ω p .(A.1)

  Note that it is possible to rewrite the scheme (A.3) in the original variables. Because of the time dependence of ω p , the latter are given by E n h =

	where E	n+ 1 2 h	= 1 2 ( E n h + E n+1 h ) and J n+ 1 2 h	= 1 2 ( J n h + J n+1 h ). For this scheme the pseudo-energy
				E n h = E n h	2 h + B n-1 2 h	2 h + J n h	2 h -∆t E n h , cR B n-1 2 h	h	(A.4)
	satisfies			E n+1 n ∆t -E n n	= S	ω n+1 p -ω n p ∆tω n+ 1 p 2	J n+ 1 2 h , J	n+ 1 2 h
	which is a centered approximation of the variation of the exact energy (21),
					d dt	E(t) =	Ω	∂ t ω p (t, x) 2µ 0 ω p (t, x)	| J(t, x)| 2 dx.	(A.5)
	allows to rewrite (A.3) as			E n h c , J n h = S (ω n p )	J n h cε 0 and B n+ 1 2 h	= µ 0 H	n+ 1 2 h	which
				1
				2
			∆t	h	= -cR t E n h	(A.3)
			J n+1 h ∆t -J n h	= S (ω n+ 1 2 p ) E n+ 1 2 h	+ S (ω c )(b∧ h ) + S	ω n+1 p -ω n p ∆tω n+ 1 p 2	J	n+ 1 2 h

  ω

								(A.6)
					n+ 1 2 p ω n+1 p )E n+ 1 2 h	+ S	ω n+1 p -ω n p ∆tω n p	J n h
						+ S (ω c )(b∧ h ) + S	ω n+1 p -ω n p ∆tω n+ 1 p 2	2 Jn+ 1 h
	where E	n+ 1 2 h	= 1 2 (E n h + E n+1 h ) and	Jn+ 1 2 h	= 1 2 S	ω n+1 p ω n p	J n h + J n+1 h

Acknowledgments

The authors acknowledge the support of ANR under contract ANR-12-BS01-0006-01. Moreover this work was carried out within the framework of the European Fusion Development Agreement and the French Research Federation for Fusion Studies. It is supported by the European 29 Communities under the contract of Association between Euratom and CEA, and has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement number 633053. IST activities also received financial support from "Fundac ¸ão para a Ciência e Tecnologia" through project Pest-OE/SADG/LA0010/2013. The views and opinions expressed herein do not necessarily reflect those of the European Commission.