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

Patrick Joly
email: patrick.joly@inria.fr

Valentin Vinoles
email: cvalentin.vinoles@epfl.ch

On the analysis of perfectly matched layers for a class of dispersive media and application to negative index metamaterials

This work deals with Perfectly Matched Layers (PMLs) in the context of dispersive media, and in particular for Negative Index Metamaterials (NIMs). We first present some properties of dispersive isotropic Maxwell equations that include NIMs. We propose and analyse the stability of very general PMLs for a large class of dispersive systems using a new change of variable. We give necessary criteria for the stability of such models that show in particular that the classical PMLs applied to NIMs are unstable and we confirm this numerically. For dispersive isotropic Maxwell equations, this analysis is completed by giving necessary and sufficient conditions of stability. Finally, we propose new PMLs that satisfy these criteria and demonstrate numerically their efficiency.

Introduction

One of the hardest difficulties to simulate wave propagation in unbounded domains is to construct artificial boundary conditions which absorb the outgoing waves without reflecting them into the computational domain. A widely used technique to do so is the so-called Perfectly Matched Layer (PML) [START_REF] Joly | An elementary introduction to the construction and the analysis of perfectly matched layers for time domain wave propagation[END_REF] first proposed by Bérenger [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] for the 3D Maxwell equations. The method consists in surrounding the computational domain by an absorbing layer (the PML) which generates no reflection inside it. Even if they are very effective in many cases, PMLs can exhibit instabilities for some problems because of backward waves, which are waves whose phase and group velocities point in "opposite" directions with respect to the interface. For non-dispersive media, the first two authors together with S. Fauqueux in [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF] established a necessary (but not sufficient) criterion of stability: classical PMLs are stable only if there is no backward wave. This is always the case for isotropic media but not for all anisotropic ones. For some of them (e.g. anisotropic acoustic equations [START_REF] Demaldent | Perfectly matched transmission problem with absorbing layers: Application to anisotropic acoustics in convex polygonal domains[END_REF], aeroacoustics [START_REF] Diaz | A time domain analysis of pml models in acoustics[END_REF]), new stable PLMs have been successfully proposed but, to our knowledge, the case of anisotropic elastic waves remains open [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF].

The case of dispersive systems has been much less studied. A very important example of such media are metamaterials, i.e. artificial composite materials having extraordinary electromagnetic properties. In particular, Negative Index Metamaterials (NIMs), also called left-handed media or double negative metamaterials, have negative permittivity and permeability at some frequencies due to microscopic resonating structures [START_REF] Smith | Metamaterials and negative refractive index[END_REF]. Since the 1990s, NIMs are the subject of active researches due to their promising applications [START_REF] Tie | Metamaterials: theory, design, and applications[END_REF]: superlens, cloaking, improved antenna, etc. The fact is that, even in isotropic media, NIMs naturally support backward waves (at least in some range of frequencies), which leads to anticipate difficulties with PMLs, by analogy with what occurs with non dispersive media.

PMLs in NIMs have already been studied by the physicists community [START_REF] Correia | 3d-fdtd-pml analysis of left-handed metamaterials[END_REF][START_REF] Steven | Perfectly matched layer behavior in negative refractive index materials[END_REF][START_REF] Dong | Perfectly matched layer-absorbing boundary condition for left-handed materials[END_REF][START_REF] Loh | Fundamental relation between phase and group velocity, and application to the failure of perfectly matched layers in backward-wave structures[END_REF][START_REF] Shi | Perfectly matched layer absorbing boundary condition for truncating the boundary of the left-handed medium[END_REF]. To our knowledge, in [START_REF] Steven | Perfectly matched layer behavior in negative refractive index materials[END_REF], Cummer was the first to notice that classical PMLs fail in NIMs and gave a physical explanation related to the presence of backward waves. He proposed stable PMLs for a particular case of Drude materials corresponding to ω e = ω m (see (2.12) in Section 2.1.3 for more details). This work was extended in [START_REF] Dong | Perfectly matched layer-absorbing boundary condition for left-handed materials[END_REF][START_REF] Shi | Perfectly matched layer absorbing boundary condition for truncating the boundary of the left-handed medium[END_REF] to the general case ω e = ω m . In all these works the instability of PMLs was observed through numerical and/or explained by arguments coming from the physics.

The goal of this paper is to bring a mathematical eye to the question of the stability of PMLs for NIMs and also to generalize the construction of stable PMLs for a much larger class of dispersive electromagnetic media than the Drude ones. Let us mention a first paper [START_REF] Bécache | Perfectly matched layers in negative index metamaterials and plasmas[END_REF] by the authors in which we have already claimed some of the results of the present article, in particular how to construct stable PMLs for the Drude model, but without any proof (we also extended our method to a simple plasma model). We also mention that this paper is an advanced version of the preliminary study presented in the last author's PhD thesis [START_REF] Vinoles | Problèmes d'interface en présence de métamatériaux: modélisation, analyse et simulations[END_REF].

The outline of the article is as follows. In Section 2, we present the class of mathematical models that we consider in this paper for isotropic dispersive electromagnetic media and give their main mathematical properties (Section 2.1) as well as an analysis of dispersion phenomena (Section 2.2), introducing in particular the definition of backward modes and negative index.

The main purpose of the rest of the paper is to propose, in a constructive manner, a generalized PML that would be stable even in the presence of backward waves, which are responsible of the instability of the classical Bérenger's PMLs. The key point is to introduce a damping parameter which depends on the frequency ω through a function χ(ω) (a rational fraction in ω) which will be cleverly chosen in order to deal with the dispersive properties of the medium and more specifically with the presence of backward waves.

Sections 3 and 4 contain the main results of this paper. In Section 3 we consider the question of PMLs for non-dissipative dispersive first order hyperbolic systems. These are presented in Section 3.1 and generalize the models of Section 2. In Section 3.2, we propose generalized PMLs for such systems and initiate their stability analysis (the main issue of this paper) in Section 3.3. In particular, we derive two necessary stability conditions in Sections 3.3.4 and 3.3.5. In Section 4, we come back to dispersive isotropic Maxwell equations and complete, in this case, the stability analysis initiated in Section 3.3. We obtain a necessary and sufficient stability condition (Theorems 4.5 and 4.10) for our generalized PMLs. Finally, in Section 5, we construct stable PMLs for Drude materials (Section 5.1, with numerical illustration), for generalized Lorentz materials (Section 5.2) and propose a procedure to construct stable PMLs for general models in Section 5.3.

Isotropic dispersive media

In this section, we consider a general class of mathematical models for the propagation of electromagnetic waves in dispersive, isotropic non-dissipative 2D media. The restriction to the 2D case is really non essential and adopted only for simplicity.

The mathematical models

The Maxwell equations

We consider the two dimensional Maxwell equations in the Transverse Electric (TE) mode in the free space [START_REF] John | Classical electrodynamics[END_REF]     

∂ t D x = ∂ y H, ∂ t D y = -∂ x H, ∂ t B = ∂ y E x -∂ x E y , ⇐⇒ ∂ t D = curl H, ∂ t B = -curl E, (2.1)
where E := (E x , E y ) T and D := (D x , D y ) T are the electric field and the electric induction, H and B the magnetic field and the magnetic induction. Here we use the differential operators curl u := ∂ x u y -∂ y u x and curl u := (∂ y u, -∂ x u) T . To close this system, one must add some constitutive relations which link E and D on the one hand, H and B on the other hand. Classically, constitutive relations for dispersive models are described in the frequency domain and relate the time Fourier transforms of the fields (we adopt for the Fourier transform the convention that transforms the time derivative ∂ t into the multiplication by iω):

D(ω, •) = ε(ω, •) E(ω, •) and B(ω, •) = µ(ω, •) H(ω, •), (2.2)
where ε(ω, x, y) and µ(ω, x, y) are respectively the permittivity and the permeability at the point (x, y) and frequency ω. In the rest of the paper, we shall restrict our presentation to homogeneous media, that is to say to the case where ε(ω, x, y) = ε(ω) and µ(ω, x, y) = µ(ω). This is not restrictive as far as the question of PMLs is concerned, PMLs are used in homogeneous regions. However, the reader will observe that a large part of this section still applies to heterogeneous media. In the frequency domain, (2.1) becomes after eliminating D and B thanks to (2.2):

iω ε(ω) E = curl H, iω µ(ω) H = -curl E. (2.3)
When ε and µ do not depend on the frequency ω, the medium is a standard dielectric one, otherwise it is a so-called dispersive medium.

Admissible forms of the permittivity and the permeability

Even though the content of this section has to be considered as rather formal from the mathematical point of view, it is worth mentioning that mathematical properties of ε(ω) and µ(ω) are usually required in order to ensure the well-posedness and physical soundness of the evolution problem corresponding to (2.3). In particular, in order to preserve the causality of the constitutive laws (2.2) and the fact that the electric and magnetic fields in time domain are real valued, one must impose that ε(ω) and µ(ω) are Fourier transforms of real and causal distributions.

(2.4)

In particular, ε(ω) and µ(ω) can be extended to the complex plane as analytic functions in the lower half-plane C -:= {z ∈ C : Im z ≤ 0}. In the following, we shall also assume that for ω ∈ R, ε(ω) and µ(ω) are real valued, (2.5) which in practice amounts to consider non dissipative media (we shall give a more precise definition in the next section). In many cases, dispersive media are also (slightly) dissipative but we are convinced that the question of PMLs is more difficult for non dissipative media. Finally we shall also assume that both ε(ω) and µ(ω) have a finite limit at high frequency: there exist two real numbers ε ∞ > 0 and µ ∞ > 0 such that

ε(ω) -→ ω→+∞ ε ∞ and µ(ω) -→ ω→+∞ µ ∞ .
(2.6)

In other words, such media behave as a standard dielectric medium at high frequency and, from the mathematical point of view, the model in time that corresponds to (2.3) can be seen as a lower order perturbation of classical Maxwell equations associated to (ε ∞ , µ ∞ ).

For the simplicity of our analysis, but also because it already covers a large class of models of the literature, we shall restrict ourselves to local dispersive media in the sense that ε(ω) and µ(ω) are rational fractions. Doing so, if one takes into account (2.4), (2.5) and (2.6), one must impose that:

ε(ω) and µ(ω) are real, even rational fractions of ω with a constant entire part and real poles.

(2.7)

The justification of the terminology "local" is due to the fact that, when coming back to the time domain from (2.2), the relationship between E and D (resp. H and B) can be expressed in terms of ordinary differential equations, modulo the introduction of auxiliary unknowns.

In what follows, we shall restrict ourselves to dispersive media satisfying (2.7). In these cases, the evolution problem corresponding to (2.3) can be rewritten as the standard Maxwell equations coupled to ordinary differential equations. Indeed, the relationships (2.2) actually have to be understood in the sense

Q e (iω) D(ω, •) = P e (iω) E(ω, •) and Q m (iω) B(ω, •) = P m (iω) H(ω, •), (2.8)
where ε(ω) = P e (iω)/Q e (iω) and µ(ω) = P m (iω)/Q m (iω), with P e and Q e real polynomials of degree 2M e , P m and Q m real polynomials of degree 2M m such that P e P m and Q e Q m are coprime.

(2.9)

We also introduce the following notations:

