Éliane Bécache
email: eliane.becache@inria.fr

Patrick Joly
email: patrick.joly@inria.fr

Maryna Kachanovska
email: maryna.kachanovska@ensta.fr

Stable Perfectly Matched Layers for a Cold Plasma in a Strong Background Magnetic Field $

Keywords: plasma, perfectly matched layer, Maxwell equations, stability, staggered grids, Fourier analysis 2000 MSC: 65M12

This work addresses the question of the construction of stable perfectly matched layers (PML) for a cold plasma in the infinitely large background magnetic field. We demonstrate that the traditional, Bérenger's perfectly matched layers are unstable when applied to this model, due to the presence of the backward propagating waves. To overcome this instability, we use a combination of two techniques presented in the article.

First of all, we consider a simplified 2D model, which incorporates some of the difficulties of the 3D case, namely, the presence of the backward propagating waves. Based on the fact that for a fixed frequency either forward or backward propagating waves are present, we stabilize the PML with the help of a frequency-dependent correction.

An extra difficulty of the 3D model compared to the 2D case is the presence of both forward and backward waves for a fixed frequency. To overcome this problem we construct a system of equations that consists of two independent systems, which are equivalent to the original model. The first of the systems behaves like the 2D plasma model, and hence the frequency-dependent correction is added to the PML for the stabilization. The second system resembles the Maxwell equations in vacuum, and hence a standard Bérenger's PML is stable for it. The systems are solved inside the perfectly matched layer, and coupled to the original Maxwell equations, which are solved in a physical domain, on a discrete level through an artificial layer.

The numerical experiments confirm the stability of the new technique.

Introduction

In this work we concentrate on the numerical modelling of wave propagation in cold strongly magnetized plasmas, in a particular case when the size of the physical domain (e.g. a tokamak) drastically exceeds the size of the domain of interest (e.g. a smaller enclosure inside a tokamak modelled during, loosely speaking, relatively short times). In this case it is reasonable to study the corresponding system posed in a free space or in a semi-infinite domain (e.g. bounded by the edge of plasma in a tokamak), see [START_REF] Jacquot | 2d and 3d modeling of wave propagation in cold magnetized plasma near the tore supra icrh antenna relying on the perfecly matched layer technique[END_REF][START_REF] Jacquot | Description non lineáeire auto-cohérente de la propagation d'ondes radiofréquences et de la périphérie d'un plasma magnétisé[END_REF] and references therein, which, in turn, requires special treatment of the unboundedness. For the application of such modelling in plasma physics we refer a reader to the work [START_REF] Heuraux | Study of wave propagation in various kinds of plasmas using adapted simulation methods, with illustrations on possible future applications[END_REF].

To deal with the unboundedness of the computational domain we wish to employ the perfectly matched layer (PML) technique suggested by Bérenger in his seminal articles [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Bérenger | Three-dimensional perfectly matched layer for the absorption of electromagnetic waves[END_REF], which consists of surrounding the domain with an artificial layer. Inside this layer the original equations are modified so that the solution decays rapidly, and the truncation of the layer with zero boundary conditions would result in a negligible reflection. Crucially, the perfectly matched property ensures that there is no reflection on the interface between the physical domain and the perfectly matched layer. However, this method is known [START_REF] Tam | Perfectly matched layer as an absorbing boundary condition for the linearized Euler equations in open and ducted domains[END_REF][START_REF] Hu | A stable, perfectly matched layer for linearized Euler equations in unsplit physical variables[END_REF][START_REF] Cummer | Perfectly matched layer behavior in negative refractive index materials[END_REF] to suffer of instabilities, connected to the presence of so-called backward propagating modes [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF], which are the consequence of the anisotropic [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF] or dispersive [START_REF] Cummer | Perfectly matched layer behavior in negative refractive index materials[END_REF][START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF] nature of the system. Although much effort had been dedicated to the construction and analysis of stable perfectly matched layers for general hyperbolic systems, see e.g. [START_REF] Appelö | A general perfectly matched layer model for hyperbolic-parabolic systems[END_REF][START_REF] Appelö | Perfectly matched layers for hyperbolic systems: general formulation, well-posedness, and stability[END_REF], so far, to our knowledge, there exists no general recipe to deal with such instabilities, and each of the problems requires a separate treatment, see e.g. [START_REF] Nataf | New constructions of perfectly matched layers for the linearized Euler equations[END_REF][START_REF] Abarbanel | Well-posed perfectly matched layers for advective acoustics[END_REF][START_REF] Diaz | Stabilized perfectly matched layer for advective acoustics[END_REF][START_REF] Demaldent | Perfectly matched transmission problem with absorbing layers: application to anisotropic acoustics in convex polygonal domains[END_REF][START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF].

In this work we consider a simplified model of cold plasma in a strong background magnetic field. The behaviour of cold plasmas is described in the frequency domain by Maxwell equations with the anisotropic frequency-dependent hermitian tensor of dielectric permittivity. For strong background magnetic fields, a simplification of this model results in a diagonal, but anisotropic and dispersive tensor of dielectric permittivity. The study of this particular case was suggested to us by L. Colas (CEA Cadarache). Already this simplified model presents difficulties for the PMLs, and in this work we present a method to overcome them. Our idea is based on exploiting a very particular structure of the dispersion relation of the Maxwell equations with the plasma term, which allows us to split the original model into two different systems, to which we can apply two different kinds of the PMLs. This idea is new in the context of plasmas, however it bears some similarities to the method of [START_REF] Barucq | New absorbing layers conditions for short water waves[END_REF], where a stable PML had been constructed for the shallow-water model. Indeed, this is not the first attempt of construction of stable perfectly matched layers in plasmas. Specifically, in [START_REF] Prokopidis | A higher-order spatial FDTD scheme with CFS PML for 3D numerical simulation of wave propagation in cold plasma[END_REF] the author considers a non-magnetized plasma model, which is dispersive, but not anisotropic. For this model the Bérenger's PML is stable (in the sense that the corresponding system does not exhibit an exponential growth in time). Similarly, in [START_REF] Sun | Perfectly matched layer absorbing boundary condition for nonlinear twofluid plasma equations[END_REF] the Bérenger's PML is applied to nonlinear nonmagnetized plasmas, where it is demonstrated to be stable. A complicated case of a general plasma model in a low-frequency regime which leads to so-called whistler modes, which are forward or backward propagating depending on the wave vector, had been studied in [START_REF] Chevalier | A pml utilizing k-vector information as applied to the whistler mode in a magnetized plasma, Antennas and Propagation[END_REF]. The authors suggested a new formulation for the attenuation of the such waves, namely, k-PML, however, its use requires an additional optimization of the reflection coefficient: this formulation is not a priori perfectly matched. A general cold plasma model (i.e. for arbitrary background magnetic field) had been considered in [START_REF] Jacquot | 2d and 3d modeling of wave propagation in cold magnetized plasma near the tore supra icrh antenna relying on the perfecly matched layer technique[END_REF]. There the difficulty of the application of the Bérenger's PML to this problem had been examined, however, no stable method had been proposed. This article is organized as follows. In Section 2 we formulate a simplified plasma model and discuss its properties. We first concentrate on a simplified two-dimensional model, and show that the conventional PMLs exhibit instabilities when applied to the cold plasma model in question. Next, we construct a stable PML based on the method of [START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF][START_REF] Cummer | Perfectly matched layer behavior in negative refractive index materials[END_REF]. This technique is well-adapted for models that have a frequency gap between forward and backward propagating modes. We demonstrate its stability with the help of theoretical and numerical arguments. This is a subject of Section 3. A part of this section had been briefly described in the conference proceedings [START_REF] Bécache | Perfectly matched layers in negative index metamaterials and plasmas[END_REF].

In Section 4 we study the three-dimensional model, which presents an additional complication, more precisely, the presence of backward and forward propagating modes for the same frequency. We demonstrate, however, that due to a special structure of the dispersion relation, such a system can be decomposed into two systems, one resembling the two-dimensional cold plasma model, and another behaving like the isotropic non-dispersive Maxwell system. As a result, we apply different perfectly matched layer techniques to each of these systems. Finally, we show how the coupling between the original model and the new systems can be done in a discrete setting (Section 5). We finish the exposition with the numerical experiments (Section 6) that confirm the stability of the new technique.

The Model of Cold Plasma in a Strong Background Magnetic Field

The Model

We consider the problem of wave propagation in cold plasmas comprised of particles of single species in a background magnetic field B 0 = (0, 0, B 0 ). This phenomenon is described by the Maxwell equations [22, Chapter 1] (where we use scaled coordinates in order that the speed of light c = 1)

curl curl E -ω 2 ε(ω)E = 0, (1)
where ε(ω) is a cold plasma dielectric tensor defined as [22, Chapter 1-2]

ε(ω) =       1 - ω 2 p ω 2 -ω 2 c i ω 2 p ωc ω(ω 2 -ω 2 c ) 0 -i ω 2 p ωc ω(ω 2 -ω 2 c ) 1 - ω 2 p ω 2 -ω 2 c 0 0 0 1 - ω 2 p ω 2       .
Here ω p is a plasma frequency, which depends on properties of particles, of which the plasma is comprised, and ω c = γB 0 is an algebraic cyclotron frequency, with γ > 0 depending on particle properties (charge and mass). We are interested in a regime when the background magnetic field becomes very large, B 0 → ∞ (thus ω c → ∞). In this case the cold plasma dielectric tensor becomes diagonal

ε(ω) =   1 0 0 0 1 0 0 0 1 - ω 2 p ω 2   . (2)
In the time domain, the model (1) associated with ε(ω) can be recast into the Maxwell equations coupled with an ODE for the current

∂ t E -curl B + Je z = 0, ∂ t B + curl E = 0, ∂ t J = ω 2 p E z , (3)
where e z = (0, 0, 1). We will assume that ω p is uniform in time and space. In a more general case, this should be true at least in a part of the physical domain close to the boundary, where we will apply the perfectly matched layer.

Properties of the Cold Plasma Model

The well-posedness of the system (3) equipped with appropriately smooth initial conditions follows from the Hille-Yosida theorem. Testing the first equation with E, the second equation with B and summing the results, we obtain, after integration by parts:

1 2 d dt R 3 |E| 2 dx + 1 2 d dt R 3 |H| 2 dx + R 3 JE z dx = 0,
or, alternatively, with the use of

E z = ω -2 p ∂ t J, d dt   1 2 R 3 |E| 2 + |H| 2 + ω -2 p J 2 dx   = 0.
One could expect that the conservation of the energy implies the non-dissipativity of the solutions of (3). In order to understand the behaviour of the solutions better, let us analyze (3) with the help of the Fourier techniques.

To do so, it would be advantageous to rewrite the system (3) in a second order formulation:

∂ tt E + ω 2 p (E • e z )e z + curl curl E = 0. ( 4
)
Plugging into the above a plane-wave solution Êe i(ωt-k•x) , with Ê ∈ R 3 , we obtain the following relation

A(k) Ê = ω 2 Ê, A(k) Ê = ω 2 p ( Ê • e z )e z -k × (k × Ê). (5)
Hence, (4) has a non-vanishing plane-wave solution if and only if ω 2 is an eigenvalue of the matrix A(k) (5), i.e. it solves the dispersion relation det ω 2 -A(k) = 0. We will adapt the following definition from ( [START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF]).

Definition 1. ( [START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF]) We will refer to a system as to a non-dissipative system if the solutions (or modes) ω j = ω j (k), j = 1, . . . , n, of its dispersion relation are real for all k ∈ R d .

From [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] we immediately obtain the following result.

Lemma 1. The system (4) is non-dissipative. Its modes are given by

ω 2 i = |k| 2 , ω 2 ± = ω 2 p + |k| 2 ± ∆(k) 2 , (6)
where the index i in ω i stands for 'isotropic'. The function ∆ : R 3 → R ≥0 is defined as

∆(k) = (ω 2 p + |k| 2 ) 2 -4k 2 z ω 2 p and satisfies (ω 2 p -|k| 2 ) 2 ≤ ∆(k) ≤ (ω 2 p + |k| 2 ) 2 for all k ∈ R 3 . (7)
Proof. Notice that the matrix A(k) is symmetric non-negatively definite:

(A(k) Ê, E ) = ω 2 p ( Ê • e z )(E • e z ) + ((k × Ê) • (k × E )),
hence all its eigenvalues are non-negative. Since the squares of the modes of the (4) are the eigenvalues of the matrix A(k), cf. ( 5), all the modes are real. Now let us obtain the eigenvalues of A(k) explicitly. Let us rewrite this matrix (using the relation k

× k × Ê = |k| 2 Ê -k(k • Ê)) A(k) Ê = |k| 2 E + B(k) Ê, B(k) Ê = -( Ê • k)k + ω 2 p ( Ê • e z )e z . (8)
From the above we see immediately that eigenvectors of A(k) coincide with the eigenvectors of B(k). Thus, the vector k × e z which belongs to Ker B(k) (since it is orthogonal to both k and e z ) is an eigenvector of A(k) associated to the eigenvalue λ = |k| 2 = ω 2 i (with index i standing for 'isotropic'). Since B(k) is symmetric, the rest of the eigenvectors Ê are orthogonal to k × e z (thus satisfy k x Êy -k y Êx = 0), and thus entirely determined by ( Ê • k) and ( Ê • e z ). Hence, using [START_REF] Hu | A stable, perfectly matched layer for linearized Euler equations in unsplit physical variables[END_REF], we obtain the following identities for such eigenvectors:

A(k) Ê • k = ω 2 ( Ê • k) = ω 2 p ( Ê • e z )k z , A(k) Ê • e z = ω 2 ( Ê • e z ) = |k| 2 ( Ê • e z ) -( Ê • k)k z + ω 2 p ( Ê • e z ).
From this, substituting (

Ê • k) in the second equation by ω -2 ω 2 p ( Ê • e z )k z , we obtain (ω 2 -|k| 2 -ω 2 p ) = -ω -2 k 2 z ω 2 p .
The solutions of the above equality are given by

ω 2 ± (k) = ω 2 p + |k| 2 ± ∆(k) 2 , ∆(k) = (ω 2 p + |k| 2 ) 2 -4k 2 z ω 2 p , with ∆(k) is a priori positive for all k ∈ R 3 , since ω ± ∈ R. One obtains immediately that ∆(k) ≤ (ω 2 p + |k| 2 ) 2 . Moreover, ∆(k) = (ω 2 p -|k| 2 ) 2 + 4(k 2 y + k 2 z )ω 2 p ≥ (ω 2 p -|k| 2 ) 2 ,
from which [START_REF] Tam | Perfectly matched layer as an absorbing boundary condition for the linearized Euler equations in open and ducted domains[END_REF] follows.

Thus, we have shown a non-dissipative character of the system (4).

Lemma 2. Assume k ∈ R 3 \ span{e z }. The eigenvector of the matrix A(k) associated to the eigenvalue ω 2 i = |k| 2 spans a one-dimensional space

Λ(k) = { Ê ∈ R 3 : k x Êx + k y Êy = 0, Êz = 0}.
The eigenvectors associated to the eigenvalues ω 2 + (k) and ω 2 -(k) span the space

Λ ⊥ (k) = { Ê ∈ R 3 : k y Êx -k x Êy = 0}. Indeed, R 3 = Λ(k) ⊕ Λ ⊥ (k).
The proof of this result is elementary and left to the reader (see also Lemma 6). Thus, we would expect plane waves associated with Λ to propagate as if they were in vacuum, and the plane waves associated with Λ ⊥ to be responsible for the dispersive effects (due to the form of the corresponding eigenvalues ω ± (k)).

The isotropic non-dispersive Maxwell system has an eigenvalue ω 2 = |k| 2 of multiplicity 2, and two modes ω = 0, thus having in total 6 modes, two out of which are non-propagative. It may seem surprising that there are actually 6 propagative and one non-propagative mode ω = 0 associated to the system (3). However, not all the modes are present in all frequency bands (we consider k ∈ R 3 ). This is detailed in the following result. Lemma 3. Let ω 2 + , ω 2 -be given by [START_REF] Bérenger | Three-dimensional perfectly matched layer for the absorption of electromagnetic waves[END_REF]. Then for all k ∈ R 3 \ {±ω p e z } it holds

ω 2 -(k) < ω 2 p < ω 2 + (k).
For k = ±ω p e z , the equality is achieved:

ω 2 -(±ω p e z ) = ω 2 p = ω 2 + (±ω p e z ).
Proof. Using ( 6) and ( 7), we obtain

ω 2 + = ω 2 p + |k| 2 + √ ∆ 2 ≥ ω 2 p + |k| 2 + |ω 2 p -|k| 2 | 2 ≥ ω 2 p , ω 2 -= ω 2 p + |k| 2 - √ ∆ 2 ≤ ω 2 p + |k| 2 -|ω 2 p -|k| 2 | 2 ≤ ω 2 p . Notice that ω 2 -(k) = ω 2 + (k) only when ∆ = 0. Also, ∆ = (ω 2 p + |k| 2 -2k z ω p )(ω 2 p + |k| 2 + 2k z ω p ) = ((ω p -k z ) 2 + k 2 x + k 2 y )((ω p + k z ) 2 + k 2 x + k 2 y ),
hence ∆ = 0 if and only if k = (0, 0, ±ω p ).

Hence, for any frequency ω < ω p and ω > ω p , there exist exactly 4 propagative modes, associated to ω 2 i and either ω 2 + or ω 2 -. For ω = ω p , the two modes ω 2 -, ω 2 + are equal.

Remark 1. The above results are valid for the first order formulation (3), which, however, has an additional mode ω = 0 compared to the second order formulation. This result can be checked by a direct computation of the dispersion relation of the system (3).

Reduced Two-Dimensional Problems

In order to construct stable perfectly matched layers for the three-dimensional problem (3), let us first start with simpler two-dimensional reduction of the original system. Similarly to the case of the Maxwell equations in vacuum, looking for solutions independent of a spatial variable decouples the original system (3) into two subsystems, namely transverse electric (TE), where the field B becomes scalar, and transverse magnetic (TM), where the field E is scalar. In particular, taking x = const, we obtain the (TE) x model

∂ t B x + ∂ y E z -∂ z E y = 0, ∂ t E y -∂ z B x = 0, ∂ t E z + ∂ y B x + J = 0, ∂ t J = ω 2 p E z . (TE) x
The corresponding (TM) x system coincides with the TM system for vacuum (as it does not involve E z ). The TE and TM equations for y = const are similar to (TE) x and (TM) x correspondingly, due to symmetry of the system (3) in (x, y).

For z = const, the (TE) z system coincides with the corresponding TE system in vacuum. The TM system for z = const, however, is different:

∂ t E z + ∂ y B x -∂ x B y + J = 0, ∂ t B x + ∂ y E z = 0, ∂ t B y -∂ x E z = 0, ∂ t J = ω 2 p E z .
(TM) z

The systems (TE) x and (TM) z differ from the 3D system (3) by the absence of the isotropic mode [START_REF] Bérenger | Three-dimensional perfectly matched layer for the absorption of electromagnetic waves[END_REF]. The following is obtained directly from Lemma 1.

Corollary 4. The dispersion relation of (TE) x is given by

F 2 (ω, k y , k z ) = ω 4 -ω 2 (ω 2 p -k 2 y -k 2 z ) + k 2 z ω 2 p = 0, (9)
and its solutions are ω 2 (k y , k z ) = ω 2 ± (0, k y , k z ), see [START_REF] Bérenger | Three-dimensional perfectly matched layer for the absorption of electromagnetic waves[END_REF]. For (TM) z the solutions are ω 2 (k x , k y ) = ω 2 ± (k x , k y , 0), more precisely ω 2 = 0 and

ω 2 = ω 2 p + k 2 x + k 2 y .
Next, we proceed as follows. First, we discuss the stability of the perfectly matched layers for the problem (TE) x , demonstrating with the help of the numerical experiment that the standard Bérenger's PML is unstable for it. We will explain this instability by applying the results of [START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF]. Next, based on the ideas of [START_REF] Cummer | Perfectly matched layer behavior in negative refractive index materials[END_REF] and [START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF], we will suggest a modified stable perfectly matched layer for the 2D system (TE) x . We will not concentrate on the problem (TM) z , since the Bérenger's PML does not lead to instabilities when applied to it. The idea of stabilizing the PML for (TE) x will serve us later as one of the components of the construction of stable perfectly matched layers in 3D.

Construction and Stability of PMLs for the System (TE) x

3.1. Bérenger's PML for the 2D Plasma System (TE) x Recall that geometrically the Bérenger's PML constitutes a layer around a box, inside which the physical unknowns are computed. The choice of the orientation of the box (and thus, of the layer) depends on a particular system of equations under consideration. For example, it is natural to align the edges of the box with e x , e y , e z axes when considering e.g. the isotropic Maxwell equations in vacuum. Here we make the same choice as well, noticing that the direction e z is 'privileged', since it coincides with the direction of the external magnetic field B z , hence it makes sense to orient one of the edges of the computational box along this direction. This choice significantly simplifies the construction of stable PMLs. To construct the Bérenger's perfectly matched layer in the direction e y , for y > 0, one performs the change of variables y → y + 1 iω y 0 σ(y )dx [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF] in the half-plane y > 0. Importantly, the absorption parameter σ(y ) > 0 for y > 0. Similarly a PML in other directions can be constructed. For a more detailed introduction into the perfectly matched layers see e.g. [START_REF] Joly | An elementary introduction to the construction and the analysis of perfectly matched layers for time domain wave propagation[END_REF].

In Figure 3.1 we demonstrate the results of the simulations obtained with the help of the perfectly matched layer based on the change of variables as per Bérenger. The discretization is done with the help of the Yee scheme [START_REF] Yee | Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media[END_REF] for dispersive Maxwell equations, which can be viewed as a coupling of the leapfrog discretization for the non-dispersive terms and the trapezoidal implicit discretization for the dispersive terms, see also [START_REF] Silva | Stable explicit coupling of the yee scheme with a linear current model in fluctuating magnetized plasmas[END_REF]. We choose zero initial conditions, and use a right-hand side perturbation f (t, y, z) = e -15(t-1.6) 2 -3y 2 -3z 2 of the first equation in the system (TE) x . The plasma frequency ω p equals 5.

In Figure 3.1 we can see that the instability occurs in the direction e y (this is clearly pronounced in the fourth figure), while in the direction e z the Cartesian PML seems to remain stable. This will be explained in Section 3.2.

Study of the Group Velocity

The questions of the stability of a PML system has been investigated in multiple works, see e.g. [START_REF] Bécache | On the analysis of Bérenger's perfectly matched layers for Maxwell's equations[END_REF][START_REF] Appelö | Perfectly matched layers for hyperbolic systems: general formulation, well-posedness, and stability[END_REF][START_REF] Abarbanel | Well-posed perfectly matched layers for advective acoustics[END_REF]. Our work is primarily based on the results of [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF] originally derived for anisotropic non-dispersive systems that has been recently generalized to dispersive non-dissipative systems in [START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF]. The results presented here are based on the analysis done for the constant absorption parameter σ(y) = σ ≥ 0, y ∈ R in [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF].

In particular, we will make use of the notion of the uniform stability, contrary to a weaker definition of stability of [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF]. After performing the PML change of variables in the frequency domain and coming back to the time domain, we obtain a perfectly matched layer system. For a moment, we do not discuss the details of passage from the frequency domain to the time domain.

Definition 2. ([11]

) The perfectly matched layer system is called uniformly stable if for all σ ≥ 0 and all k ∈ R 2 all its modes ω j (k) satisfy Im ω j (k) ≥ 0.

Given a mode ω(k), we call V ph (ω(k)) = ω(k) |k| k |k| the phase velocity, and V g (ω(k)) = ∇ k ω(k) the group velocity. Definition 3. (see [START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF][START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF]

) Given k ∈ R d , a mode ω(k) is called backward in the direction n ∈ R d if, for some k ∈ R d , (V ph (ω(k)) • n) (V g (ω(k)) • n) < 0.
If the contrary holds true (i.e. for all k ∈ R 3 the above product is nonnegative), such a mode is called forward. If the above product vanishes for all k ∈ R d , such a (forward) mode is called non-propagative.

The following result shows that Bérenger's (or traditional) perfectly matched layers are not uniformly stable in the presence of backward propagating waves.

Theorem 5. ([11, Proposition 4.8]) Let the system in the time domain associated to the PML system be uniformly stable in the sense of Definition 2. Then all the modes of the original system are forward (as in Definition 3) in e x .

Hence, we will consider the stability of the PMLs in two selected directions e y , e z . The following result shows that all the modes of the two-dimensional system (TE) x are forward propagating in the direction e z collinear to B 0 , and backward propagating in the direction e y . Lemma 6. The modes ±ω -(0, k y , k z ) of the system (TE) x are forward propagating in the direction e z and backward propagating in the direction e y for all k y , k z ∈ R. The modes ±ω + (0, k y , k z ) of the system (TE) x are forward propagating in the directions e y , e z .

Proof. Recall that ω 2 ± (0, k y , k z ) = ω 2 p + k 2 y + k 2 z ± ∆(0, k y , k z ) /2
, see Lemma 1. First we consider the direction e y . A straightforward computation gives:

∂ ky ω ± (0, k y , k z ) = 1 ω ± (0, k y , k z ) k y ± (ω 2 p + k 2 y + k 2 z )k y ∆(0, k y , k z )
.

With the use of the second inequality in [START_REF] Tam | Perfectly matched layer as an absorbing boundary condition for the linearized Euler equations in open and ducted domains[END_REF], we can see that

V ph,y (ω -)V g,y (ω -) = k y ω - k 2 y + k 2 z ∂ ky ω -(0, k y , k z ) = k 2 y 1 - ω 2 p + k 2 y + k 2 z ∆(0, k y , k z ) ≤ 0,
for all k y , k z ∈ R, and hence the mode ω -(0, k y , k z ) is backward propagating in the direction e y (and so is -ω -(0, k y , k z )). Similarly, by definition of the forward propagating mode, one sees that ω + (0, k y , k z ) is forward in the direction e y (and so is -ω + (0, k y , k z )). To consider the stability in the direction e z , we compute

∂ kz ω ± (0, k y , k z ) = 1 ω ± k z ± 4(ω 2 p + k 2 y + k 2 z )k z -8k z ω 2 p 4 ∆(0, k y , k z ) = k z ω ± ∆(0, k y , k z ) ∆(0, k y , k z ) ± (k 2 y + k 2 z -ω 2 p ) .
With the use of the first inequality in [START_REF] Tam | Perfectly matched layer as an absorbing boundary condition for the linearized Euler equations in open and ducted domains[END_REF], we can see that the modes ω ± (0, k y , k z ), -ω ± (0, k y , k z ) are forward propagating in the direction e z .

The above result, combined with Theorem 5, explains the results of the numerical experiments in Section 3.1, in particular, the fact that the instability seem to occur only in the direction e y , cf. Figure 3.1. However, forward and backward propagating modes are not present simultaneously, i.e. for the same frequency ω. This is a direct corollary of Lemmas 3 and 6. Namely, the modes ±ω + , which are forward propagating in the direction e y , and the modes ±ω -, which are backward propagating in e y , satisfy

ω 2 -(0, k y , k z ) < ω 2 p < ω 2 + (0, k y , k z ), for all (k y , k z ) ∈ R : (k y , k z ) = (0, ±ω p )
, and in the points (0, ±ω p ), it holds that ω 2 -= ω 2 p = ω 2 + . Thus, the following holds true for the modes ω of (TE) x :

V ph,y (ω)V g,y (ω) ≤ 0, ω < ω p , V ph,y (ω)V g,y (ω) ≥ 0, ω ≥ ω p . (11)
An illustration to this statement is shown in Figure 3.2. In order to understand the figure better, consider the dispersion relation ( 9) rewritten in the following form:

k 2 y 1 - ω 2 p ω 2 -1 + k 2 z -ω 2 = 0.
For ω < ω p , the level curves ω = const in the (k y , k z )-plane, as seen from the above, are hyperbolaes, and for ω > ω p , they are ellipses. The group velocity V g = ∇ k ω(k) is orthogonal to curves ω = const, and its actual direction can be computed as in Lemma 6. We can see that for ω < ω p (in the left) the projections of V g , V ph on this direction are of opposite signs, and for ω ≥ ω p they are of the same sign. In the following sections we discuss how to stabilize the perfectly matched layers for the system (TE) x , first considering the construction of the PMLs in different directions, and then proving their stability for the constant absorption parameter.

k y k z V g V ph k y k z V g V ph ω(k y , k z ) = ω * , ω * < ω p ω(k y , k z ) = ω, ω > ω p .

Construction of Stable PMLs

In Lemma 6 we have shown that all the modes are forward propagating in the direction e z , which is a necessary condition of the stability, as suggested by Theorem 5. In practice this condition is often sufficient, hence we would expect the Bérenger's PML to be stable in the direction e z . Hence, we use the Bérenger's change of variables in the direction e z (assuming that we need to construct the layer is located in the half-plane z > 0), namely

z → z + 1 iω z 0 σ(z )dz .
As for the direction e y , we will make use of the ideas of Bécache et al. [START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF], who themselves extended the work of Cummer [START_REF] Cummer | Perfectly matched layer behavior in negative refractive index materials[END_REF] on the Drude model. In [START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF] it was suggested to use a new change of variables that would stabilize the PML:

y → y + ψ(ω) iω y 0 σ(y )dy , (12)
cf. [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF]. Here, like before, we suppose that the PML is to be constructed in the half-plane y ≥ 0. It is suggested to choose the function ψ(ω) so that

ψ(ω)V ph,y (ω)V g,y (ω) ≥ 0. ( 13
)
Combining the above idea with [START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF] we immediately notice that choosing

ψ(ω) = 1 - ω 2 p ω 2 -1 (14)
in ( 12) may stabilize the PML, since this choice ensures that the product ( 13) is always non-negative.

The New PML System and its Stability

Let us apply the new PML change of variables to (TE) x . We perform the change of variables as described in Section 3.3, both in directions e z with the absorption parameter σ z and in e y with the absorption parameter σ y . There exists no single way to write the PML system in the time domain. The corresponding PDE formulations depend on the choice of the auxiliary unknowns, which can be done somewhat arbitrarily. This choice does not affect the original fields, however, impacts the complexity of the system. The choice we have done below is the one that appears to us as leading to a more tractable PML system. A detailed derivation of the corresponding equations can be found in Appendix A. Here we present the final system only (we introduced a new unknown j p = ω -2 p J, ∂ t j p = E z (this is done for convenience), and two PMLrelated unknowns B *

x , K * x ):

∂ t B x + σ z B x + B * x + ∂ y (E z + σ z j p ) -∂ z E y = 0, ∂ t E y + σ z E y -∂ z B x = 0, ∂ t E z + σ y E z + ω 2 p j p + ∂ y B x = 0, ∂ t j p = E z , ∂ t B * x + ω 2 p K * x + σ y B * x + σ y ∂ y (E z + σ z j p ) = 0, ∂ t K * x = B * x . (15)
Setting σ y = 0 (σ z = 0) while keeping σ z (σ y ) positive constant, we obtain from [START_REF] Abarbanel | Well-posed perfectly matched layers for advective acoustics[END_REF] the formulation of the PML in one direction, and the choice σ z , σ y > 0 corresponds to the corner PML.

Let us demonstrate that the PML system ( 15) is uniformly stable in the sense of Definition 2. Crucially, from now on we will assume that σ z , σ y are nonnegative constants. The dispersion relation of the system [START_REF] Abarbanel | Well-posed perfectly matched layers for advective acoustics[END_REF] reads (here we omit some tedious computations, presenting the final result)

1 + σ z iω 2 (ω 2 p -ω 2 ) 2 1 + σ y ψ(ω) iω 2 F pml (ω, k x , k y ) = 0, F pml (ω, k x , k y ) = k 2 z 1 + σ z iω -2 + k 2 y ψ(ω) 1 + σ y iω ψ(ω) -2 -ω 2 ,
with ψ(ω) as in [START_REF] Nataf | New constructions of perfectly matched layers for the linearized Euler equations[END_REF]. Notice that F pml can be obtained from F 2 , see ( 9) by a change of variables k y → k y 1 + σ y ψ(ω)(iω) -1 -1 and k z → k z 1 + σ z (iω) -1 -1 . We will show that the above equation cannot have solutions ω, s.t. Im ω < 0, for any σ z , σ y ≥ 0, thus answering altogether the question of the stability of the PML in one direction (when either σ z or σ y vanish) and of the corner PML (σ z , σ y > 0). Before formulating the stability result, let us present the following auxiliary lemma. Let us denote C -:= {ω ∈ C| Im ω < 0}.

Lemma 7. Let α(ω) : C -→ C. Assume Im (ωα(ω)) > 0 for ω ∈ C -.
Then, for all σ ≥ 0 and all ω ∈ C -, it holds

Im ωα(ω) 1 + σ iω α(ω) -2 > 0. ( 16
)
Proof. To see [START_REF] Diaz | Stabilized perfectly matched layer for advective acoustics[END_REF], notice that Im ω α(ω)z -1 = |z| -2 Im (ω α(ω)z) for z ∈ C, and hence

sign Im ω α(ω) 1 + σ iω α(ω) -2 = sign Im ω α(ω) 1 - σ iω α(ω) 2 = sign Im ω α(ω) + 2iσ|α(ω)| 2 -σ 2 α(ω)|α(ω)| 2 ω .
The term Im (ω α(ω)) > 0, for all ω ∈ C -. The term Im 2iσ|α(ω)| 2 ≥ 0 (since σ ≥ 0). Finally, for all ω ∈ C -, the expression sign Im -ω -1 α(ω) = sign Im α(ω)ω -1 = sign Im (ωα(ω)), from which [START_REF] Diaz | Stabilized perfectly matched layer for advective acoustics[END_REF] follows.

Now we have all the means to formulate the result of the stability of the PML.

Theorem 8. For all σ y , σ z ≥ 0, k y , k z ∈ R, the solutions of

1 + (iω) -1 σ z 1 + (iω) -1 ψ(ω)σ y = 0 (a) F pml (ω, k x , k y ) = 0 (b) lie in C \ C -.
As a consequence, the PML system ( 15) is uniformly stable.

Proof. Let us first prove the statement for the equation (a). The solutions of (a) are either ω = iσ z , with Im ω ≥ 0 as soon as σ z ≥ 0, or satisfy 1 + σ y ψ(ω)(iω) -1 = 0. Instead of using an explicit expression for ψ(ω) to show that the solutions to the above equation are located in the upper half-plane, we suggest to employ a certain property of ψ(ω). Let us show that ψ(ω) satisfies

Im (ωψ(ω)) > 0, for all ω ∈ C -. (17)
First of all, sign Im (ωψ

(ω)) = sign Im ωψ(ω) -1 |ψ(ω)| 2 = sign Im ωψ(ω) -1 .
Thus, using an explicit expression of ψ(ω), see [START_REF] Appelö | A general perfectly matched layer model for hyperbolic-parabolic systems[END_REF],

sign Im (ωψ(ω)) = sign Im ω 1 - ω 2 p ω 2 = sign Im ω -ω 2 p ω|ω| -2 > 0
for all ω ∈ C -. Thus, Im 1 + (iω) -1 ψ(ω) does not vanish in C -:

Im 1 + σ y ψ(ω) ω = σ y |ω| -2 Im (ψ(ω)ω) > 0, σ y > 0.
Hence, (a) has no solutions in C -. To prove the second part of the statement, we multiply (b) by ω. This gives

F 1 (ω) + F 2 (ω) + F 3 (ω) = 0, F 1 (ω) = -|ω| 2 ω, (18)
F 2 (ω) = ωk 2 z (1 + (iω) -1 σ z ) 2 , F 3 (ω) = ωk 2 y ψ(ω) (1 + (iω) -1 σ y ψ(ω))
2 .

We will show that [START_REF] Barucq | New absorbing layers conditions for short water waves[END_REF] has no solutions in C -, since for all ω ∈ C -the imaginary part of the left-hand side of the above expression is always strictly positive, which is a consequence of the fact Im F j > 0, j = 1, 2, 3. More precisely,

• F 1 (ω) = -Im(ω|ω| 2 ) > 0, for ω ∈ C -;
• to see that Im F 2 (ω) > 0 for ω ∈ C -, use Lemma 7 with α(ω) = 1.

• to see that Im F 3 (ω) > 0 in C -, use Lemma 7 with α(ω) = ψ(ω). Indeed, it satisfies the condition of Lemma 7, as shown in [START_REF] Demaldent | Perfectly matched transmission problem with absorbing layers: application to anisotropic acoustics in convex polygonal domains[END_REF]. With the use of this result, we see that Im F 3 (ω) > 0 in C -, for all σ y ≥ 0.

The above lemma proves the uniform stability of the PML system [START_REF] Abarbanel | Well-posed perfectly matched layers for advective acoustics[END_REF]. An extension of this result to a more general class of systems can be found in [START_REF] Bécache | Stable perfectly matched layers for a class of anisotropic dispersive models[END_REF], while the analysis of the stability of the PML system with the help of the energy techniques has been performed in [START_REF] Kachanovska | Stable perfectly matched layers for a class of anisotropic dispersive models. Part II: Energy Estimates and Stability in Corners[END_REF].

Remark 2. Unfortunately, we were not able to obtain the result for varying σ y , σ z (which is not known in general, cf. e.g. [START_REF] Bécache | On the analysis of Bérenger's perfectly matched layers for Maxwell's equations[END_REF][START_REF] Diaz | A time domain analysis of PML models in acoustics[END_REF][START_REF] Halpern | The analysis of matched layers[END_REF][START_REF] Halpern | Bérenger/Maxwell with discontinous absorptions: existence, perfection, and no loss[END_REF][START_REF] Chen | Long-time stability and convergence of the uniaxial perfectly matched layer method for time-domain acoustic scattering problems[END_REF]).

Illustration of the Stability of the New PML System

To demonstrate that in practice the suggested change of variables leads to a stable PML formulation, we conduct a numerical experiment with the same parameters as in Figure 3

Plasma in Three Dimensions

To analyze the stability of the perfectly matched layers for the system (3) (or, equivalently, (4)), let us extend the result of Lemma 6 to the 3D system, making use of Lemmas 3 and 1. The proof of the result that follows repeats almost verbatim the proof of Lemma 6. ). With arrows we show computed directions of group and phase velocities. We can see that unlike in the 2D case, now, for a fixed frequency ω < ω p , there are both forward and backward modes present. Hence, it is impossible to apply the results of the work [START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF] directly, see Section 3.3. More precisely, choosing ψ(ω) to ensure that the product (13) for the mode ω -is positive would cause the corresponding product for the mode ω i to be negative, and vice versa. In this section we suggest a strategy to overcome this issue. Namely, the dispersion relation of ( 4) is a product of the dispersion relations of two systems. One of these systems resembles 3D Maxwell equations in vacuum (and hence, the Bérenger's PML is stable for this problem), and another one behaves like the 2D plasma model, thus we can apply to it the PML suggested in Section 3.

k ⊥ k z V g V ph V ph V g k ⊥ k z V ph , V g V g V ph ω(k) = ω * , ω * < ω p ω(k) = ω, ω > ω p .

Splitting of the Maxwell System

Decomposition of Fields

To derive a new, split system, we first analyze the eigenspaces corresponding to the modes of the system. Recall that a plane wave Ê e i(ωt-k•x) solves (4) if and only if Ê ∈ C 3 is an eigenvector corresponding to the eigenvalue ω 2 of the matrix A(k), see [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]. As already discussed in Lemma 2, the eigenvectors of the matrix A(k) define a decomposition of the space R 3 = Λ(k) ⊕ Λ ⊥ (k). Here we would like to extend this construction to arbitrary square-integrable fields.

Let us introduce the following notation. Given a vector a = (a x , a y , a z ) T ∈ R 3 , a scalar-valued function f : R 3 → R, and a vector-valued function F : R 3 → R 2 , let us introduce the following operators:

a ⊥ = a x a y , curl ⊥ f = ∂ y f -∂ x f , ∇ ⊥ f = ∂ x f ∂ y f , curl ⊥ F = ∂ x F y -∂ y F x , ∇ ⊥ • F = ∂ x F x + ∂ y F y .
We will need the following auxiliary spaces:

H d,0 = {f ∈ L 2 (R 3 ) 2 : ∇ ⊥ • f = 0}, H r,0 = {f ∈ L 2 (R 3 ) 2 : curl ⊥ f = 0}. ( 19
)
Lemma 10. The space L 2 (R 3 ) 2 can be decomposed as follows:

L 2 (R 3 ) 2 = H d,0 ⊕ H r,0 .
The spaces H d,0 , H r,0 are orthogonal with respect to the L 2 (R 3 ) 2 scalar product.

A proof of this result follows trivially from Lemma 2 (it suffices to consider the corresponding projectors P Λ , P ⊥ Λ in the Fourier space L 2 (R 3 ) 2 , and then apply the result to the first two components of the vector from L 2 (R 3 ) 2 ).

Remark 3. Notice that the decomposition of Lemma 10 is different from the Helmholtz decomposition in 2D or 3D, but can be viewed as an anisotropic two-dimensional Helmholtz decomposition applied to fields R 3 → R 2 .

Derivation of the Split System

Now we have all the necessary ingredients to introduce the splitting of the original system (3) into two subsystems, one of which will behave as the Maxwell equations in vacuum, and another one will resemble the 2D system (TE) x .

First of all, let us rewrite (3) with the help of the new notation introduced in the previous section. For this let us introduce an auxiliary matrix

α = 0 1 -1 0 . ( 20
)
Remark 4. Notice that for f : R 3 → R and F : R 3 → R 2 :

α∇ ⊥ f = -curl ⊥ f, α curl ⊥ f = ∇ ⊥ f, ∇ ⊥ (αF) = curl ⊥ F, curl ⊥ (αF) = -∇ ⊥ • F.
Then, the original Maxwell system (3) reads

∂ t E ⊥ + α ∂ z B ⊥ -curl ⊥ B z = 0, ∂ t E z -curl ⊥ B ⊥ + J = 0, ∂ t J = ω 2 p E z , ∂ t B ⊥ -α ∂ z E ⊥ + curl ⊥ E z = 0, ∂ t B z + curl ⊥ B ⊥ = 0. (M)
Let us introduce

E ⊥ = E d + E r , B ⊥ = B d + B r , with E r , B r ∈ H r,0 , E d , B d ∈ H d,0 , (21)
where H d,0 , H r,0 are defined in [START_REF] Prokopidis | A higher-order spatial FDTD scheme with CFS PML for 3D numerical simulation of wave propagation in cold plasma[END_REF]. Inserting the above into (M) and using curl ⊥ B r = curl ⊥ E r = 0,

(∂ t E d + α ∂ z B r -curl ⊥ B z ) + (∂ t E r + α ∂ z B d ) = 0, ( 22a
)
∂ t E z -curl ⊥ B d + J = 0, ( 22b
)
∂ t J = ω 2 p E z , (22c)
(∂ t B d -α ∂ z E r + curl ⊥ B z ) + (∂ t B r -α ∂ z E d ) = 0, ( 22d
)
∂ t B z + curl ⊥ E d = 0. ( 22e
)
One notices that provided E d , B d ∈ H d,0 and E r , B r ∈ H r,0 , the first term in (22a) belongs to the space H d,0 , and the second term stays in H r,0 . The same is valid for (22d). Thanks to the orthogonality of these spaces, the above leads to two systems of equations:

∂ t E d + α ∂ z B r -curl ⊥ B z = 0, ∂ t B r -α ∂ z E d = 0, ∂ t B z + curl ⊥ E d = 0. (TE)
and

∂ t B d -α ∂ z E r + curl ⊥ E z = 0, ∂ t E r + α ∂ z B d = 0, ∂ t E z -curl ⊥ B d + J = 0, ∂ t J = ω 2 p E z .
(TM)

We call the latter system (TM), since it does not have a B z -component of the magnetic field. It is easy to see that summing (TE) and (TM), we obtain (M) provided that curl ⊥ B r = curl ⊥ E r = 0.

Remark 5. In particular, for the system (TM), it is easy to notice that provided the initial data (E d0 , B r0 , E z0 , J 0 ), where

E d0 ∈ H d,0 , B r0 ∈ H r,0 , the solution E r (t) ∈ H r,0 , B d (t) ∈ H d,0 for all t ≥ 0.
To see this, it is sufficient to take ∇ ⊥ • of the first equation, and curl ⊥ of the second equation, which would give, thank to the choice of the right-hand side data:

∂ t ∇ ⊥ • B d -∂ z curl ⊥ E g = 0, ∂ t curl ⊥ E r -∂ z ∇ ⊥ • B d = 0.
This gives a 1D wave equation with zero initial and right-hand side data, hence the conclusion. A similar result holds for the system (TE).

We will discuss the well-posedness of (TE) and (TM) slightly later, first formulating a simplified result of the equivalence of these systems to the original Maxwell system, which follows straightforwardly from the above construction.

Theorem 11. Let (E, B, J) solve (M) with sufficiently smooth initial data (E 0 , B 0 , J 0 ). Let the decomposition of the initial data as per Lemma 10 be given by

E 0⊥ = E d0 + E r0 , B 0⊥ = B d0 + B r0 , E d0 , B d0 ∈ H d,0 , E r0 , B r0 ∈ H r,0 , (23)
Then

E = E r + E d Ẽz , B = B r + B d Bz , J = J, (24)
where (B d , E r , Ẽz , J) solves (TM) with the initial data (B d0 , E r0 , E z0 , J 0 ), and (E d , B r , Bz ) solves (TE) with the initial data (E d0 , B r0 , B z0 ).

Reciprocally, suppose that (B d , E r , Ẽz , J) solves (TM) with the initial data (B d0 , E r0 , E z0 , J 0 ) and (E d , B r , Bz ) solve (TE) with the initial data (E d0 , B r0 , B z0 ), and let (23) hold true. Then (E, B, J) defined as in ( 24) solves (M) with the initial data defined by [START_REF] Joly | An elementary introduction to the construction and the analysis of perfectly matched layers for time domain wave propagation[END_REF] and identities E z0 = Ẽz0 , J 0 = J0 , B z0 = Bz0 .

The well-posedness of the systems (TM), (TE) follows from the Hille-Yosida theory. Moreover, it is possible to show a conservation of a certain energy. Namely, after the integration by parts, we obtain the following energy identities for the systems (TE) and (TM) correspondingly:

1 2 d dt R 3 |B d | 2 + |E z | 2 + |E r | 2 + ω -2 p J 2 dx = 0, 1 2
d dt R 3 |E d | 2 + |B z | 2 + |B r | 2 dx = 0.

Dispersive Properties of the Systems (TE) and (TM)

Here we would like to show that the main goal of the splitting of the original Maxwell system had been successfully achieved: (TE) has only forward propagating modes, and the system (TM) has both backward and forward propagating modes, but they are not present simultaneously for a fixed frequency.

Let us first perform a plane-wave analysis of the system (TE). One can verify that the dispersion relation of this system reads:

F T E (ω, k) = ω(ω 2 -k 2 z )(ω 2 -|k| 2 ) = 0. ( 25
)
The corresponding modes are ω ≡ 0, ω 2 z (k) = k 2 z , ω 2 i (k). They are forward (or non-propagative, cf. Definition 3) in the directions e x , e y , e z .

Concerning the system (TM), its dispersion relation reads

F T M (ω, k) = ω 2 -k 2 z (ω 2 -ω 2 p )k 2 z + ω 2 (k 2 x + k 2 y ) -(ω 2 -ω 2 p )ω 2 = 0. ( 26
)
The solutions of ( 26) are ω 2 z (k) = k 2 z , ω 2 ± (k). Using Lemmas 6, 3, 9, we observe:

1. for |ω| < |ω p | there are 4 modes present: ω 2 z and ω 2 -. The modes ±ω z (k) are non-propagative in the directions e x , e y , and the modes ±ω -(k) are backward for all k ∈ R 3 . All these modes are forward in e z .

2. for |ω| ≥ |ω p | there are 4 modes present: ω 2 z (k) and ω 2 + (k). These modes are either non-propagative or forward in the directions e x , e y , e z .

Construction of Stable PMLs for the Systems (TE) and (TM)

The system (TE) has only forward modes, therefore, we will apply to it the Bérenger's perfectly matched layer change of variables. Like before, we will apply the changes of variables in all directions; to obtain the PML system in one direction it suffices to set the rest of the absorption parameters to zero. us present the corresponding PML system whose derivation is given in Appendix B. For convenience, we introduce the PML-related unknowns E * d , K z and B * zx , B * zy that serve as a perturbation of the original system

∂ t E d + σ z E d + E * d -curl ⊥ (B z + σ z K z ) + α ∂ z B r = 0, ∂ t B r + σ z B r -α∂ z E d = 0, ∂ t B z + B * zx + B * zy + curl ⊥ E d = 0,
∂ t E * dx + σ y E * dx -σ y ∂ y (B z + σ z K z ) = 0, ∂ t E * dy + σ x E * dy + σ x ∂ x (B z + σ z K z ) = 0, ∂ t B * zx + σ y B * zx + σ y ∂ y E dx = 0, ∂ t B * zy + σ x B * zy -σ x ∂ x E dy = 0, ∂ t K z = B z .
The stability of the above system can be proved as it was done in Section 3.4. As for the system (TM), the analysis in Section 4.1.3 shows that it behaves similarly to the 2D plasma model, see Section 3. Summarizing the said-above, for the modes ω of (TE) the following holds true:

V ph,β (ω)V g,β (ω) ≤ 0, ω < ω p , V ph,β (ω)V g,β (ω) ≥ 0, ω ≥ ω p , β ∈ {x, y}, V ph,z (ω)V g,z (ω) ≥ 0, for all ω. (27)
Hence, to stabilize the PML in the directions e x , e y we employ the same idea as in Section 3.3, and in the direction e z we use the Bérenger's PML. More precisely, we use the following changes of variables:

x → x + ψ(ω) iω x 0 σ(x )dx , y → y + ψ(ω) iω y 0 σ(y )dy , where ψ(ω) = 1 - ω 2 p ω 2 -1 = 1 + ω 2 p ω 2 -ω 2 p , z → z + 1 iω z 0 σ(z )dz . (28)
Let us now formulate the PML system and write the dispersion relation corresponding to it. Like before, for convenience, we introduced a new variable j p = ω -2 p J, ∂ t j p = E z , as well as several unknowns defined only inside the PML layer, namely B * d = (B * dx , B * dy ), E * zx , E * zy , and that can be viewed as a perturbation of the original system:

∂ t B d + σ z B d + B * d + curl ⊥ (E z + σ z j p ) -α ∂ z E r = 0, ∂ t E r + σ z r + α∂ z B d = 0, ∂ t E z + E * zx + E * zy -curl ⊥ B d + ω 2 p j p = 0, ∂ t j p = E z .
Indeed, we couple the above equations with the relations for the PML unknowns

B * d = (B * dx , B * dy ), E * zx , E * zy .
For convenience, we additionally introduce K * dx , K * dy , j * px , j * py , which are connected to the dispersive nature of the new PML change of variables (28):

∂ t B * dx + ω 2 p K * dx + σ y B dx + σ y ∂ y (E z + σ z j p ) = 0, ∂ t K * dx = B * dx , ∂ t B * dy + ω 2 p K * dy + σ x B dy -σ x ∂ x (E z + σ z j p ) = 0, ∂ t K * dy = B * dy , ∂ t E * zx + ω 2 p j * px + σ y E * zx -σ y ∂ y B dx = 0, ∂ t j * px = E * zx , ∂ t E * zy + ω 2 p j * py + σ x E * zy + σ x ∂ x B dy = 0, ∂ t j * py = E * zy .
The derivation of this system is given in Appendix B. To show its well-posedness and stability, we will use the modal approach, see e.g. [START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF], which consists of showing that the PML system is uniformly stable in the sense of Definition 2.

We essentially repeat the procedure of Section 3.4.

The dispersion relation of the new PML system reads, see also ( 26)

ω 4 ψ(ω) 4 1 + σ x ψ(ω) iω 2 1 + σ y ψ(ω) iω 2 (iω + σ z ) 4 F T M pml (ω, k) = 0, F T M pml (ω, k) = F T M ω, iωk x iω + σ x ψ(ω) , iωk y iω + σ y ψ(ω) , iωk z iω + σ z . (29)
Repeating verbatim the arguments of Lemma 8, we see that all the solutions to the above lie in the upper half-plane C \ C -.

Lemma 12. The equation ( 29) has no solutions ω ∈ C -.

A General Formulation of the Problem

While the splitting procedure we have introduced before indeed allows to extend the original PML construction for the system (TE) x , it is heavily based on the fact that the original Maxwell system (3) is homogeneous. The application we have in mind is the wave propagation in the media that is not necessarily homogeneous (cf. e.g. [START_REF] Silva | Stable explicit coupling of the yee scheme with a linear current model in fluctuating magnetized plasmas[END_REF]), but all its inhomogeneities are contained in a domain Ω 0 ⊂ Ω, where Ω is a cuboid. Thus, the Cartesian PML can be applied in R 3 \ Ω. It is therefore reasonable to use the Maxwell system in Ω, while outside Ω Ω s B 0 = (0, 0, B 0 ) of this domain perform the splitting, in order to apply the PML. For homogeneous applications, computationally it is cheaper to use the original system (3), which has seven unknowns, rather than the split system (TE)-(TM) for eleven unknowns.

y x z Ω Ω s Ω s Ω pml x z
We surround the cuboid Ω (which we will refer to as a physical domain) by a layer Ω s (the index s stands for 'split'), which, in turn, is surrounded by Ω pml . These domains are defined as follows:

Ω = (-L x , L x ) × (-L y , L y ) × (-L z , L z ), (30)
Ω s = (-L s x , L s x ) × (-L s y , L s y ) × (-L z , L z ) \ Ω, Ω pml = (-L pml x , L pml x ) × (-L pml y , L pml y ) × (-L pml z , L pml z ) \ Ω ∪ Ω s ,
where , L pml β > L s β > L β > 0, for β ∈ {x, y} and L pml z > L z > 0. This geometrical setting is shown in Figure 4.2.

Inside the physical domain Ω we solve the Maxwell equations (necessarily, we assume that the right hand side and initial data are supported inside Ω), and inside Ω s we solve only the split system. These formulations are coupled through an interface between Ω and Ω s . To the split formulation which we solve in Ω s we apply the PML in Ω pml . Notice that we do not perform splitting in the direction z, since the Bérenger's PML is stable in this case, and we simply apply it to the original Maxwell system.

The derivation of a variational formulation for such a coupled system is not entirely trivial. The transmission conditions through the interface Γ s between Ω s and Ω read (denoting by n the normal pointing from Ω into Ω s , and by

E s r , E s d , E s z , B s r , B s d , B s
z the solution of the split system):

E × n = E s r + E s d E s z × n, B × n = B s r + B s d B s z × n on Γ s . (TC)
For instance, on parts of Γ s where n = (1, 0, 0):

E z = E s z , B z = B s z , E y = E dy + E ry , B y = B dy + B ry .
The difficulty in the derivation of the variational formulation for (M, TE, TM) with the transmission conditions (TC) lies in the fact that there is some asymmetry in the information required by (TE,TM) and (M). More precisely, the equation for E z in (M) requires knowing tangential components of B d + B r , while the equation for E z in (TM) requires the knowledge of the tangential components of B d . We were not able to derive a stable variational formulation in a continuous setting for the system (M, TE, TM)-(TC). That is why we chose to use the Yee scheme, for which we were able to derive a stable coupling of these systems in Section 5.2.

Discretization and Numerical Implementation

The Yee Scheme: Application to the Split System

What follows is essentially standard, and constitutes the most natural extension of the Yee scheme [START_REF] Yee | Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media[END_REF] for the dispersive Maxwell equations. Moreover, we introduce some new notation to allow for a compact presentation of the scheme. Without loss of generality, assume that the spatial mesh size is the same, h, in all directions. For any i ∈ Z, we set î = i + 1/2. We introduce the spaces of discrete vector-valued 2D fields (denoting v ijk = v(ih, jh, kh) and so on, ...):

V e⊥ = E ⊥ = (E x ) î,j,k , (E y ) i,,k , i, j, k ∈ Z 3 , V b⊥ = B ⊥ = (B x ) i,, k , (B y ) î,j, k , i, j, k ∈ Z 3 , (31)
and discrete scalar-valued functions

V ez = V i,j, k, i, j, k ∈ Z 3 , V bz = V î,,k , i, j, k ∈ Z 3 . ( 32
)
In the following, concerning the electric fields

• the total field E is discretized in V e⊥ × V ez • the fields E d , E r are discretized in V e⊥ .
Similarly, concerning the magnetic fields

• the total field B is discretized in V b⊥ × V bz • the fields B d , B r are discretized in V b⊥ .
This is illustrated in Figure 5.1. Naturally, the current J is discretized in the same space V ez as E z , which follows naturally from the original Maxwell system with the plasma term [START_REF] Jacquot | Description non lineáeire auto-cohérente de la propagation d'ondes radiofréquences et de la périphérie d'un plasma magnétisé[END_REF].

Let us now introduce discrete spatial differential operators. The discrete curl h ⊥ operator approximates the continuous operator curl ⊥ , and is defined by Similarly it's defined for fields U ∈ V b⊥ . In the same way, the discrete curl h ⊥ operator is defined as an approximation curl ⊥

curl h ⊥ ∈ L (V e⊥ ; V bz ) and curl h ⊥ ∈ L (V b⊥ ; V ez ) . For instance, for V ∈ V e⊥ , curl h ⊥ V î,,k = (V y ) i+1,,k -(V y ) i,,k h - (V x ) î,j+1,k -(V x ) î,j,k h . E y E ry E dy E x E rx E dx E z E z B z B z B x B dx B rx B y B dy B ry y x z
curl h ⊥ ∈ L (V bz ; V e⊥ ) and curl h ⊥ ∈ L (V ez ; V b⊥ ) .
We omit its natural definition. Finally, the discrete operator α ∂ h z approximates α ∂ z and is defined for the fields

α ∂ h z ∈ L(V b⊥ , V e⊥ ) and α ∂ h z ∈ L(V e⊥ , V b⊥ ).
For convenience, we use the same notation for operators defined on spaces associated to fields E (e.g. V e⊥ ) and B.

Concerning the time discretization, we use a constant time step ∆t > 0, discretize the electric field and the magnetic field at instants (n + 1/2)∆t. The discrete time differentiation operators are defined as follows:

∂ ∆t t E n+ 1 2 = E n+1 -E n ∆t , ∂ ∆t t B n = B n+ 1 2 -B n-1 2 ∆t .
As for the current J, we choose to discretize it at instants n∆t. We will need the following averaging operator

Jn+ 1 2 = 1 2 J n+1 + J n .
A natural generalization of the Yee scheme to the dispersive Maxwell system reads (notice that the case ω p = 0 reduces the scheme to the Yee scheme for the Maxwell system in vacuum):

∂ ∆t t E ⊥ n+ 1 2 + α ∂ h z B n+ 1 2 ⊥ -curl h ⊥ B n+ 1 2 z = 0, (33a)
∂ ∆t t E z n+ 1 2 -curl h ⊥ B n+ 1 2 ⊥ + J n+ 1 2 = 0, ∂ ∆t t J n+ 1 2 = ω 2 p E n+ 1 2 , (33b)
∂ ∆t t B ⊥ n+1 -α ∂ h z E n+1 ⊥ + curl h ⊥ E n+1 z = 0, (33c)
∂ ∆t t B z n+1 + curl h ⊥ E n ⊥ = 0. ( 33d
)
Notice that in the above scheme we use the (implicit) trapezoid rule to discretize the dispersive term, which allows to preserve the CFL stability condition of the Yee scheme for classical Maxwell equations, see Remark 7. Note, however, that the scheme remains in practice fully explicit, see [START_REF] Silva | Stable explicit coupling of the yee scheme with a linear current model in fluctuating magnetized plasmas[END_REF]. Similarly, the discretization of (TE) reads

∂ ∆t t E d n+ 1 2 + α ∂ h z B n+ 1 2 r -curl h ⊥ B n+ 1 2 z = 0, (34a)
∂ ∆t t B r n+1 -α ∂ h z E n+1 d = 0, (34b)
∂ ∆t t B z n+1 + curl h ⊥ E n d = 0. (34c)
Finally, we discretize (TM):

∂ ∆t t B d n+1 -α ∂ h z E n+1 r + curl h ⊥ E n+1 z = 0, (35a)
∂ ∆t t E r n+ 1 2 + α ∂ h z B n+ 1 2 d = 0, (35b)
∂ ∆t t E z n+ 1 2 -curl h ⊥ B n+ 1 2 d + J n+ 1 2 = 0, ∂ ∆t t J n+ 1 2 = ω 2 p E n+ 1 2 . (35c)

A Discrete Coupling Procedure

As we have already discussed in Section 4.3, the coupling of the split system (TE, TM) and (M) is done on a discrete rather than continuous level. In this section we will study this question in more detail. We will not present the details of the coupling between the split system and the perfectly matched layer, and refer an interested reader to e.g. [START_REF] Joly | An elementary introduction to the construction and the analysis of perfectly matched layers for time domain wave propagation[END_REF].

In the domain Ω, we solve the Maxwell system of equations, while in the domain Ω s we discretize (T E, T M ). An additional coupling is needed at the interface for the unknowns E z and B z . This can be seen in Figure 5.2: indeed, in order to compute i = Nx , (B z )

n+ 1 2
Nx,,k using the equation of the split system (TE), one requires the values

(E dy ) n Nx+1,,k , ( E dy )
n Nx,,k , the latter of which is not defined. The unknowns E, B, J are defined in the spaces associated to the mesh on the set Ω, and the unknowns E r , E d , E z , J, as well as B r , B d , B z are defined on the spaces associated to the mesh inside Ω s . For the case of finite index sets, we additionally impose the boundary data

E × n = 0, B • n = 0 on ∂Ω \ (∂Ω ∩ ∂Ω s ) , (E r , E z ) t × n = 0, (E d , 0) t × n = 0, (B r , B z ) t • n = 0, (B d , 0) t • n = 0 on ∂Ω s \ (∂Ω ∩ ∂Ω s ) . ( 36
)
The discretization is done with the help of [START_REF] Bidégaray-Fesquet | Stability of FD-TD schemes for Maxwell-Debye and Maxwell-Lorentz Equations[END_REF][START_REF] Collino | Optimizing the perfectly matched layer[END_REF]35). The interface Γ s is composed of the edges of the box Ω that have normals perpendicular to e z . It can be characterized by the index set In order to obtain the coupling conditions for B z , we surround Ω by a thick 1-cell interface (cf. î,,k , for a special case when |j| < N y and i > 0:

I s = (i, j) : |i| = N x , |j| ≤ N y ∪ (i, j) : |i| ≤ N x , |j| = N y . E y E y E x E x E z E z E z E dy E dx E dx E ry E rx E rx B n+ 1 2 z = B n-1 2 z -∆t ∂ h x E n y -∂ h y E n x B n+ 1 2 z = B n-1 2 z -∆t ∂ h x E n dy -∂ h y E n dx B n+ 1 2 z = B n-1 2 z -∆t ∂ h x (E n dy + E n ry ) -∂ h y (E n dx + E n rx ) Ω Ω s Thick interface
• for i > N x , inside Ω s , use (34c); • for i < N x , inside Ω, use (33d);
• for i = N x , use a modified equation (see the large arrow in Figure 5.2):

∂ ∆t t B z n Nx,,k + curl h ⊥ Ẽn ⊥ Nx,,k = 0, (37)
where Ẽ⊥ = ( Ẽx , Ẽy ) with Ẽx = E rx + E dx and Ẽy i,,k

= (E ry + E dy ) i,,k , i = N x + 1, (E y ) i,,k , i = N x .
Similarly, to compute E z n+1 ij k , we use the following procedure (described for the case when |j| < N y and i > 0), cf. Figure 5.3:

• for i > N x , inside Ω s , use (35c); • for i < N x , inside Ω, use (33b);
• for i = N x , i.e. on the interface, use a modified equation (see the thick arrow in Figure 5.3)

∂ ∆t t E z n+ 1 2 Nx,j, k + Jn+ 1 2 Nx,j, k = curl h ⊥ Bn+ 1 2 ⊥ Nx,j, k , (38)
where B⊥ = ( Bx , By ), with Bx = B x and In the rest of the cells that are close to the interface Γ s , the unknowns E z and B z are computed similarly. Moreover, as we can see in Figure 5.3, no special coupling equations should be introduced for unknowns E ⊥ (and B ⊥ ).

By Nx,j, k = (B ry + B dy ) î,j, k , i = N x , (B y ) î,j, k , i = N x -1. B x B x B dy B ry B y E y B z B z Ω Ω s E n z + Jn+ 1 2 = E n z + ∆t ∂ h x (B n dy + B n ry ) -∂ h y B n x E n z + Jn+ 1 2 = E n z + ∆t curl h B n+ 1 2 d E n z + Jn+ 1 2 = E n z + ∆t curl h B
The following result shows that the above scheme for the coupled Maxwellsplit system is consistent with the conventional Yee scheme for the Maxwell's equations (33) (i.e. gives the same result) in Ω ∪ Ω s . Lemma 13. Consider a discrete system [START_REF] Bidégaray-Fesquet | Stability of FD-TD schemes for Maxwell-Debye and Maxwell-Lorentz Equations[END_REF], defined in Ω ∪ Ω s with the initial and right-hand side data supported inside Ω, and with the boundary condition E × n = 0, B • n = 0 on the boundary of Ω ∪ Ω s . Let the solution of this system at a time step n be E n m , B

n+ 1 2 m , J n m . Let E n ⊥ , B n+ 1 2 ⊥ , E n z , B n+ 1 2 z , J n , B n+ 1 2 , E n r , E n d , B n+ 1 2 r , B n+ 1 2 d
be the corresponding solution at a time step n to the coupled system, namely, [START_REF] Bidégaray-Fesquet | Stability of FD-TD schemes for Maxwell-Debye and Maxwell-Lorentz Equations[END_REF] on Ω, with the same initial and right-hand side data, and [START_REF] Collino | Optimizing the perfectly matched layer[END_REF]35) on Ω s , with vanishing initial and right-hand side data, coupled through the interface with the help of the coupling scheme (37, 38), with the boundary conditions (36). We define the extension of E n ⊥ and B

n+ 1 2 ⊥
to the domain Ω s as

E n ⊥ | Ω s := E n r + E n d , B n ⊥ | Ω s := B n r + B n d .
Then, for all n ≥ 0,

E n m -E n = 0, B n+ 1 2 m -B n+ 1 2 = 0, J n m -J n = 0.
Proof. The proof of this result can be done similarly to the proof of Theorem 11. More precisely, it suffices to show that E r , B r are discrete 2D-curl-free and E d , B d are discrete 2D-divergence-free, see Remark 5, which is a consequence of the fact that the Yee scheme preserves the discrete divergence of a vector.

An immediate consequence of the above result is the conservation of a discrete energy associated with the Yee scheme applied to the Maxwell dispersive system, where the components of the electric and magnetic field in Ω s are defined like in the above lemma. Remark 6. For the separate components of the field, namely B r , B d , E r , E d , the proof of the stability of the formulation is not as trivial. One could formulate the stability in terms of the components of the 2D Helmholtz decomposition of the fields E ⊥ and B ⊥ , by performing the discrete decomposition in the whole space and demonstrating that these components satisfy (TE) and (TM). Then the corresponding energies will again be inherited from the Yee scheme.

Remark 7. Let us remark that the CFL condition for the discretization [START_REF] Bidégaray-Fesquet | Stability of FD-TD schemes for Maxwell-Debye and Maxwell-Lorentz Equations[END_REF] (and thus for [START_REF] Collino | Optimizing the perfectly matched layer[END_REF] and (35)) is not worse than the one in the vacuum, cf. [START_REF] Bidégaray-Fesquet | Stability of FD-TD schemes for Maxwell-Debye and Maxwell-Lorentz Equations[END_REF].

The above result shows that the suggested coupling is exact on the discrete level. In practice, however, there are some roundoff errors, which do not affect the stability of the algorithm though, see Section 6.1.

Numerical Experiments

To demonstrate the stability of the new perfectly matched layer technique, as well as of the newly introduced coupling procedure, we perform several numerical experiments. We split the system for |x| > L x or |y| > L y , where L x = L y = L ⊥ . In all the experiments, we choose the spatial discretization ∆x = ∆y = ∆z = h. The time step in most of the experiments is chosen as ∆t = h 2 , see Remark 7. The solution is computed on the interval [0, T ].

In the experiments we use the quadratic profile of the functions σ x,y,z (.), namely, for the PML in the direction x for x > L s x = L s ,

σ x (x) = σ x0 (x -L s x ) 2 , x ≥ L s x , 0, x < L s x .
, σ x0 > 0.

In all the experiments we choose σ x0 = σ y0 = σ z0 = σ 0 . One could alternatively use piecewise-constant function σ x,y,z (.); for a correct treatment of such functions in the context of the Yee scheme see e.g. [START_REF] Collino | Optimizing the perfectly matched layer[END_REF]. In all the directions we choose the same width of the PML. The width of the PML layer we denote by W := L pml β -L s β = L pml -L z , β ∈ {x, y}, and the width of the layer where the split system is solved but no PML is applied we denote by W s := L s x -L x = L s y -L y . In all the experiments we use zero initial conditions E 0 = B 1 2 = 0, J 0 = 0, and non-zero right hand side data f (x, y, z, t) s.t. f (x, y, z, 0) ≈ 0.

The Validity of the Splitting (Illustration of the Statement of Lemma 13)

As we have shown in Lemma 13, we should expect the coupling procedure of Section 5.2 to be exact on the discrete level. To support this claim, in this section we compare the solution to the problem (M) with zero Dirichlet boundary conditions computed with the help of [START_REF] Bidégaray-Fesquet | Stability of FD-TD schemes for Maxwell-Debye and Maxwell-Lorentz Equations[END_REF] to the solution of the coupled Maxwell-split system computed with the help of the procedure of Section 5.2. This can be viewed as a demonstration of Lemma 13. We choose non-zero right hand side data in the equations for B x , E x and E z of the original Maxwell system (M), supported inside the domain Ω (the exact form is not important for the discussion). The rest of the parameters are given in Table 6 a continuous energy functional (notice that it does not coincide with a discrete energy of the Yee scheme), which we will use to measure errors and quantify the stability:

E n (E, J, B) = E n 2 + ω -2 p J n 2 + B n 2 . ( 39
)
Here . is an L 2 -discrete norm, defined as:

u ijk 2 = h 3 i j k u 2 ijk . (40)
Given E m , B m , J m the solution of the discretized Maxwell system [START_REF] Bidégaray-Fesquet | Stability of FD-TD schemes for Maxwell-Debye and Maxwell-Lorentz Equations[END_REF], and E, B, J the solution to the coupled system computed as in Lemma 13, we define the absolute and relative errors at a time step n e a n = E

1 2 n (E -E m , J -J m , B -B m ), e r n = e a n E -1 2 n (E m , J m , B m ). (41)
In Figure 6.1 we plot the relative error e r n depending on a time t = n∆t. Notice that the error is measured in the whole computational domain Ω ∪ Ω s . We can see that e r n grows linearly on time, however, very slowly; moreover, it remains of the order of a roundoff error. r on time for the solution computed in the experiment without the PML with the parameters in Table 6.1.

Stability of the New PML technique

To demonstrate the stability of the new PML technique, we perform a numerical experiment with parameters given in Table 6.2. The non-zero components of the right-hand side data are, with r = x 2 + y 2 + z 2 , F Ex = 50e -500(t-0.2) 2 -300r 2 , F Ez = 20(t -0.1)e -10 4 (t-0.1) 2 -400r 2 .

In Figure 6.2 (resp. 6.3) we show solutions E x (resp. E z ) obtained with the help of the new PML technique, and the blow-up of solutions obtained with the old (Bérenger's) PML technique. In order to quantify the stability, we measure the energy (39) inside the whole computational domain Ω ∪ Ω s ∪ Ω pml . The results in Figure 6.4 show that the old PML exhibits the instability, while the new perfectly matched layer remains stable. We show the approximation of a continuous energy E (39) for stable and unstable (old) PMLs, computed for the experiment of Section 6.2, inside the full domain Ω ∪ Ω s ∪Ω pml .

ω p L ⊥ L z T h ∆t W s W σ 0 40 0.

The Validity of the Solution

In the following experiments, we compare the solution obtained with the help of the PML in a smaller ('physical') domain to the solution to the Maxwell equations in plasmas with zero Dirichlet boundary conditions computed on a larger domain Ω D , whose size is chosen so that the reflections from the exterior boundary do not reach the physical domain during the time of the simulation. We choose a non-vanishing right hand side for the equations for E x , E z (the actual value is not of importance for the present discussion), whose L ∞ -norm is less than 10 -6 outside of the interval t ∈ [0.01, 0.38], and study how the error depends on the parameters of the PML. The choice of the PML absorption parameter, as well as the PML width, are illustrated in Figure 6.3. To quantify the error, let us notice that the energy (39) measured in the large domain Ω D stabilizes after some moment t = t s ; we will denote this stabilized value E n (E m , J m , B m ) by E * . For large t the energy inside the domain almost vanishes, while the PML error, strictly speaking, vanishes with a different rate (which depends on the rate of reflection/absorption of the PML). Thus, it makes sense to measure the scaled (rather than relative) error of the ith experiment

C = Ω ∪ Ω s No. ω p L x,y L z T h ∆t W s W σ 0 1 60 0. 9
ẽr i (n∆t) = E 1 2 n ( (E -E m )| C , (J -J m )| C , (B -B m )| C ) E -1 2 * . (42)
As before, to compute the actual electromagnetic field, inside the region where we solve the split system, we employ the identities layer width, which is larger in the third experiment. We can clearly see that the increase of the width allows to decrease the error for t ≤ 8, however, for t = 8 the error magnitude becomes almost the same. This may be connected to the fact that at this moment of time the dominant error is the discrete PML reflection error (which comes from the fact that on the discrete level the absorption parameter σ(x) is no longer continuous), while the continuous PML error becomes negligible. For the discussion of the discrete PML error see [START_REF] Collino | Optimizing the perfectly matched layer[END_REF], while its continuous part had been analyzed in [START_REF] Diaz | A time domain analysis of PML models in acoustics[END_REF], see also [START_REF] Joly | An elementary introduction to the construction and the analysis of perfectly matched layers for time domain wave propagation[END_REF] for a summary of those. As for the left plot of Figure 6.3, we observe that the increase of the absolute value of the absorption parameter allows to decrease the error, however, up to a certain moment. This is likely to be related to the fact that at some moment of time the discrete PML reflection error dominates, and is only very slightly absorbed by the PML (see the time t ≥ 5). This is not surprising if one compares the profiles of σ(x) for different experiments.

E ⊥ = E c + E g and B ⊥ = B c + B g .

Conclusions

In this article we have presented a method to construct stable perfectly matched layers for a cold plasma in an infinitely large background magnetic field. The difficulty of the application of the PMLs to such a model is two-fold.

First, it is the presence of the dispersion and anisotropy, which leads to the instability of the traditional, Bérenger's PML, both in 2D and 3D. In this work it has been overcome by the use of the method of [START_REF] Bécache | On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials[END_REF]. It is demonstrated that in the free space and for constant absorption parameters, the energy of the resulting PML system does not grow; the stability of the PML is confirmed by the numerical experiments.

Second, it is the presence of the backward and forward propagating modes at the same time, a difficulty inherent to the 3D plasma model. We have suggested to construct an equivalent system, which essentially consists of two systems, one 2D-plasma like, and another similar to the Maxwell equations in vacuum. Then different types of perfectly matched layer techniques were applied to each of these systems. We have shown that the constructed PMLs are uniformly stable.

The coupling between the Maxwell system, which we solve in the physical domain, and the split system, which we solve in the PML domain, is done through an artificial layer on a discrete level. In this work we have suggested a procedure that allows to couple these systems exactly in the discrete setting. Our numerical experiments support the claims of the stability of the PMLs.

While the procedure of splitting described in this work has a potential for the extension to a class of slightly more general models, it fails even in a simple physical case of a one-species cold plasma model under arbitrary magnetic field. In this case no local splitting is possible, and such models require a separate treatment, probably with techniques different from the PML. This is a subject of a future research.

Appendix A. Derivation of the System [START_REF] Abarbanel | Well-posed perfectly matched layers for advective acoustics[END_REF] Let us now discuss how the system (15) had been derived. For convenience, we rewrite it here:

∂ t B x + σ z B x + B *
x + ∂ y (E z + σ z j p ) -∂ z E y = 0, (A.1a) The first equation in the above system was obtained by performing the PML change of variables in z and y in the first equation of (TE) x :

∂ t E y + σ z E y -∂ z B x =
iω Bx + 1 + σ y ψ(ω) iω -1

∂ y Êz -1 + σ z iω -1
∂ z Êy = 0.

Thanks to the fact that σ z does not depend on y, the above can be rewritten:

iω Bx + σ z Bx + 1 + σ y ψ(ω) iω In the time-domain this is indeed (A.1e, A.1f).

As for the rest of the equations, the expression (A.1b) is a result of a trivial application of the 1D PML in the direction z. The equation (A.1c) is obtained from the corresponding expression of (TE) x , which can be rewritten by substituting  = (iω) -1 Êz and ψ(ω) = (1 + (iω) -2 ω 2 p ) -1 :

iω + ω 2 p (iω) -1 Êz = 1 + σ y ψ(ω) iω In the time domain this is indeed (A.1c, A.1d).

Appendix B. Derivation of the PMLs for Systems (TE) and (TM)

First of all, let us concentrate on a slightly more complicated case of the PML system (TM), which reads ∂ z Êy = 0.

∂ t B d + σ z B d + B * d + curl ⊥ (E z + σ z j p ) -α ∂ z E r =
We multiply it by 1 + σz iω , which, together with p = (iω) -1 Êz gives iω Bdx + 1 + σ y ψ(ω) iω In the time domain this indeed corresponds to (B.1f, B.1g).

The derivation of the PML system for (TE) follows the same path, we the exception that the change of the variables is as per Bérenger, i.e. as if [START_REF] Kachanovska | Stable perfectly matched layers for a class of anisotropic dispersive models. Part II: Energy Estimates and Stability in Corners[END_REF] one took ψ(ω) = 1. Notice that (A. Let us now show an example of how to derive the first equation in

∂ t E d + σ z E d + E * d -curl ⊥ (B z + σ z K z ) + α ∂ z B r =

Figure 3 . 1 :

 31 Figure 3.1: The snapshots of the field Ez computed with the help of (TE)x combined with the conventional, Bérenger's change of variables, at t = 10, 20, 25, 40.

Figure 3 . 2 :

 32 Figure 3.2: In the left figure we schematically depict level curves of the modes ±ω -(0, ky, kz) ≤ ωp, and in the right the level curves of the modes ±ω -(0, ky, kz) ≥ ωp.

  .1, see Section 3.1, however, with the help of the new perfectly matched layer. The results are demonstrated in Figure 3.3.

Figure 3 . 3 :

 33 Figure 3.3: The snapshots of the field Ez computed with the help of (TE)x with the help of the new PML model (15) at t = 10, 20, 25, 40, with the same data as in Figure 3.1.

Lemma 9 .

 9 The solutions of the dispersion relation of the system (4) satisfy: 1. for |ω| < |ω p | there are four modes present: ω 2 = ω 2 -, and ω 2 = ω 2 i . The modes ±ω -(k) are backward in the directions e x , e y for all k ∈ R d and are forward in e z . The modes ω = 0 and ±ω i (k) are forward e x , e y , e z . 2. for |ω| ≥ |ω p | there are four (not counting multiplicity) modes present: ω 2 = ω 2 + and ω 2 = ω 2 i . All these modes are forward in e x , e y , e z . An illustration to this statement is shown in Figure 4.1, where, due to rotational invariance of ω(k) in coordinates k x , k y , we show the curves ω(k) = const in coordinates (|k ⊥ |, k z ), where k ⊥ = (k y , k z

Figure 4 . 1 :

 41 Figure 4.1: The level sets of ω(k). The dashed circle corresponds to the mode ω i (k).

Figure 4 . 2 :

 42 Figure 4.2: In the left figure we depict the domains Ω and Ω s . In the right figure we show the cut y = const of the domains Ω, Ω s and Ω pml .

Figure 5 . 1 :

 51 Figure 5.1: The discretization of the fields in the Maxwell and split systems.

Figure 5 . 2 :

 52 Figure 5.2: Coupling procedure for Bz. In blue we denote the boundary Γs between the regions where the split and non-split systems are solved.

Figure 5 . 2 )

 52 , inside which we perform coupling. Let us now explain how we compute B z n+ 1 2

n+ 1 2 ⊥Figure 5 . 3 :

 253 Figure 5.3: Coupling procedure for Ez. In blue we denote the boundary Γs between the regions where the split and non-split systems are solved.

Figure 6 . 1 :

 61 Figure 6.1: The dependence of the error e n r on time for the solution computed in the experiment without the PML with the parameters in Table 6.1.

  9 1.01 25 0.01 0.005 0.11 0.19 8000

Figure 6 . 2 :Figure 6 . 3 : 4 tE

 62634 Figure 6.2: The field Ex(x, y, 0, t) computed for the experiment of Section 6.2. By a dashed line we denote a PML boundary. All figures use the same color scale.

Figure 6 . 5 :

 65 Figure 6.5: Absorption parameters for different experiments.

  The results of the experiment are shown in Figure6.3. We can see that varying the absorption parameter σ(x) and the width of the PML we can control the error of the PML solution. Let us first address the plot in the right in Figure6.3. The difference between the PML of the experiment 2 and 3

1 ∂

 1 expression follows from ψ(ω) = (1 -ω 2 p ω -2 ) -1 . One chooses to introduce the auxiliary unknown B * by comparing (A.2) to (A.1a), and with the use p = (iω) y Êz + σ z p -∂ y ( Êz + σ z p ) = -σ y iω + ω 2 p (iω) -1 + σ y -1 ∂ y Êz + σ z p ,with the use of (A.3). Alternatively, iω + ω 2 p (iω) -1 + σ y B * x = ∂ y Êz + σ z p .

- 1 ∂ 1 - 1 ∂

 111 y Bx = 1 + σ y iω + ω 2 p (iω) -y Bx , or, equivalently, iω + ω 2 p (iω) -1 + σ y Êz = ∂ y Bx .

- 1 ∂ 1 ∂ 1 ∂ y Bdx + 1 + σ x ψ iω - 1 ∂ 1 - 1 1 - 1

 111111111 y Êz + σ z p -∂ z Êy = 0.Comparing the above to (B.1a), we get, with the use of (A.y Êz + σ z p -∂ y Êz + σ z p = -σ y iω + ω 2 p (iω) -1 + σ y ∂ y Êz + σ z p ,or iω + ω 2 p (iω) -1 + σ y B * dx + σ y ∂ y Êz + σ z p = 0. One sees that in the time domain this coincides with (B.1d). The equation (B.1a) for Bdy , and thus, the identities (B.1e) are obtained in the same manner.The equation (B.1b) follows trivially by applying the PML in the direction z and multiplying both sides by (1 + (iω) -1 σ z ). Finally, the equation (B.1c) can be obtained as follows. Applying the PML change of variables[START_REF] Kachanovska | Stable perfectly matched layers for a class of anisotropic dispersive models. Part II: Energy Estimates and Stability in Corners[END_REF] to the third equation of (TM) in the frequency domain, and using  = ω -2 p Ĵ, we obtainiω Êz -1 + σ y ψ iω x Bdy + ω 2 p  = 0.Like before, comparing the above to (B.1c), we introduceÊ * zx = -1 + σ y ψ iω -∂ y Bdx (A.3) = σ y iω + ω 2 p (iω) -1 + σ y ∂ y Bdx , ∂ x Bdy (A.3) = -σ x iω + ω 2 p (iω) -1 + σ x ∂ x Bdy .

∂

  t K z = B z . (B.3h)Consider the equation of (TE) for E dx with the Bérenger's change of variables:iω Êdx -(1 + σ y (iω) -1 )∂ y Bz + (1 + σ z (iω) -1 )∂ z By = 0.We multiply both sides by (1 + σ z (iω) -1 ), and use the fact that ∂ y σ z ≡ 0:iω Êdx + σ z Êdx -(1 + σ y (iω) -1 ) -1 ∂ y Bz + σ z (iω) -1Bz + ∂ z Bry = 0.Comparing the above to (B.1a) and introducing Kz = (iω) -1 Bz , see (B.3h), we obtain the identityÊ * dx = -(1 + σ y (iω) -1 ) -1 -1 ∂ y Bz + σ z Kz (B.2) = σ y iω + σ y ∂ y Bz + σ z Kz .In the time domain, this is indeed (B.3d). The rest of the equations in (B.1) are derived similarly to (B.1).

  coupled with the equations for the PML unknowns E *

d = (E * dx , E * dy ), B * zx , B * zy and K z , which read

Table 6 .

 6 1: Parameters for the experiment of Section 6.1 without the PML.

	.1. Let us define

Table 6 . 2 :

 62 Parameters for the experiment of Section 6.2.

Table 6 .

 6 3: Parameters for the experiments in Section 6.3.

	300				σ(x), Exp. 1
					σ(x), Exp. 2
	200				σ(x), Exp. 3
	100
	0	0	0.1	0.2	0.3	0.4
				x

  ∂ t E z + σ y E z + ω 2 p j p + ∂ y B x = 0, (A.1c) ∂ t j p = E z , (A.1d) ∂ t B * x + ω 2 p K * x + σ y B * x + σ y ∂ y (E z + σ z j p ) = 0, (A.1e)

	0,	(A.1b)
	∂ t K * x = B * x .	(A.1f)

  0, (B.1a)∂ t E r + σ z E r + α∂ z B d = 0, (B.1b) ∂ t E z + E * zx + E * zy -curl ⊥ B d + ω 2 p j p = 0, ∂ t j p = E z , (B.1c) ∂ t B * dx + ω 2 p K * dx + σ y B dx + σ y ∂ y (E z + σ z j p ) = 0, ∂ t K * dx = B * dx , (B.1d) ∂ t B * dy + ω 2 p K * dy + σ x B dy -σ x ∂ x (E z + σ z j p ) = 0, ∂ t K * dy = B * dy , (B.1e) ∂ t E * zx + ω 2 p j * px + σ y E * zx -σ y ∂ y B dx = 0, ∂ t j * px = E * zx , (B.1f) ∂ t E * zy + ω 2 p j * py + σ x E * zy + σ x ∂ x B dy = 0, ∂ t j * py = E * zy . (B.1g)Let us show how the equation (B.1a) had been derived. First of all, applying the PML change of variables[START_REF] Kachanovska | Stable perfectly matched layers for a class of anisotropic dispersive models. Part II: Energy Estimates and Stability in Corners[END_REF] to the first equation of (TM) in the Fourier domain, we obtain, for Bdx :

	iω Bdx + 1 +	σ y ψ(ω) iω	-1	∂ y Êz -1 +	iω σ z	-1

  ∂ t B r + σ z B r -α∂ z E d = 0, (B.3b) ∂ t B z + B * zx + B * zy + curl ⊥ E d = 0, (B.3c) ∂ t E * dx + σ y E * dx -σ y ∂ y (B z + σ z K z ) = 0, (B.3d) ∂ t E * dy + σ x E * dy + σ x ∂ x (B z + σ z K z ) = 0, (B.3e) ∂ t B * zx + σ y B * zx + σ y ∂ y E dx = 0, (B.3f) ∂ t B * zy + σ x B * zy -σ x ∂ x E dy = 0, (B.3g)

	0,	(B.3a)

The work of the third author was supported by the ANR CHROME (ANR-