P em := P e P m and Q em := Q e Q m .
(2.10)

In the time domain, (2.8) becomes a system of ordinary differential equations

Q e (∂ t ) D = P e (∂ t ) E and Q m (∂ t ) B = P m (∂ t ) H. (2.11)
Remark 2.1. Condition (2.9) is not restrictive since it always can be satisfied up to a change of unknowns.

We give below two examples of such models.

Figure 1 -Plot of ε(ω) and µ(ω) defined by the Drude model (2.12)

The Drude model

One of the simplest local dispersive medium is the Drude one that corresponds to

ε(ω) = ε 0 1 - ω 2 e ω 2 and µ(ω) = µ 0 1 - ω 2 m ω 2 , (2.12)
where ε 0 > 0 and µ 0 > 0 are the permittivity and the permeability of the vacuum, and ω e > 0 and ω m > 0 are the electric and the magnetic plasma frequencies. Both functions ε(ω) and µ(ω) have a changing sign as illustrated by Figure 1. The Drude model is widely used to describe the permittivity in metals [16, §7.5]. For the permeability, the Drude model can be thought as a simplification of more complicated models which describe artificial magnetism in resonating structures [START_REF] Brian | Negative refraction makes a perfect lens[END_REF][START_REF] Richard | Wave propagation in media having negative permittivity and permeability[END_REF]. A derivation of this model can be found in [START_REF] Stephen | Photonic band-gap effects and magnetic activity in dielectric composites[END_REF] (see also [START_REF] Bouchitté | Homogenization of the 3d maxwell system near resonances and artificial magnetism[END_REF][START_REF] Bouchitté | Homogenization of maxwell's equations in a split ring geometry[END_REF] for mathematical justifications based on high-contrast homogenization). Let us emphasize that, despite its simplicity, the Drude model (2.12) already contains all the difficulties to construct stable PMLs. That is why we consider this model as our privileged toy problem. To write the Drude model in the time domain, we first substitute (2.12) into (2.3):

       ε 0 iω E + ω 2 e iω E = curl H, µ 0 iω H + ω 2 m iω H = -curl E.
(2.13)

We now introduce auxiliary fields J := (J x , J y ) T and K (called the induced electric and magnetic currents) which are the primitives in time of the fields E and H :

iω J = E and iω K = H. (2.14)
We obtain the TE first-order system for the Drude model in the time domain

         ∂ t E + ω 2 e J = ε -1 0 curl H, ∂ t J -E = 0, ∂ t H + ω 2 m K = -µ -1 0 curl E, ∂ t K -H = 0.
(2.15)

Generalized Lorentz models

A generalization of the Drude model (2.12) is the generalized Lorentz materials corresponding to :

ε(ω) = ε 0   1 - Me j=1 ω 2 e,j ω 2 -Ω 2 e,j   and µ(ω) = µ 0 1 - Mm k=1 ω 2 m,k ω 2 -Ω 2 m,k .
(2.16)

The reader will easily check that the Drude model (2.12) is a particular case of (2.16) with M e = M m = 1 and Ω e,1 = Ω m,1 = 0. The time domain model corresponding to (2.16) can be written (using auxiliary unknowns as we did for the Drude model, we omit the details) as a first order system of dimension m = 2(M m + M e + 1) + 2M e + 1:

                       ∂ t E + Me j=1 ω 2 e,j J j = ε -1 0 curl H, ∂ t J j + Ω 2 e, P j -E = 0, ∂ t P j -J j = 0, (1 ≤ j ≤ M e ), ∂ t H + Mm k=1 ω 2 m,k K k = -µ -1 0 curl E, ∂ t K k + Ω 2 m,k R k -H = 0, ∂ t R k -K k = 0, (1 ≤ k ≤ M m ).
(2.17)

Note that the spatial differential operators only apply to the electric and magnetic fields E and H and not to the auxiliary fields Remark 2.2. For the generalized Lorentz models, the fields D and B can be recovered by

D = ε 0   E + Me j=1 ω 2 e,j P   and B = µ 0 H + Mm k=1 ω 2 m,k R k . (2.18)
The well-posedness of the system (2.17) is a simple application of one of the classical theories for linear evolution equations like the Hille-Yosida's theory [START_REF] Lawrence | Partial differential equations[END_REF] (the details are left to the reader). Moreover, one can easily establish the following energy identity

d dt (E e + E m ) = 0, (2.19)
where

             E e = 1 2 R 2 ε 0 |E| 2 + 1 2 Me j=1 R 2 ε 0 ω 2 e,j |J j | 2 + Ω 2 m,j |P j | 2 E m = 1 2 R 2 µ 0 |H| 2 + 1 2 Mm k=1 R 2 µ 0 ω 2 m,k |K k | 2 + Ω 2 m,k |R k | 2 .
(2.20)

Remark 2.3. It has been shown (see e.g. [START_REF] Tip | Linear dispersive dielectrics as limits of drude-lorentz systems[END_REF]) that, modulo an adequate generalization that would consist in the fact that the indexes j and k in (2.16) and (2.17) become continuous variables, the sum becoming an integral with respect to a given positive measure, the generalized Lorentz models are representative of all causal and passive electromagnetic materials. Note that such materials are in general no longer local dispersive materials in the sense of Definition (2.7).

Growing property

One has a very useful property of the generalized Lorentz models (2.16) that we call "growing property" in the following.

to conclude that d dω ω ε(ω) = ε 0 Me j=1 Ω 2 e ω 2 e,j + ω 2 ω 2 -Ω 2 e,j 2 > 0.
(2.23)

The same holds for µ(ω).

Analysis of dispersive properties

Modal analysis

The dispersion properties of (2.1)-(2.2) are classically analysed through the notion of harmonic plane waves. Definition 2.6. For a given wave vector k ∈ R 2 , a plane wave is a solution u = (E, H, B, D) T of (2.1)-(2.2) (or more rigorously of (2.1) and (2.11)) under the form u(x, t) = U e i(ωt-k•x) , (2.24) where U ∈ C 6 \ {0} is the amplitude vector and ω ∈ C the frequency.

Introducing the set F of forbidden frequencies defined as

F := ω ∈ R : ω is a pole of ε(ω) µ(ω) , (2.25)
it is easy to see using (2.9) that the existence of a plane wave implies ω ∈ O where O is the set of admissible frequencies:

O := C \ F. (2.26)
For ω ∈ O, one sees from (2.3) that (ω, k) should satisfy the dispersion relation

ω ε(ω) ω 2 ε(ω)µ(ω) -|k| 2 = 0, (2.27)
therefore either ω ε(ω) = 0, which corresponds to non-propagative (or static) waves associated to curl-free electric fields, whose 2M s + 1 solutions (for instance M s = M e for the Lorentz model (2.17)) will be denoted

ω = ω s j , -M s ≤ j ≤ M s (with ω s -j = -ω s j ), (2.28)
or (ω, k) should satisfy the reduced dispersion relation

ω 2 ε(ω)µ(ω) = |k| 2 .
(2.29) whose 2 M solutions (for instance M = M e + M m + 1 for the Lorentz model (2.17)) will be denoted

ω = ± ω j (k) = ± ω j |k| , 1 ≤ j ≤ M , (2.30)
where the functions |k| → ω j |k| are analytic functions. Remark 2.7. The zeroes of ε(ω) play here a particular role because we consider here a 2D model. In 3D, ε(ω) and µ(ω) would play a symmetric role.

Non dissipative models

Definition 2.8. The dispersive model associated to the constitutive laws (2.2) is said to be non dissipative if all solutions ω j (k) of (2.29) are real valued. By opposition, the medium is called dissipative if Im ω j (k) > 0 for some j and k.

In the sequel, we will consider only non dissipative local dispersive models. In other words, the set O of admissible frequencies (2.26) satisfies

O ⊂ R. (2.31)
For the Drude model (2.12), we can check explicitly that it is non dissipative. Indeed, injecting (2.12) in (2.29) shows that, for a given wave vector k, ω has to be a zero of the polynomial

P k (ω) := ω 4 -(ω 2 e + ω 2 m + c 2 |k| 2 ) ω 2 + ω 2 e ω 2 m .
(2.32)

where c := 1/ √ ε 0 µ 0 is the speed of light in the vacuum.

Proposition 2.9. For all wave vectors k ∈ R 2 , the polynomial P k given by (2.32) has 4 real distinct zeros ± ω j (k), j = 1, 2 depending only on |k| and given by Note that, since we consider isotropic media, phase and group velocities are collinear. In a dielectric material, in which ε and µ do not depend on ω, they are always equal but for more general models these quantities are different: this characterizes dispersive effects. Moreover, the scalar product of the two vectors can be positive (as in a non dispersive medium) or negative. This leads us to distinguish two cases: Definition 2.12. For a given wave vector k, a plane wave associated to ω(k) is called forward if

       ω 1 (k) := √ 2 2 ω 2 e + ω 2 m + c 2 |k| 2 + Σ(|k|), ω 2 (k) := √ 2 2 ω 2 e + ω 2 m + c 2 |k| 2 -Σ(|k|), ( 2
v p ω(k) • v g ω(k) > 0, (2.36)
or equivalently ω |k| ω (|k|) > 0, and backward if

v p ω(k) • v g ω(k) < 0, (2.37)
or equivalently ω |k| ω (|k|) < 0.

A mode associated to a branch ω(k) will be called forward (resp. backward) if for all k, the plane wave associated to ω(k) is forward (resp. backward).

Negative index

Thanks to Definition 2.12, we can give a precise definition of a Negative Index Metamaterial (NIM). We refer to Remark 2.18 for a justification of the denomination and further discussion. Definition 2.13. A Negative Index Metamaterial (NIM) is a local dispersive medium (in the sense of (2.7)) for which backward modes exist.

One can easily characterize NIMs only in terms of properties of the functions ε(ω) and µ(ω). More precisely let us introduce, for ω ∈ O (the set of admissible frequencies, see (2.31)) the following function:

D(ω) := 1 2ω d dω ω 2 ε(ω)µ(ω) = µ(ω) d dω ω ε(ω) + ε(ω) d dω ω µ(ω) . (2.38)
Note that, from the properties (2.7) of ε and µ, D(ω) is an even function, strictly positive as soon as ω is large enough.

Lemma 2.14. A mode is forward (resp. backward) when D ω(k) > 0 (resp. D ω(k) < 0).
Proof. The differentiation of (2.29) with respect to k leads to the fact that for any branch ω(k) one has (as soon as

D ω(k)) = 0) v g (ω(k)) = k ω |k| 1 D ω(k) , (2.39) which implies v p ω(k) • v g ω(k) = D ω(k) -1 .
(2.40)

Let us now introduce the spectrum S of the medium as the set of propagative frequencies, namely the range of admissible frequencies ω which are attained by all the branches of solution ω j (k) of (2.29):

S := M j=1 ω j (k) : k ∈ R 2 , (2.41)
and G, the set of non propagative frequencies (also called gap):

G := O \ S. (2.42)
If ω ∈ O and ε(ω) µ(ω) > 0, one can find k such that (2.29) holds. Thus ω = ω j (k) for some j. This leads to the following characterizations:

S = ω ∈ O : ε(ω)µ(ω) ≥ 0 and G = ω ∈ O : ε(ω)µ(ω) < 0 . (2.43) Note that, since ε(ω)µ(ω) is continuous in R \ F, S (resp. G) is an open (resp. closed) subset of R \ F.
N := S ∩ {ω ∈ O : D(ω) < 0} = ∅. ( 2
O = R \ {± ω e , ± ω m } and G = (-ω + , -ω -) ∪ (ω -, ω + ).
(2.45)

On the other hand, one easily computes that

D(ω) = ω 4 -ω 2 e ω 2 m c 2 ω 4 , ( 2.46)
so that N = (-ω -, ω -) = ∅. We refer to [START_REF] Vinoles | Problèmes d'interface en présence de métamatériaux: modélisation, analyse et simulations[END_REF] for numerical illustrations of the presence of forward and backward waves for the Drude model. One can wonder whether a NIM is nothing but a material for which ε(ω) and µ(ω) become both negative for some ω. In fact, as far as our definition of NIM is retained, the two things are decorrelated: a dispersive medium may be a NIM even though ε(ω) and µ(ω) remain positive. Let us take an example of such a medium (that however escapes the class of generalized Lorentz materials (2.16)):

ε(ω) = 1 - ω 2 0 ω 2 2 and µ(ω) = 1. (2.47)
The spectrum of this material is S = R while the function D is the same as for the Drude model with ω e = ω m = ω 0 . As a consequence, N = (-ω 0 , ω 0 ) = ∅ thus this material is a NIM. However, for generalized Lorentz materials, or more generally for any material that satisfies the growing property (2.21), we have the following proposition. Proposition 2.16. For any local dispersive material (in the sense of (2.7)) possessing the growing property (2.21), one has for all ω ∈ S ε(ω) D(ω) > 0 and µ(ω) D(ω) > 0.

(2.48)

Proof. Observe that ε(ω) D(ω) = [ε(ω)µ(ω] d dω (ωε(ω)) + [ε(ω) 2 ] d dω (ωµ(ω)) , (2.49)
and that all quantities inside brackets are positive when ω belongs to S. The same holds for µ.

As a consequence of this property, we have the result:

Theorem 2.17. For any local dispersive material (in the sense of (2.7)) possessing the growing property (2.21), one has

N = ω ∈ S : ε(ω) < 0 and µ(ω) < 0 . (2.50)
In other words, this medium is a NIM if, and only if, ε(ω) and µ(ω) are both negative for some ω.

Proof. From Theorem 2.15, the material is a NIM if and only if N = ∅, i.e. there exists some ω ∈ O such that ε(ω)µ(ω) > 0 and D(ω) < 0. From Proposition 2.16, a material possessing the growing property (2.21) satisfies (2.48), which implies that ε(ω) and µ(ω) become negative simultaneously.

Remark 2.18. The terminology "negative index" comes from the physics. The refraction index n of a dispersive medium is defined as n(ω) = ε(ω)µ(ω) where this quantity is to be understood as the limit of ε(ω)µ(ω) when Im ω tends to 0. One can show that n(ω) is actually negative when both ε(ω) and µ(ω) are negative. As a consequence, if one is looking at the reflection-transmission of a plane wave across a plane interface between the vacuum and a negative material, the usual Snell-Descartes law is inverted. We refer to [START_REF] Viktor | The electrodynamics of substances with simultaneously negative values epsilon and mu[END_REF][START_REF] Richard | Wave propagation in media having negative permittivity and permeability[END_REF] for more details and to [START_REF] Bécache | Perfectly matched layers in negative index metamaterials and plasmas[END_REF][START_REF] Vinoles | Problèmes d'interface en présence de métamatériaux: modélisation, analyse et simulations[END_REF] for illustrations of the effects of the negativity of refraction indexes.

Construction and analysis of PML models for a class of dispersive systems

We now want to explain why the classical PMLs are unstable for the Drude model. For linear wave propagation in non dispersive -but possibly anisotropic -media, it has been shown in [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF] that the presence of backward waves, due to anisotropy, gives rise to instabilities. It seems that it is still the case for NIMs as it was pointed out in several papers from the physicist community [START_REF] Correia | 3d-fdtd-pml analysis of left-handed metamaterials[END_REF][START_REF] Steven | Perfectly matched layer behavior in negative refractive index materials[END_REF][START_REF] Dong | Perfectly matched layer-absorbing boundary condition for left-handed materials[END_REF][START_REF] Loh | Fundamental relation between phase and group velocity, and application to the failure of perfectly matched layers in backward-wave structures[END_REF][START_REF] Shi | Perfectly matched layer absorbing boundary condition for truncating the boundary of the left-handed medium[END_REF] (see also our previous works [START_REF] Bécache | Perfectly matched layers in negative index metamaterials and plasmas[END_REF][START_REF] Vinoles | Problèmes d'interface en présence de métamatériaux: modélisation, analyse et simulations[END_REF] and Section 3.3.6).

In this section, we aim at extending the result of [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF] (valid for non dispersive systems) to a class of dispersive systems, described in Section 3.1. These models contain in particular the isotropic dispersive Maxwell models considered in Section 2. We shall consider -this is not restrictive of course -a PML in the x-direction, involving a single damping coefficient denoted σ. In fact, we propose to introduce and analyse a new class of PML models (Section 3.2) obtained by a generalization of the complex change of variable used to obtain the classical PMLs [START_REF] Joly | An elementary introduction to the construction and the analysis of perfectly matched layers for time domain wave propagation[END_REF]: in addition to the damping coefficient σ, the new change of variables involves a new ingredient, namely a real valued function of the frequency χ(ω) to be determined in practice (the choice χ(ω) = 1 giving the classical PML model). Such a generalization appears to be useful to propose an alternative to classical PMLs when these are unstable.

Finally, in Section 3.3, we provide a preliminary stability analysis for these new PML models. In particular, we establish a necessary stability condition (Proposition 3.12) which generalizes the one given in [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF]. This allows us to explain the instabilities of classical PMLs and will help in designing new stable PMLs in the next Section 4 for the isotropic dispersive Maxwell models.

We have to emphasize that our stability analysis will be performed through a modal analysis, thus limited to constant coefficients. Rigorously speaking, our analysis cannot be applied for situations of practical interest where the damping coefficient σ in the PMLs must depend on the space variable (in particular it has to be 0 in the physical domain and strictly positive inside he PML region). Nevertheless, we think that this relative weakness does not affect the interest of our analysis for the following reasons:

-The case of constant coefficients is a particular case of non constant coefficients and thus deserves to be studied anyway. Moreover, one can hardly expect that a model that leads to instability in the constant coefficients case becomes stable when used with variable coefficients. -There are numerical evidences, as we shall see in Section 5.1, that our analysis provides the good answers even in the non constant coefficients case. -Even in the case of non dispersive media, there are very few results on the stability (we speak here of stability and not only well-posedness for which there are some results, e.g. [START_REF] Halpern | The analysis of matched layers[END_REF]) of PML models with non constant σ (at the exception of [START_REF] Diaz | A time domain analysis of pml models in acoustics[END_REF] for the wave equation). Remark 3.1. Let us mention the paper [START_REF] Duru | A well-posed and discretely stable perfectly matched layer for elastic wave equations in second order formulation[END_REF] which concerns anisotropic elasticity where a specific stretching of the coordinates that are orthogonal to the absorbing boundary is performed in order to stabilize numerically the PMLs. It is not clear for us that it would work in presence of dispersion effects like backward waves. Moreover it seems that this stabilizing effect only occurs at the discrete level whereas we are interested in stability at the continuous level.

A class of first order dispersive hyperbolic systems

Consider a general 2D hyperbolic system with a zero order perturbation of the form

∂ t u + Bu = A x ∂ x u + A y ∂ y u, (x, y) ∈ R 2 , t ≥ 0, (3.1)
where the unknown u takes values in R m and A x , A y and B are real constant m × m matrices. The hyperbolic nature of the unperturbed system (B = 0) is equivalent to assuming that:

for all (k x , k y ) ∈ R 2 , k x A x + k y A y has only real eigenvalues. (3.2)
Next, we express the fact that the perturbed system is non dissipative.

Definition 3.2. The dispersive system (3.1) is said to be non dissipative if for all (k x , k y ) ∈ R 2 , k x A x + k y A y -iB has only real eigenvalues. (3.3)
Looking at |k| tends to +∞, (3.3) implies (3.2), i.e. the hyperbolicity of the unperturbed system. Definition (3.2) is justified by the Fourier analysis that permits to study the stability in terms of plane waves that are solutions of (3.1) of the form

u(x, y, t) = U e i(ωt-k•x) , ( 3.4)
where U ∈ C m is the amplitude vector, k = (k x , k y ) ∈ R 2 the wave vector and ω ∈ C the angular frequency. It is easy to check that ω and k must satisfy the so-called dispersion relation

F (ω, k x , k y ) = 0, (3.5)
where F is a polynomial in ω of order m defined by

F (ω, k x , k y ) = det(iωI + B + ik x A x + ik y A y ), (3.6)
which is nothing, up to a multiplicative coefficient, but the characteristic polynomial of

k x A x + k y A y -iB.
As a consequence, (3.5) admits m branches of solutions ω j (k), j = 1, . . . , m that is to say the eigenvalues of k x A x + k y A y -iB repeated with their multiplicity which are continuous and piecewise analytic functions of k x and k y respectively (see for instance [18, chapter II, §1]). Considering a family of plane wave solutions (a mode) associated to a smooth branch ω(k) = ω j (k), as an extension of the Definition 2.11, its group and phase velocities are defined by

v p (ω(k)) := ω k |k| k |k| and v g ω(k) := ∇ k ω(k). (3.7)
The main difference with the isotropic case is that these two vectors are not necessarily collinear.

A particular case of non dissipative systems is the symmetrisable systems which corresponds to the case where it exists a symmetric positive definite matrix M ∈ M m (R) such that M A x and M A y are symmetric and M B is skew symmetric. Under condition (3.3), using Fourier analysis, it is then possible, for instance, to obtain polynomial (in time) bounds on the L 2 -norm (in space) of the solution of the Cauchy problem associated to (3.1) (see [START_REF] Kreiss | Initial-boundary value problems and the navier-stokes equations[END_REF]). These bounds are even uniform in time in the case of symmetrisable systems. Indeed, one has the conservation in time of the energy:

E(t) := 1 2 R 2 (M u, u), (3.8)
where (•, •) denotes the inner product in R m . In this case, the perturbation theory for self-adjoint analytic families of matrices (see for instance [18, chapter II, §6]) even asserts that the functions ω j (k) can be chosen in such a way that they are analytic functions of k x and k y (separately) everywhere.

One can see that the Drude model (2.15) is a particular case of such systems and is nondissipative. Indeed, it can be written under the form (3.1) with u = (E x , J x , E y , J y , H, K) T ,

A x = -         0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ε 0 -1 0 0 0 0 0 0 0 0 0 µ 0 -1 0 0 0 0 0 0 0 0 0         , A y =         0 0 0 0 ε 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 µ 0 -1 0 0 0 0 0 0 0 0 0 0 0         and B =         0 ω 2 e 0 0 0 0 -1 0 0 0 0 0 0 0 0 ω 2 e 0 0 0 0 -1 0 0 0 0 0 0 0 0 ω 2 m 0 0 0 0 -1 0         . (3.9)
One easily checks that this system is symmetrisable with the diagonal matrix

M = diag(ε 0 , ε 0 ω 2 e , ε 0 , ε 0 ω 2 e , µ 0 , µ 0 ω 2 m ), (3.10)
and the energy (3.8) is nothing but E e + E m in (2.20). More generally, it is not difficult, although tedious, to show that the generalized Lorentz model (2.17) is also symmetrisable and non-dissipative.

Generalized PMLs for dispersive models

From now on, we will only consider dispersive systems of the form (3.1) that are non-dissipative, i.e. they satisfy (3.3). What follows only concerns PMLs in the x-direction, but everything can easily be adapted to PMLs in the y-direction. As announced in the introduction of this section, we consider a generalized complex change of variable under the form (at a given frequency ω)

x -→ X(ω, x) := x + χ(ω) iω x 0 σ(s) ds, (3.11)
where the damping function σ is non-negative and χ is a real-valued function of ω. The classical Bérenger's PML corresponds to the choice χ(ω) = 1 (see for instance [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 frequency domain, the corresponding PML model, i.e. the equation satisfied by

u(x, y, ω) := Fu (X(ω, x), y, ω) , (3.12)
where Fu is the time Fourier transform of u, is simply obtained by changing

∂ x -→ 1 + σχ(ω) iω -1 ∂ x , (3.13)
in the frequency version of the original system (3.1). Here, we assume that the Fourier transform of u can be extended analytically to complex values of x (which could be proven at the price of technical developments which are outside the scope of this paper). We obtain

iω u + B u = 1 + σχ(ω) iω -1 A x ∂ x u + A y ∂ y u. (3.14)
Note that X(ω, x) = x for x < 0 implies that u(x, y, ω) := Fu (x, y, ω) which is nothing but the perfect matching property: by looking at u(x, y, ω) instead of Fu (x, y, ω) we only modify the solution in the half-space x > 0 (see also Remark 3.4). That is why, by abuse of notation, we shall denote in what follows u the inverse time Fourier transform of u(x, y, ω), which coincides with the original solution u in the half-space x < 0.

In what follows, we shall restrict ourselves to functions χ(ω) that mimic (2.16): they are of the form

χ(ω) = 1 + N =1 b a 2 -ω 2 , ( 3.15)
where a and b are non-zero real constants such that 0

< a 2 1 < a 2 2 < • • • < a 2 N .
Notice that χ is chosen to fulfil the following criteria:

-At high frequencies, the system becoming non dispersive (in the sense that the functions ω j (k) behave more and more as functions which are homogeneous of degree 1), there is no reason that χ varies with ω for large ω so a natural criterion is that χ(ω) tends to 1 when ω tends to +∞. -For practical reasons linked to numerical simulations, the function χ(ω) should be the symbol of a local operator in time (i.e. expressed with the help of differential operators and inverse of differential operators) which means that χ(ω) should be a rational fraction in ω. -In order to ensure that the solutions of the time dependent problems remain real-valued, the operator of symbol χ(ω) should be real (in the sense that it transforms real valued functions into real valued functions). -The poles of χ(ω) are real, i.e. the coefficients a are purely real. This appears as a necessary condition for the stability of the corresponding evolution problem (see Remark 3.11). We shall also use in the following analysis the equivalent expression of χ :

χ(ω) = 1 + P χ (ω) Q χ (ω) , ( 3.16)
where

P χ (ω) = N =1 b p = (a 2 p -ω 2 ) and Q χ (ω) = N =1 (a 2 -ω 2 ), (3.17)
and we will denote by ±z , = 1, . . . , N , the zeros of χ. Using (3.11) and coming back to the time domain, we are going to show that (3.14) can be rewritten as an augmented first order system of the form

∂ t u a + B u a + σ S u a = A x ∂ x u a + A y ∂ y u a , (x, y) ∈ R 2 , t ≥ 0, (3.18)
in which the new "augmented" unknown u a takes values in R ma , and B, S, A x and A y are m a × m a real matrices. The dimension m a can be decomposed as m a = m + m ex so u a can be seen as u a = (u, u ex ) T where u, with values in R m , is the inverse time Fourier transform of û defined by (3.12): it corresponds to the original solution in the region where σ is 0. u ex , with values in R mex , is a vector of extra auxiliary variables. More precisely, defining

iω d := σ χ(ω) u, (3.19)
the equation (3.14) can be rewritten

(iω + B) ( u + d) = A x ∂ x u + A y ∂ y ( u + d). (3.20)
From the expression (3.15) of χ, we deduce that

χ(ω) u = u + N =1 b u , ( 3.21)
where the u are defined by (a 2 -ω 2 ) u = u, = 1, . . . , N. (3.22) Using the above, (3.19) gives

i ω d = σ u + N =1 b u , ( 3.23)
so that (3.20) can be rewritten as

i ω u + B( u + d) + σ u + N =1 b u = A x ∂ x u + A y ∂ y u + A y ∂ y d. ( 3
                 i ω u + B( u + d) + σ u + N =1 b u = A x ∂ x u + A y ∂ y u + A y ∂ y d, i ω d -σ u + N =1 b u = 0, i ω u = v , i ω v + a 2 u -u = 0, ( = 1, . . . , N ), (3.25)
where (3.22) was decomposed into the two last equations of (3.24) thanks to the introduction of the v in order to get a first-order system in time. Indeed, (3.25) gives in the time domain the following system: Remark 3.3. Let us emphasize that there is not a unique way to write the augmented system (3.18) from (3.14)-(3.15) (even not a unique choice for the dimension m a ) because there are several possible choices for the additional auxiliary unknowns. In particular, this type of formulation is an alternative to the split-form of PMLs, as it is done in the original work of Bérenger [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]. The one chosen here in (3.26) appears more convenient to us for the mathematical analysis but the conclusions are valid for any equivalent reformulation of this system.

                 ∂ t u + B(u + d) + σ u + N =1 b u = A x ∂ x u + A y ∂ y u + A y ∂ y d, ∂ t d -σ u + N =1 b u = 0, ∂ t u -v = 0, ∂ t v +
Remark 3.4. In [START_REF] Appelö | Perfectly matched layers for hyperbolic systems: general formulation, well-posedness, and stability[END_REF], the question of constructing PMLs for the system (3.1) is attacked in a different way. The authors propose a very general PML model under the form of an augmented PDE system, similar to (3.26), involving ad hoc auxiliary variables and prove a posteriori that the model is perfectly matched (showing in short that the layer does not produce any reflection). They observe a posteriori that their system can be seen as coming from a complex change of variables similar to (3.11), where the function χ would also depend on the tangential Fourier variable k y , implying the apparition of additional y-derivatives in the augmented system.

Modal stability analysis of the generalized PML model

Definition of stability

We are interested in the stability analysis of the system (3.18) (or equivalently (3.26)) when σ is constant (i.e. independent of x). Here again, we use a modal approach to the notion of stability, looking at plane waves

u a (x, y, t) = U a e i(ωt-kxx-kyy) , U a ∈ C ma , k = (k x , k y ) ∈ R 2 , ω ∈ C, (3.27
ω j (k, σ) ∈ C + := z ∈ C : Im z ≥ 0 . (3.29)
For a plane wave (3.27), the condition (3.29) means that it is (in modulus) non increasing in time, and even exponentially decreasing as soon as Im ω j (k, σ) > 0. Since in practice one wants stability for all σ ≥ 0, we shall adopt a slightly more general notion of stability:

Definition 3.6. The family of PML models (3.18) is uniformly stable if, for all σ ≥ 0, for all j ∈ {1, . . . , m a } and for all k ∈ R 2 , one has

ω j (k, σ) ∈ C + .
(3.30) Remark 3.7. In [START_REF] Appelö | Perfectly matched layers for hyperbolic systems: general formulation, well-posedness, and stability[END_REF], the question of the stability of PMLs for (3.1) is also considered from the point of view of the modal analysis. In particular, they provide a criterion to ensure (3.29): this is based on checking that a finite number of functions that appear in the continued fraction expansion of a particular function -coming from the plane wave analysis of the PML system -have a certain sign. It is not clear at all for us that this criterion can be easily used in our case. We have chosen to follow an alternative approach, of perturbative nature.

Reduction of the analysis

Given k ∈ R 2 and writing

U a = U, D, U 1 , . . . , U N , V 1 , . . . , V N T , ( 3.31)
it is easy to see that, after elimination of ( V 1 , . . . , V N ), we are reduced to looking for values of ω ∈ C for which there exists ( U, D, U 1 , . . . , U N ) = 0 such that

                         i ω U + B( U + D) + σ U + N =1 b U + i k x A x U + i k y A y ( U + D) = 0, i ω D -σ U + N =1 b U = 0, -ω 2 U + a 2 U = U, ( = 1, . . . , N ). (3.32) Let us denote Ω(k, σ) := {ω j (k, σ), 1 ≤ j ≤ m a }. (3.33)
This is also the set of values of ω for which (3.32) has a non-trivial solution. We now describe the structure of Ω(k, σ) and give a new form of the stability criterion (3.30) that will help us to decompose the stability analysis into several steps. First we need to introduce several sets: A := {±a 1 , . . . , ±a N } the set of the poles of χ, Ω c (σ) the set of critical frequencies defined as

Ω c (σ) := ω ∈ C * \ A : f (ω) := 1 + σχ(ω) i ω = 0 , ( 3.34)
and finally Ω PML (k, σ) the set of the solutions ω / ∈ A∪Ω c (σ) of the PML dispersion relation defined by

(ω, k) ∈ Ω PML (k, σ) ⇐⇒ F PML (ω, k, σ) = 0, (3.35)
where F PML is defined by

F PML (ω, k, σ) := F ω, k x f (ω) -1 , k y = det(iωI + B + ik x f (ω) -1 A x + ik y A y ), (3.36)
with F given by (3.6) (coming from the dispersion relation (3.5) of the original system (3.1)) and f is defined in (3.34).

Lemma 3.8. The set Ω(k, σ) can be decomposed as the disjoin union

Ω(k, σ) = Ω PML (k, σ) ∪ Ω c (σ) ∪ E(k) if Ker k x A x = {0} Ω PML (k, σ) ∪ E(k) if Ker k x A x = {0} (3.37)
where Ω PML (k, σ) and Ω c (σ) are defined by (3.35) and (3.34), and where E(k) ⊂ R is a set of real additional modes (described in the proof).

Proof. Adding the first two equations of (3.32) leads to the equivalent system:

             i ω + B ( U + D) + i k x A x U + i k y A y U + i k y A y D = 0, i ω D -σ U + N =1 b U = 0, -ω 2 U + a 2 U = U, ( = 1, . . . , N ).
(3.38)

We next distinguish three cases:

1. ω / ∈ A ∪ {0}. The last two equations of (3.38) give U = U/(a 2 -ω 2 ) and i ω D = σ χ(ω) U. Substituting this into the first equation of (3.38), we obtain using f (ω) defined in (3.34)

i ω + B f (ω) U + i k x A x U + i k y f (ω)A y U = 0.
(3.39)

Since we look for non trivial solutions, For the application to the dispersive models considered in Section 2, the reader will easily check that we are in the case Ker A x = {0}, so the set Ω c (σ) has to be taken into account. We will restrict ourselves to this case in the following. Since E(k) ⊂ R, according to Lemma 3.8 and Definition 3.6, the analysis can be reduced to the following criterion: the family of PML models (3.18)) is uniformly stable, if and only if, for all σ ≥ 0 one has

U = 0. If ω / ∈ Ω c (σ), f (ω) = 0 thus i ω + B U = -i k x f (ω) -1 A x U -i k y A y U, i.e. ω ∈ Ω PML (k, σ). Now if ω ∈ Ω c (σ), then k x A x U = 0, i.e. Ker k x A x = {0}.
Ω c (σ) ∪ Ω PML (σ) ⊂ C + , ( 3.43)
where Ω PML (σ) :=

k∈R 2 Ω PML (k, σ).
(3.44)

Guideline for the stability analysis

Before going into details of the analysis of these sets, we first describe the type of arguments that will be used several times in the rest of the paper. Each branch ω(k, σ) of solutions of the augmented dispersion relation (3.28) composing the sets of (3.43) is continuous with respect to σ and verifies ω(k, 0) ∈ R for all k ∈ R 2 . In order to prove uniform stability, we aim to show one of the following assertions holds:

1. ω(k, σ) ∈ R for all k ∈ R 2 and for all σ ≥ 0 (therefore it belongs to C + ); 2. ω(•, σ) is not real for σ > 0 and belongs to C + for small σ > 0.

In the last case, since the imaginary part of ω(•, σ) cannot change sign because of the continuity of ω(•, σ) with respect to σ, we can conclude that ω(k, σ) ∈ C + for all k ∈ R 2 and for all σ ≥ 0.

Study of the set Ω c (σ)

The set Ω c (σ) is completely independent of the original physical domain and only depends on χ(ω). We now look for the necessary and sufficient condition for the inclusion Ω c (σ) ⊂ C + . The first thing to notice is that Definition (3.34) of Ω c (σ) can be rewritten, using the expression (3.16) of χ, as a polynomial equation in ω of degree 2N + 1, namely

(iω + σ) Q χ (ω) + σP χ (ω) = 0.
(3.45)

It admits 2N + 1 solutions, denoted by ω c (σ) ( = -N, . . . , N ), that are smooth (in particular continuous) functions of σ. We now follow the guideline of Section 3.3.3.

Lemma 3.9. For σ > 0, the solutions of (3.45) cannot be real, except for the case χ(0) = 0 for which ω = 0 is the only real root of multiplicity 1.

Proof. Assume that (iω + σ)Q χ (ω) + σP χ (ω) has a real root ω. Then, since P χ and Q χ have real coefficients, taking the imaginary part leads to ωQ χ (ω) = 0, so ω = 0 or Q χ (ω) = 0 (note that Q χ (0) = 0 since the a do not vanish). Now taking the real part gives (Q χ + P χ )(ω) = 0. If Q χ (ω) = 0, then ω = a or -a for one ∈ {1, . . . , N }, but P χ (±a ) = 0, which contradicts (Q χ + P χ )(ω) = 0. If now ω = 0, then (Q χ + P χ )(0) = 0 that can be written, since Q χ (0) = 0, as

Q χ (0)[1 + P χ (0)/Q χ (0)] = Q χ (0)χ(0) = 0.
So this is absurd, except when χ(0) = 0. But in this case it is easy to see that the constant term of (iω + σ)Q χ (ω) + σP χ (ω) is σQ χ (0)χ(0) = 0 and the first order term in ω is i Q χ (0) = 0, i.e. ω = 0 is a root of multiplicity 1.

It is then sufficient to study the sign of the imaginary part of ω c (σ), -N ≤ ≤ N , for small σ > 0 to conclude. For each branch of solution ω c (σ) (we omit here the subscript), using a Taylor expansion and denoting ω * := ω c (0) ∈ R, one gets

ω c (σ) = ω * + (ω c ) (0) σ + O(σ 2 ), (3.48) which implies that Im ω c (σ) = σ Im(ω c ) (0) + O(σ 2
). Injecting (3.48) into (3.45) and using i ω * Q χ (ω * ) = 0 leads, by identifying the first order term in σ, to

i ω * Q χ (ω * ) + i Q χ (ω * ) (ω c ) (0) + P χ (ω * ) + Q χ (ω * ) = 0. (3.49) -For ω c (σ) = ω c (σ), = 0, one has ω * = ±a , therefore Q χ (ω * ) = 0 and (3.49) gives (ω c ) (0) = -P χ (±a ) ±i a Q χ (±a ) = -i b 2a 2 , (3.50)
whose imaginary part is positive, if and only if, b < 0.

-For

ω c (σ) = ω c 0 (σ) one has ω * = 0 so, since Q χ (0) = 0, (3.49) becomes (ω c 0 ) (0) = i (1 + P χ (0)/Q χ (0)) = i χ(0). If χ(0) = 0, it is clear that Im ω j (0) > 0 if,
and only if, χ(0) > 0. If χ(0) = 0, Lemma 3.9 showed that ω 0 (σ) = 0 ∈ C + for all σ ≥ 0. Remark 3.11. Note that, because of (3.47), the fact that the a are real appears as a necessary condition for the inclusion Ω c (σ) ⊂ C + . More generally, if χ is a rational fraction in iω with real coefficients, there will be some solutions that approach the poles of χ when σ tends to 0. Thus it is necessary that these poles have a positive imaginary part (in the case of even rational fractions this is equivalent to the fact that all the poles are real).

Study of the set Ω PML (σ) : a necessary stability analysis

We are now interested in the set Ω PML (σ), i.e. the solutions of (3.35). More precisely we will study a particular subset of Ω PML (σ): observing that, for fixed ω /

∈ Ω c (σ) and k ∈ R 2 , F PML (ω, k, σ) tends to F (ω, k x , k y ) when σ tends to 0 (where F is defined by (3.6) and F PML by (3.36)). By continuity, there exists m branches of solutions of (3.35), let us say by convention the m first ones ω j (k, σ), 1 ≤ j ≤ m, the so-called "physical branches", such that lim

σ 0 ω j (k, σ) = ω j (k), 1 ≤ j ≤ m. (3.51)
where the ω j (•) are the solutions of the original dispersion relation (3.5). In what follows we make the assumption that these are simple:

ω j (k) = ω m (k)
, for all m = j and for all k ∈ R 2 .

(3.52) However, we think that the following result remains valid without this assumption.

Proposition 3.12. Under assumption (3.52), a necessary condition for Ω PML (σ) ⊂ C + for all σ > 0 (thus for the uniform stability of the PML model (3.18)) is that for all j ∈ {1, . . . , m} and for all k ∈ R 2 , one has

v g,x ω j (k) v p,x ω j (k) χ ω j (k) ≥ 0. (3.53)
where v p (ω j (k)) and v g (ω j (k)) are the phase and group velocities (3.7) and the subscript x refers to their x component.

Proof. We work with any branch ω j (k, σ) satisfying (3.51). Using ω j (k, 0) = ω j (k) ∈ R and the Taylor expansion for small values of σ > 0

ω j (k, σ) = ω j (k, 0) + σ ∂ω j ∂σ (k, 0) + O(σ 2 ), (3.54)
the uniform stability condition (3.30) implies in particular that Im ∂ω j ∂σ (k, 0) ≥ 0, for all k ∈ R 2 .

(3.55)

The quantity ∂ω j /∂σ can be computed using the implicit function theorem. More precisely, injecting (3.54) into (3.35), using F (ω j (k), k) = 0 and looking at the O(σ) term leads to

∂ω j ∂σ (k, 0) ∂F ∂ω (ω j (k), k) + ik x χ ω j (k) ω j (k) ∂F ∂k x (ω j (k), k) = 0, (3.56)
that is to say, since, by (3.52)

∂F ∂ω ω j (k), k = 0, ∂ω j ∂σ (k, 0) = -i ∂F ∂ω (ω j (k), k) -1 k x χ ω j (k) ω j (k) ∂F ∂k x (ω j (k), k). (3.57)
By definition of the phase and group velocities (3.7), we have

v p,x (ω j (k)) = k x |k| ω j (k) |k| , (3.58) and v g,x (ω j (k)) = - ∂F ∂ω (ω j (k), k) -1 ∂F ∂k x (ω j (k), k), (3.59)
where the second identity is obtained by using the implicit function theorem to the dispersion relation (3.5) seen as an equation in ω. Finally, we obtain

∂ω j ∂σ (k, 0) = i v p,x (ω j (k)) v g,x (ω j (k)) χ ω j (k) |k| 2 ω j (k) 2 , (3.60)
from which it is easy to conclude.

For classical PMLs corresponding to χ(ω) = 1, the condition (3.53) can be interpreted through a directional version of backward mode (that differs from Definition 2.12 only for anisotropic media):

Definition 3.13. For a given vector n ∈ R 2 , a mode ω(k) ∈ R is a forward mode in the direction n if (v p (ω(k))•n)(v g (ω(k))•n) > 0 and a backward mode in the direction n if (v p (ω(k))•n)(v g (ω(k))• n) < 0.
Thus, Proposition 3.12 means classical PML models can be uniformly stable only if the original model does not support backward modes in the direction x. In particular, classical PMLs can not be stable for the Drude model. This generalizes the result of [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF] to dispersive models.

Numerical results

Let us confirm this by performing a FDTD numerical simulation using the Drude model. The computational domain is the square Ω = [-20, 20] 2 and the physical domain is [-17, 17] 2 (i. e. the PML has width 3). The coefficient σ is variable and grows quadratically from 0 to 1 in the layer (from inside to outside). We choose ε 0 = µ 0 = 1 and take ω e = ω m = 2. We add at the right hand side of (2.15) (the one for H) a source term f (x, y, t) = g(|x -y|)h(t) with g(r) = e -5r 2 and h(t) = -20(t -1)e -10(t-1) 2 .

(3.61)

The Figure 3 shows snapshots of the field H at different times. We clearly observe instabilities. One can notice that the backward waves are responsible of the instabilities. Indeed, the forward waves -the faster ones -are well absorbed by the PML but as soon as the backward waves reach the PML, the solution grows exponentially in time. The Figure 4 illustrates this statement by plotting in log-log scale the energy E e + E m defined by (2.20). More precisely, one can deduce from the Proposition 2.9 that the group velocity of the forward modes v g (±ω 1 ) is bounded by c = 1 but the group velocity of the backward modes v g (±ω 2 ) is bounded by c/2 = 0.5. That is why instabilities can only occur after a time t = 17/0.5 = 34 s.

Remark 3.14. The reader will easily realize that, when |k| tends to +∞, one recovers the stability criterion of classical PMLs for the non dispersive hyperbolic system corresponding to B = 0. In other words, the introduction of the function χ(ω) cannot produce a stable PML when the classical PML for the non dispersive system already fails to be stable, due to particular anisotropy effects as it occurs for instance in elastodynamics [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF]. It may help to only solve the stability problems due to dispersive effects (i.e. the presence of B).

Summary of the results

Let us now recap the stability results that we obtained:

1. a necessary and sufficient condition for the uniform stability of the PML system (3.18) is that Ω c (σ) ∪ Ω PML (σ) ⊂ C + for all σ ≥ 0;

2. the condition Ω c (σ) ⊂ C + is equivalent to (3.46); 3. a necessary condition to ensure that Ω PML (σ) ⊂ C + is given by (3.53).

It is now natural to look for sufficient conditions to ensure that Ω PML (σ) ⊂ C + . This is possible to do in the case of isotropic dispersive Maxwell equations, which was our original motivation (see Section 2) and is precisely the object of the next section.

Application to isotropic dispersive Maxwell equations

In this section, we pursue our analysis of the generalized PML models introduced in Section 3.2 in the particular case of isotropic dispersive Maxwell equations of the Section 2 that enter the framework of Section 3.1. We have seen in Section 3.3 that (3.46) is a necessary and sufficient condition for Ω c (σ) ⊂ C + . According to (3.43), it remains to seek for the conditions to ensure Ω PML (σ) ⊂ C + . Proposition 3.12 gives a necessary condition (3.53) for Ω PML (σ) ⊂ C + , by analysing a particular subset of Ω PML (σ), the m "physical" branches ω j (k, σ), for small values of σ (see Section 3.3.5). In order to obtain a necessary and sufficient stability condition, it remains to:

1. describe all the branches ω j (k, σ) of Ω PML (σ), for j ∈ {1, . . . , m a } (in particular for j ∈ {m + 1, . . . , m a })

2. show that each branch ω j (k, σ) for j ∈ {1, . . . , m a }, stays in the upper complex plane C + for all values of σ (not only small).

In the two next sections (Sections 4.1 and 4.2), we describe the set Ω PML (σ) and express the necessary stability condition obtained in Section 3.2 in this context. In Section 4.3, we establish necessary and sufficient conditions for the uniform stability of PMLs in the case of isotropic dispersive Maxwell models satisfying the growing property (2.21). We finally extend in Section 4.4 the necessary and sufficient conditions to general isotropic dispersive Maxwell models.

Description of the set Ω PML (σ).

We start here by identifying the m a branches ω j (k, σ) of Ω PML (σ), for j ∈ {1, . . . , m a } and k ∈ R 2 . According to (3.36), the PML dispersion equation is obtained from the one of the original system (2.27), as follows:

ωε(ω) ω 2 ε(ω)µ(ω) -1 + σχ(ω) iω -2 k 2 x -k 2 y = 0. ( 4.1)
Let us remind that the size of the original Maxwell system is m = 2 M + 2M s + 1 (section 2.2.1), and that the corresponding modes are the 2 M propagative modes associated solutions of the dispersion relation (2.29), relabeled as (2.30), and 2M e + 1 non-propagative modes corresponding to the zeros of ωε(ω) (relabeled as (2.28)). We easily see on (4.1) that these last modes are not changed by the presence of the PML and thus do not play any role. Therefore, there only remains to study the branches satisfying

σχ(ω) + iω 2 ω 2 ε(ω)µ(ω) -k 2 y = -ω 2 k 2 x . (4.2)
This equation must be read as a polynomial equation since ε(ω

) = P e (ω)/Q e (ω), µ(ω) = P m (ω)/Q m (ω)
and χ(ω) = 1 + P χ (ω)/Q χ (ω) are rational fractions. Thus, denoting

P em = P e P m , Q em = Q e Q m , then (4.2)
can be rewritten as

[σ(P χ + Q χ )(ω) + iωQ χ (ω)] 2 ω 2 P em (ω) -k 2 y Q em (ω) + ω 2 k 2 x Q 2 χ Q em (ω) = 0. (4.3)
This is a polynomial equation in ω of degree 2( M + 2N + 1). Seen as a polynomial equation in iω (remember that ε(ω), µ(ω) and χ(ω) are even functions of ω), (4.3) has real coefficients so ω is a solution of (4.3) if, and only if, -ω is a solution of (4.3). When σ tends to 0, (4.3) becomes

ω 2 Q 2 χ (ω) ω 2 P em (ω) -Q em (ω)(k 2 x + k 2 y ) = 0, (4.4)
that we can write formally under the form

ω 2 Q 2 χ (ω) ω 2 ε(ω)µ(ω) -|k| 2 = 0, (4.5)
where one recognizes the dispersion relation (2.29) of the original equation. Equations (4.4) and (4.5) show that the 2( M + 2N + 1) branches of solutions of (4.3) can be labelled as follows:

-We first denote ω j (k, σ) and -ω j (k, σ) where 1 ≤ j ≤ M the M = 2 M "physical" branches of solutions such that lim

σ 0 ω j (k, σ) = lim σ 0 ω j (k, σ) = ω j (k), (1 ≤ j ≤ M ), (4.6)
where ±ω j (k), 1 ≤ j ≤ M are the M solutions of the reduced dispersion relation (2.29).

-We denote by ω ± (k, σ) and -ω ± (k, σ), where 1 ≤ ≤ N , the 4N branches such that:

lim σ 0 ω ± (k, σ) = lim σ 0 ω ± (k, σ) = ±a , (1 ≤ ≤ N ). (4.7)
-The 2 remaining branches are denoted ω 0 (k, σ) and -ω 0 (k, σ); they verify

lim σ 0 ω 0 (k, σ) = lim σ 0 ω 0 (k, σ) = 0. ( 4.8)
Using the fact that z and -z have the same imaginary part, proving that Ω PML (σ) ⊂ C + amounts to proving that for all σ > 0 one has (1) ω j (k, σ) ∈ C + for all k ∈ R 2 and for all j ∈ {1, . . . , M };

(2) ω ± (k, σ) ∈ C + for all k ∈ R 2 and for all ∈ {1, . . . , N };

(3) ω 0 (k, σ) ∈ C + .

(4.9)

Reinterpretation of the necessary condition (3.53) for isotropic dispersive Maxwell models satisfying the growing property

From now on we assume the growing property (2.21) to be satisfied unless otherwise specified. First let us notice that, the assumption (3.52) on the simplicity of the solutions of the dispersion relation is automatically satisfied. The necessary condition (3.53) is equivalent to say that the m branches of solutions ω j (k, σ) go to C + for small σ > 0, for j = 1, . . . , m. Since the non-propagative modes do not play any role, this is equivalent to restrict this condition to j = 1, . . . , M (or to j = 1, . . . , M , see the previous section). Proof. First notice that, as a consequence of (2.48) where the growing property (2.21) is used, D(ω) = 0 for all ω ∈ S. Therefore, from (2.39), the necessary condition (3.53) is equivalent to

k 2 x D(ω j (k)) χ ω j (k) ≥ 0, (4.12)
We know (Section 3.3.4) that (4.16b) is a necessary and sufficient condition to ensure Ω c (σ) ⊂ C + for all σ ≥ 0. Moreover in Section 4.2 we saw that (4.16a) is a necessary condition. Thus it remains to show that (4.16a)-(4.16b) are sufficient.

In the rest of this section, we assume that (2.21), (4.16a) and (4.16b) hold. The proof reduces to show that all the branches of Ω PML (σ) have a positive imaginary part. To do so, we want to follow the guideline given in Section 3.3.3 but as we will see the cases k x = 0 and k y = 0 play particular roles. We will do the analysis for k belonging to each of the four quadrants Λ ± ± := R * ± × R * ± delimited by k x = 0 and k y = 0. More precisely we will follow these steps: (step 1) we seek for real solutions when σ > 0 (Lemma 4.6): it cannot happens except when k x k y = 0 or χ(0)Q em (0) = 0; (step 2) we show that for all σ > 0 and for all k in each of the four quadrants Λ ± ± , either a branch ω(k, σ) is equal to zero (corresponding to the case χ(0)Q em (0) = 0) or Im ω(k, σ) = 0 and keeps a constant sign (Lemma 4.7); (step 3) an asymptotic analysis (Lemmas 4.8 and 4.9) shows that Im ω(k, σ) is positive for small σ ≥ 0 and for all k in each of the four quadrants Λ ± ± ; (step 4) we conclude by a continuity argument: since the imaginary parts of the all the branches ω(k, σ) are positive for all σ ≥ 0 and for all k in each of the four quadrants Λ ± ± , by continuity with respect to k, it implies that these are positive for all σ ≥ 0 and for all k ∈ R 2 . We start by seeking the real solutions of (4.3) (step 1): Lemma 4.6. Let σ > 0 be fixed. Then (4.3) has a real solution if and only if k y Q em (0)χ(0) = 0 (that corresponds to the solution ω = 0) or k x = 0.

Proof. First, we start to check if ω = 0 can be a solution of (4.3). Taking ω = 0 directly gives -σ 2 k 2 y (P χ + Q χ )(0)Q em (0). Since σ = 0 and Q χ (0) = 0 (because the a are not zero), it shows that ω = 0 is a solution only if Q em (0) = 0, k y = 0 or χ(0) = Q χ (0) -1 (P χ + Q χ )(0) = 0. Now, assume that ω ∈ R is a real solution of (4.3). Since all the polynomials involved have real coefficients, identifying real and imaginary part leads to

(ω 2 P em (ω) -k 2 y Q em (ω)) σ 2 (P χ + Q χ ) 2 (ω) -ω 2 Q 2 χ (ω) + ω 2 k 2 x Q 2 χ Q em (ω) = 0, (4.17)
and ωQ χ (ω)

ω 2 P em (ω) -k 2 y Q em (ω) Q χ (ω) + P χ (ω) = 0. ( 4.18)
We need to distinguish two cases: (a) ω 2 k 2 x Q 2 χ Q em (ω) = 0. Then (4.17) implies ω 2 P em (ω) -k 2 y Q em (ω) = 0. From (4.18) one has Q χ (ω) + P χ (ω) = 0, because ωQ χ (ω) = 0. This means that χ(ω) = 0, therefore Lemma 4.4 implies that ω ∈ G. On the other hand, (4.17) implies in this case

ω 2 P em (ω) -Q em (ω)[k 2 y + k 2 y ] = 0, i.e. ω ∈ S. This case is thus impossible, since G ∩ S = ∅. (b) ω 2 k 2 x Q 2 χ (ω)Q em (ω) = 0. Then either (i) k x = 0, (ii) ω = 0 (which is possible only if k y Q em (0)χ(0) = 0) or (iii) Q 2 χ (ω)Q em (ω) = 0.
It is easy to check (using the same kind of arguments than in (a)) that the last case (iii) is impossible.

Thanks to this results, we can deduce that the imaginary parts of the solutions do not change sign on each of the four quadrants Λ ± ± for all σ > 0, thus the sign is given by an asymptotic analysis for small σ > 0 (step 2): Lemma 4.7. Let k → ω(k, σ) be any branches of solution of (4.3). Then, for all σ > 0 and for all k in one of the four quadrants Λ ± ± , either ω(k, σ) = 0 (possible only if χ(0)Q em (0) = 0) or Im ω(k, σ) keeps a non-zero constant sign.

Proof. By continuity of ω(k, σ) with respect to k and σ, if its imaginary part changes sign on one of the Λ ± ± it means that there exists k ∈ Λ ± ± and σ such that ω(k , σ ) ∈ R. From Lemma 4.6, this is possible only if Q em (0)χ(0) = 0. If Q em (0)χ(0) = 0, we immediately deduce that Im ω(k, σ) keeps a non-zero constant sign. Assume now that Q em (0)χ(0) = 0. Actually, we detail the proof for χ(0) = 0, the proof being similar for Q em (0) = 0. Since χ only depends on ω 2 , one can factorise χ(ω) = ω 2 χ(ω) and we assume that χ(0) = 0. It is easy to show in this case that if ω satisfies (4.3), then it cannot be a zero of Q χ : from Lemma 4.4 it would imply that ω ∈ G, on the other hand (4.3) would imply ω 2 P em (ω) -k 2 y Q em (ω) = 0 which implies that ω ∈ S. But ω ∈ S ∩ G = ∅ is impossible. Therefore, we can divide by Q χ (ω) in (4.3):

ω 2 [σω χ(ω) + i] 2 ω 2 P em (ω) -k 2 y Q em (ω) + k 2 x Q em (ω) = 0. (4.19)
We deduce that ω(k, σ) = 0 is a solution of order 2. The other branches satisfy

[σω χ(ω) + i] 2 ω 2 P em (ω) -k 2 y Q em (ω) + k 2 x Q em (ω) = 0. (4.20)
With similar arguments that used in the proof of Lemma 4.6, we can easily show that this equation cannot have real solution, therefore the other branches do not cross the real axis.

We thus only need to perform an asymptotic analysis for small σ to conclude (step 3). We already know from Proposition 4.2 that (4.16a) ensures that the "physical " branches (4.9)-( 1) have a positive imaginary part for small σ > 0. Thus it remains to analyse the other branches (4.9)-( 2) and (4.9)-(3) in the next two lemmas. Proof. First notice that, according to the proof of Lemma 4.7, if χ(0) = 0, we have ω 0 (k, σ) = 0 for all k ∈ R 2 and σ > 0. We assume now that χ(0) = 0. Injecting the Taylor expansion

ω 0 (k, σ) = ω 0 (k, σ) =0 +σ ∂ ω 0 ∂σ (k, 0) + O(σ 2 ), (4.21)
for small values of σ > 0 into (4.3) gives (after easy computation)

∂ ω 0 ∂σ (k, 0) = χ(0)k y |k| 2 (ik y ± k x ). ( 4

.22)

Notice that |k| = 0 on Λ ± ± . We deduce

Im ω 0 (k, σ) ∼ σ→0 σ k 2 y χ(0) |k| 2 , (4.23)
which is clearly a positive quantity since k y = 0 and (4.16b) is satisfied.

Lemma 4.9. Assume (4.16b). Then for small values of σ > 0 and for all k in one of the four quadrants Λ ± ± , one has ω ± (k, σ) ∈ C + for all ∈ {1, . . . , N }.

Proof. We only treat the case where ω (k, σ) tends to a when σ tends to 0, the case where it tends to a -is similar. From Lemma 4.7, we know that, for σ > 0, ω ± (k, σ) is not real, thus this is not a pole of χ(ω) or ε(ω)µ(ω): we can work with (4.2) instead of (4.3). The idea is then to get an approximation of this equation when ω is close to a , and to deduce from this approximate equation the behaviour of ω (k, σ) when σ is small. When ω tends to a , ω 2 ε(ω) µ(ω) tends to D := a 2 ε(a )µ(a ) and χ(ω) ∼ b /2a (a -ω). Thus, an approximate equation for (4.2) is

σ b 2a (a -ω) + ia 2 D -k 2 y + a 2 k 2 x = 0. ( 4
σ b 2a (a -ω) + ia = ± X with X = a 2 k 2 x k 2 y -D , ( 4.25)
that is to say

ω = a - σ b 2 a (± X -ia ) . ( 4.26)
In particular, since X is real we have

Im ω ∼ ω→a - σ b 2(X 2 + a 2 ) . ( 4.27)
Then we conclude that Im ω > 0 for small σ > 0 since b < 0.

We can now conclude the proof of Theorem 4.5 (step 4). Combining Lemma 4.7 with Proposition 4.2 and Lemmas 4.8-4.9 shows that the imaginary parts of the branches ω(k, σ) of solutions of (4.3) are positive for all σ ≥ 0 and for all k in each of the four quadrants Λ ± ± . By continuity, this implies that these are positive for all σ ≥ 0 and for all k ∈ R 2 . Theorem 4.5 is proven.

Extension of the necessary and sufficient stability conditions for general isotropic dispersive Maxwell models

We will now discuss briefly what happens to Theorem 4.5 when the growing property (2.21) is not satisfied any more. Two of the key points of the previous proof are Lemma 4.4 and the fact that D(ω) = 0 for all ω ∈ S. If the growing property does not hold, these two points might not be true any more so one must add additional assumptions to ensure stability. Doing so, we have the following result that extends Theorem 4.5: ± with k y small enough and σ > 0 small enough. In particular, it is negative if (4.28d) holds.

Proof. The beginning of the proof is the same as in the proof of Lemma 4.9. In particular, (4.24) still holds. Now if a / ∈ G, it means that ε(a )µ(a ) > 0, i.e. D > 0. Taking k y small enough such that D -k 2 y > 0, (4.25) becomes

X = i X , X = |k x a | (D -k 2 y ) 1 2 . (4.29)
At the end, we obtain the equivalent of (4.27) in our case:

Im ω ∼ ω→a - σ b 2a (a ± X ) . ( 4

.30)

Since a < X , one of the two last quantities has the same sign of b .

Lemma 4.12. Let σ > 0 be fixed. If there exists ∈ {1, . . . , N } such that z / ∈ G, then there exists j ∈ {1, . . . , M } such that the imaginary part of ω j (•, σ) changes sign.

Proof. Recall that z /

∈ G means that ε(z )µ(z ) > 0. We only treat the case z = 0, the case z = 0 is similar. Setting

k * x = z ε(z )µ(z ), we have z 2 ε(z )µ(z ) = k * x . It shows that z is a solution of (4.3) (since χ(z ) = 0) for k * = (k *
x , 0). In particular, there exists j ∈ {1, . . . , M } such that ω j (k * , σ) = z . Let study the behaviour of ω(k x ) := ω j (k x , 0, σ) near k *

x (we omit the dependence in j and σ). One has

ω(k x ) = z + dω dk x (k * x )(k x -k * x ) + O[(k x -k * x ) 2 ]. ( 4

.31)

To show that the imaginary part of ω changes sign when k x -k * x does, it is sufficient to show that Im(dω/dk x )(k *

x ) = 0. Injecting the above expansion into (4.3) and using χ(z ) = 0 and

z 2 ε(z )µ(z ) = k *
x gives, after some computations,

Im dω dk x (k * x ) = 2 σ χ (z ) z (k * x ) 3 |Z(z )| 2 (4.32)
where

Z(ω) := ω D(ω) -2i σ (k * x ) 2 χ (ω) ω . ( 4

Design of Stable PMLs for isotropic dispersive media

We consider in this last section the construction of a function χ(ω) fulfilling (4.16) (or (4.28)). Such a function is certainly non unique, as we shall see, and one should privilege, for computational reasons, the ones that minimize N in (3.15), that is to say the number of additional unknowns.

The case of the Drude model

The Drude model is a particular case of an isotropic dispersive Maxwell model satisfying the growing property (2.21). According to the Theorem 4.5, we propose to construct a stable PML by choosing a function χ satisfying D(ω)χ(ω) > 0 for all ω ∈ S, where we have seen S = S + ∪ S -, S -= [-ω -, ω -], S + = R \ (-ω + , ω + ).

(5.1)

For the Drude model, D(ω), given by (2.46) satisfies D(ω) > 0 in S + and D(ω) < 0 in S -. We thus need χ(ω) > 0 in S + and χ(ω) < 0 in S -which is achieved by picking ω * ∈ [ω -, ω + ] and taking

χ(ω) = 1 - ω 2 * ω 2 -1 . ( 5.2)
Finally, we simply have to check (4.16b): one has χ(ω) = 1 -ω 2 * /(ω 2 * -ω 2 ) that corresponds to (3.15) with N = 1, a 1 = ω * and b 1 = -ω 2 * < 0. Notice that we are in the limit case χ(0) = 0. Remark 5.1. One particular choice consists in choosing ω * = ω e or ω * = ω m which amounts to take χ(ω) = ε 0 /ε(ω) or χ(ω) = µ 0 /µ(ω) (this choice is used for instance in [START_REF] Steven | Perfectly matched layer behavior in negative refractive index materials[END_REF][START_REF] Dong | Perfectly matched layer-absorbing boundary condition for left-handed materials[END_REF][START_REF] Shi | Perfectly matched layer absorbing boundary condition for truncating the boundary of the left-handed medium[END_REF]). This is in fact a particular case of a more general result (see Section 5.2). Such a choice is imposed automatically when ω e = ω m .

We performed a numerical simulation to illustrate the stability of the new PML. More precisely, we reproduced exactly the same experiment as in Section 3.3.6 (same parameters, same source, etc.) and choose here ω * = ω e (however the same conclusions as below were also observed with other choices of ω * between ω e and ω m ). As expected, and contrary to what was observed in Section 3.3.6, the solution does not blow up stable during the whole computation with a final time of T = 250 s. The stability is confirmed by the Figure 6 where we show the evolution of the energy E e + E m defined by (2.20). The Figure 5 shows some snapshots of the field H at different times which the shows that the PMLs works well: both forward and backward waves are well absorbed. This is confirmed by a comparison with numerical computations done with a larger computational domain.

The case of generalized Lorentz materials

The following result is a generalization of the Remark 5.1 for generalized Lorentz materials. Between the two possibilities offered by (5.3), one should obviously choose the one that minimizes N , that is to say the number of additional unknown functions for the time domain PML model. Moreover, there is no reason a priori to state that one of the two choices (5.3) is optimal in terms of the minimization the number N ; it is not difficult, by applying the construction process that we describe in the next section, to find particular examples for which an optimal choice is none of them.

Remark 5.3. Noticing that the set of rational functions χ of the form (3.15) satisfying (4.16) is convex, a stable PML is also obtained with

χ(ω) = θ ε 0 ε(ω) + (1 -θ) µ 0 µ(ω)
, for any θ ∈ [0, 1].

(5.4)

A general construction process for optimal stable PMLs

We propose a procedure for the construction of a stable PML for a general dispersive medium based on the properties of the function D(ω), in the sense that χ verifies (4.28). This construction is moreover supposed to be minimal in the sense that the number of poles 2N of the rational function χ is minimal.

5.3.1

Case where the number of changes of sign of D(ω) in S ∩ R + is even.

We suppose that the number of changes of sign of D(ω) inside the set S ∩ R + is equal to 2M . Since D(ω) > 0 for large ω, we assume that there exists 2M numbers 0 < ω 1 < • • • < ω 2M such that (with the convention ω 0 = 0 and ω 2M +1 = +∞) (-1) D(ω) > 0 in (ω 2 , ω 2 +1 ) ∩ S, (0 ≤ ≤ 2M ).

(5.5)

The reader will easily realize if the number M is fixed, there may exist some flexibility in the choice of the ω if there are gaps (G = ∅). A good choice for the function χ(ω) is then: given by (3.15) is stable. Moreover, this additional condition is automatically satisfied for materials possessing the growing property (2.21).

χ(ω) = M =1 1 - ω 2 2 -1 ω 2 M =1 1 - ω 2 2 ω 2 -1 . ( 5

5.3.2

Case where the number of changes of sign of D(ω) in S ∩ R + is odd.

In this case, there exists 2M + 1 numbers 0 < ω 1 < • • • < ω 2M +1 such that (with the convention ω 0 = 0 and ω 2M +2 = +∞) (-1) ( +1) D(ω) > 0 in (ω 2 , ω 2 +1 ) ∩ S, (0 ≤ ≤ 2M + 1).

(5.11)

Then, the reader will easily check that a good choice, in the sense that a result similar to the Theorem 5.4 still holds, for the function χ(ω) is then :

χ(ω) = M =1 1 - ω 2 2 ω 2 M =0 1 - ω 2 2 +1 ω 2 -1
.

(5.12)

Figure 2 -

 2 Figure 2 -the modes ω 1 and ω 2 in terms of |k|

. 24 )

 24 With (3.21), (3.22), (3.23) and (3.24), we finally obtain the system:

a 2 u

 2 u = 0, ( = 1, . . . , N ).(3.26) This is an augmented system of the form(3.18) with m a = m + (2N + 1)m (thus m ex = (2N + 1)m) and u a = (u, u ex ) T where u ex = (d, u 1 , . . . , u N , v 1 , . . . , v N ) T .

2 .

 2 ω ∈ A: ω 2 = a 2 for one ∈ {1, . . . , N }. One easily checks from (3.38) that U = 0, U k = 0 for k = and (±a I + iB -k y A y ) D = 0, U = ±ia σb D. (3.40) This is possible only if ±a is an eigenvalue of iB -k y A y (otherwise D = U = 0). In this case, ±a ∈ E(k) = ∅. 3. ω 2 = 0. From (3.38), U = U/a 2 thus 0) = 0, then U = 0 and we get a non trivial solution only if ker(B + ik y A y ) = {0} (so in this case 0 ∈ E(k)). If now χ(0) = 0, then the second equation of (3.38) is just 0 = 0 and the first equation reads B ( U + D) + i k x A x U + i k y A y U + i k y A y D = 0. (3.42) So if this equation has non trivial solutions, then 0 ∈ E(k) = ∅.

Proposition 3 .

 3 10 (necessary and sufficient condition for Ω c (σ) ⊂ C + ). The two following statements are equivalent: 1. for all σ > 0, one has Ω c (σ) ⊂ C + ; 2. the function χ defined by (3.15) verifies χ(0) = 1 + N =1 b a 2 ≥ 0, and b < 0 for all ∈ {1, . . . , N }. (3.46) Proof. For σ = 0, (3.45) becomes iω Q χ (ω) = 0 whose solutions are ω = 0 and ω = ±a , 1 ≤ ≤ N . By continuity, we can number the ω c such that lim

Figure 3 -

 3 Figure 3 -Snapshots of the field H at different times t.

Figure 4 -

 4 Figure 4 -Logarithm of the energy E e + E m defined by (2.20).

Lemma 4 . 1 .

 41 For a material satisfying the growing property(2.21), all solutions of the dispersion relation (2.29) are simple.Proof. It suffices to show that for any ω = ω j (k), D(ω) = 0. HoweverD(ω) = (ω ε) (ω) ω µ(ω) + (ω µ) (ω) ω ε(ω) .(4.10)Since ω 2 ε(ω) µ(ω) = |k| 2 , ω µ(ω) and ω ε(ω) are non-zero real numbers with the same sign. Since (ω ε) (ω) and (ω µ) (ω) are strictly positive, using the growing property (2.21) it follows from (4.10) that D(ω) has the same sign than ω µ(ω) and ω ε(ω).

Proposition 4 . 2 .

 42 For isotropic dispersive Maxwell models satisfying the growing property (2.21), the necessary stability condition (3.53) is equivalent to D(ω) χ(ω) > 0, for all ω ∈ S. (4.11)

Lemma 4 . 8 .

 48 Assume (4.16b). Then for small values of σ > 0 and for all k in one of the four quadrants Λ ± ± , one has ω 0 (k, σ) ∈ C + .

. 24 )

 24 If (4.16b) holds then a ∈ G (using Lemma 4.4) i.e. ε(a )µ(a ) ≤ 0 or equivalently D ≤ 0. Since D -k 2 y < 0 (because k y = 0), we get

Theorem 4 . 10 .

 410 Consider an isotropic dispersive Maxwell model of Section 2 such that D(ω) = 0 for all ω ∈ S \ A. A necessary and sufficient condition for the uniform stability (Definition 3.6) of the PML model (3.18) associated to the function χ(ω) given by (3.15) is D(ω) χ(ω) > 0 for all ω ∈ S \ A, (4.28a) a ± ∈ G for all ∈ {1, . . . , N }, (4.28b) z ± ∈ G for all ∈ {1, . . . , N }, (4.28c) χ(0) ≥ 0 and b < 0 for all ∈ {1, . . . , N }. (4.28d) Note that recover Theorem 4.5 if the growing property (2.21) holds, since in this case D(ω) = 0 for all ω ∈ O while, by Lemma 4.4, (4.28a) implies (4.28b)-(4.28c) . The additional conditions (4.28b) and (4.28c) allow us to "bypass" Lemma 4.4. Doing so, we can follow exactly the same steps of the proof of Theorem 4.5, using (4.28b) and (4.28c) instead of Lemma 4.4 when needed. That shows that the four conditions (4.28) are sufficient and the conditions (4.28a) and (4.28d) are necessary. It remains to show that (4.28b) and (4.28c) are necessary. This is the object of the two following lemmas that conclude the proof of Theorem 4.10. Lemma 4.11. If it exists ∈ {1, . . . , N } such that a / ∈ G, then the imaginary part of one of the corresponding branches ω ± (k x , k y , σ) has the sign of b for all k in each quadrant Λ ±

  .33) Since z = 0, one has k * x = z 2 ε(z )µ(z ) = 0. Moreover, since the zeros of χ are simple (in the case z = 0), χ (z ) = 0. So |Z(z )| = 0 and Im(dω/dk x )(k * x ) = 0.

Figure 5 -

 5 Figure 5 -Some snapshots of the field H at different times t for the case ω e = ω m .

Figure 6 -Theorem 5 . 2 .

 652 Figure 6 -the energy E e + E m defined by(2.20). As soon as the waves reach the PML, they are well absorbed and the energy is decaying, both for forward and backward waves.

. 6 )Theorem 5 . 4 .

 654 By construction χ satisfies (4.28a). It is a rational fraction in ω 2 , with M simple zeros (the ω2 2 -1 ) and M simple poles (ω 2 2 ). Moreover it tends to 1 at infinity. Consequently, it admits a decomposition of the form χ(ω) . a = ω 2 referring to(3.15). The stability conditions (4.28d) is satisfied, sinceχ(0) to conclude that b < 0 since #{k : ω 2 2 < ω 2 2k-1 } = M -and #{k : ω 2 2 < ω 2 2k } = M -. (5.10) It also shows that (4.28b) is satisfied. The remaining condition (4.28c) is less obvious. That is why we shall restrict ourselves to the following result: Assume that the ω , 1 < < 2M , are chosen in such a way that (5.5) is satisfied. If in addition ε(ω 2 ) µ(ω 2 ) ≤ 0, 1 ≤ ≤ M, the PML model associated to the function χ(ω)

.2.3 Phase and group velocities Definition 2.10. A family

  ω j (k) ≡ ω j (|k|) of solutions of (2.29) indexed by k, where ω j (•) is smooth, is called a branch. The associated plane waves (2.24) are called modes.

								.33)
	where Σ(|k|) := (ω 2 e -ω 2 m + c 2 |k| 2 ) 2 + 4 ω 2 m c 2 |k| 2 .
	Proof. By defining W = ω 2 , the dispersion relation is the quadratic equation W 2 -(ω 2 e + ω 2 m + c 2 |k| 2 )W + ω 2 e ω 2 m = 0. The discriminant of this equation is actually Σ(|k|) 2 . Thus the two solutions
	are real and given by
	W 1 =	ω 2 e + ω 2 m + c 2 |k| 2 + Σ(|k|) 2	and W 2 =	ω 2 e + ω 2 m + c 2 |k| 2 -Σ(|k|) 2	.	(2.34)
	It is clear that W 1 > 0. Notice that the product W 1 W 2 of the zeros of this quadratic equation is
	also equal to ω 2 e ω 2 m > 0, then necessarily W 2 > 0. Taking the square roots of W 1 and W 2 gives the
	result.
	2Definition 2.11. For a mode associated to a branch ω(k), its phase and group velocities v p (ω(k))
	and v g (ω(k)) are defined as
	v p (ω(k)) :=	ω |k| |k|	k |k|	and v g (ω(k)) := ∇ k ω(k) = ω |k|	k |k|	.	(2.35)

  ) which are a solution of (3.18) if, and only if, (ω, k x , k y ) satisfies the augmented dispersion relationdet i ωI + B + σ S + ik x A x + ik y A y = 0.(3.28)Seen as an equation in ω, (3.28) defines m a branches of solutions ω j (k, σ), 1 ≤ j ≤ m a where the functions ω

j (k, σ) are continuous and piecewise smooth in k and σ (they are piecewise analytic in k x , k y and σ separately [18, chapter II, §2]). Definition 3.5. Given σ > 0, the PML model (3.18) is stable if, for all j ∈ {1, . . . , m a } and for all k ∈ R 2 , one has

The third author was partially supported by the ANR project METAMATH (ANR-11-MONU-0016)

for all j = 1, . . . , M . We can show in the following Lemma 4.3 that in the case k x = 0, the M physical branches are real for all σ and do not play any role in the necessary stability condition. Therefore (4.12) becomes equivalent to D(ω) χ(ω) ≥ 0, for all ω ∈ S.

(4.13)

Since D(ω) = 0, the left hand side could vanish only if χ could have a zero in S. We use Lemma 4.4 to show that this is not possible and conclude the proof.

We now give and prove the two lemmas used in the previous proof.

Lemma 4.3. In the case k x = 0, the solutions of equation (4.3) for all k y ∈ R and σ > 0, are given by:

1. ω j (0, k y , σ) = ω j (0, k y ) ∈ R for all j ∈ {1, . . . , M }, where the ±ω j are the solutions of the original dispersion relation (2.29);

The equation (i) clearly corresponds to the original dispersion relation for k x = 0, whose solutions are real (this gives (1)). We recognize in the equation (ii) the characterization (3.45) of Ω c (σ) whose solutions are the ω c (σ). We conclude using Proposition 3.10 (this gives (2)). Proof. We first notice that S is an open set of R composed of the union of open intervals I j of R.

On each of them, (2.48) (satisfied since we assume (2.21)) shows that ε, µ and D have the same constant sign. In particular, D does not change sign in each I j . Since all the zeros and poles of χ are simple, if one of them was inside one of the I j , it would mean that χ changes sign inside this interval, which contradicts (4.11).

Necessary and sufficient stability conditions for isotropic dispersive Maxwell models satisfying the growing property

We are able to establish necessary and sufficient stability conditions for the PMLs introduced in Section 3.2. This is the main result of this section: