Jonathan Viquerat

Maciej Klemm

Stéphane Lanteri

Claire Scheid

Theoretical and numerical analysis of local dispersion models coupled to a discontinuous Galerkin time-domain method for Maxwell's equations

Keywords: Discontinuous Galerkin method, Maxwell's equations, numerical electromagnetism, dispersive media, nanophotonics, Drude model

This report focuses on a centered-fluxes discontinuous Galerkin method coupled to a second-order Leap-Frog time scheme for the propagation of electromagnetic waves in dispersive media. After a presentation of the physical phenomenon and the classical dispersion models (particularly the Drude one), a generalized dispersive model is introduced. An a priori stability and convergence study is lead for the Drude model, as well as in the generalized dispersive case. Eventually, numerical results are presented for various test-cases, highlighting the interest of a proper description of the dispersion phenomenon in metals at the nanoscale.

Theoretical and numerical analysis of local dispersion models coupled to a discontinuous Galerkin time-domain method for Maxwell's equations

Résumé : Ce rapport présente une méthode de Galerkin discontinue à flux centrés couplée à un schéma d'avancée en temps de type Leap-Frog d'ordre deux pour la propagation des ondes électromagnétiques dans les milieux dispersifs. Après une présentation du phénomène physique ainsi que des modèles de dispersion les plus classiques (notamment celui de Drude), un modèle de dispersion généralisé est introduit. Une étude de stabilité et de convergence a priori est conduite dans le cas du modèle de Drude, ainsi que dans le cas généralisé. Enfin, des résultats numériques sont présentés pour différents cas-tests, mettant en lumière l'intérêt d'une bonne description des phénomènes de dispersion des métaux à l'échelle nanoscopique.

Mots-clés : Méthode de Galerkin discontinue, équations de Maxwell, électromagnétisme numérique, milieux dispersifs, nanophotonique, modèle de Drude.

List of Figures

Photonics and dispersive media

The increasing need in electronic signals processing during the last century came along with an advanced miniaturization of all the associated components. The revolution that followed the advent of the transistor deeply modified the common use of information. In the seventies, the development of the optical fiber led the way to fastest and largest data exchanges by exploiting light flow instead of electron flow. A precise control of the electromagnetic waves (confinment, transmitted frequencies, propagation direction) using wavelength-size systems would permit to replace the electronic-based management of the information. For many fields (e.g. optics, informatics, medicine), the benefits withdrawned from such an outcome would be invaluable.

The particular properties required by such devices are artificial, i.e. they are not displayed by any natural material. The macroscopic behavior of these "metamaterials" usually results from their nanoscopic structure, which induces particular interactions with propagating waves. Under this vocable are included the negative-index materials [START_REF] Veselago | Negative refractive index materials[END_REF], as well as photonic cristals [START_REF] Joannopoulosm | Photonic Crystals, Molding the Flow of Light[END_REF], for example.

Then, the fast-growing metamaterials development is due to their innovative features at the macroscopic scale. However, one must not forget that these properties rely on those of "regular" materials (dielectrics, metals, etc.) assembled at the nanoscopic scale. Hence, poorly modelized attributes for the latter would result in imprecisely predicted properties for the resulting metamaterial. Among others, physical dispersion has a great impact on the material/electromagnetic wave interaction; therefore, attention has been called on the various ways of modeling physical dispersion, especially in metals.

Physical dispersion 1.Definition

Dispersion is a common phenomenon to all kinds of waves traveling through a material medium : it results from the way the latter reacts to the presence of the wave, therefore affecting its propagation. For a polychromatic wave, it often happens that all the frequencies do not travel at the same speed through the medium : this phenomenon is called dispersion. In order to characterize it, one generaly tries to express a dispersion relation that bounds the angular frequency ω to the wave number k.

Dispersion relation

Since it links a temporal parameter ω to a spatial one k, this relation ω = f (k) describes the allowed modes inside the considered domain. It is often written as :

ω(k) = v(k)k, (1)
where v(k) is the speed of the monochromatic wave of wavenumber k. This relation is often explicit in the most simple cases, but it can become implicit in more complicated situations.

Let us now consider a polychromatic wave that can be expanded into a sum of monochromatic modes. In the case of a non-dispersive medium, v(k) would remain constant for every value of k, and all the frequencies would travel at the same speed. On the contrary, if the medium is dispersive, (1) indicates that each frequency travels at a different speed. Therefore, a signal consisting of frequencies that are initially spatially localized will sprawl during its propagation. The description of this phenomenon calls for the definition of two different speeds :

• The speed of the wavefront, usually called phase velocity : v φ = ω(k) k ,

• The speed of the wave envelope, called group velocity : v g = ∂ω ∂k

. These two velocities are equal in the case of a non-dispersive medium. Otherwise, the group velocity is the physical speed of the problem, given it is linked to the speed at which the energy moves through the medium. On the other hand, the phase velocity is more artificial, since it represents the speed at which a point of constant phase moves through the medium. One could imagine situations for which the phase velocity exceeds the speed of light c, without violating any fundamental principle of physics.

We now restrain our study to the case of propagating electromagnetic waves (EMW) through dispersive metals, which is at the heart of to the numerical study that will follow.

Physical origin of the dispersion in metals for EMW

In the presence of a constant electric field, the electrons of a metal are subjected to a Coulomb force which brings them, in a given characteristic time τ c , to an equilibrium position. This leads to a general electric polarization of the metal, which is usually expressed with the polarization vector P. The latter constitutes an additional term to the electric displacement field D : D = ε 0 E + P. Moreover, P can be related to E in homogeneous isotropic media through its susceptibility χ such that P = χE. One should now grasp the importance of taking the dispersion effects into account when P cannot be neglected, since it has a great influence on the permittivity ε of the considered medium, and hence on its optical index.

If one is to consider a variable electric field of given angular frequency ω, the frequency dependence of P can be intuitively understood : for low enough frequencies, the electrons relaxation time τ is negligible compared to 1 ω . Therefore, the electrons dispose of a sufficient amount of time to adapt to the variations of the electric field. However, at higher frequencies, the field varies significantly during the time τ required by the electrons to reach a stable state. Then, the higher the frequency, the shorter the distance traveled by the electrons from their steady state equilibrium, and the lower the polarization. This explains the observed transparency of the metals for very high frequencies EMW. Moreover, one can now easily picture that the electrons contributing to the polarization mainly belong to the conduction band, since they are less bound to the nuclei, and therefore more movable.

It is now possible to introduce the most common modelisations of the disperion phenomenon in metals. Among others, the well-known Drude and Drude-Lorentz models will be presented.

Drude model

The discovery of the electron in 1897 by Thomson was followed, three years later, by the Drude theory based on the kinetic theory of gases [START_REF] Drude | Zur elektronentheorie der metalle[END_REF]. This particularly simple theory successfully accounts for the optical and thermical properties of some metals. In this Inria model, the metal is considered as a static lattice of positive ions immersed in a free electrons gas. Those electrons are considered to be the valence electrons of each metallic atom, that got delocalized when put into contact with the potential produced by the rest of the lattice atoms. They are often called conduction electrons, and wander freely in the conduction band, the core electrons being considered to be strongly bound to the ion core. As a consequence, the Drude model is often described as an intraband model. Moreover, several hypothesis are made :

• The electrons description is non-relativistic;

• The only considered interactions are the electron/wave and the electron/ion ones;

• The electron/ion collisions are instantaneous and random events, and their probability of happening during a dt amount of time is equal to dt τ f1 ;

• After an electron/ion collision, the new velocity and direction of the electron are independant of those before the collision. Under these hypotheses, the frequency dependence of the medium permittivity can be deduced from the equations of motion.

Relative permittivity

In the absence of relativistic effects (as the Laplace force), the velocity equation of an electron can be written as follows :

∂v ∂t + γ d v = - e m e E(t), (2)
where m e represents the electron mass, e the electronic charge, and γ d a coefficient linked to the electron/ion collisions. One could notice that [γ d ] = T -1 , and therefore γ d matches the definition of the inverse of the mean free path τ f . For an harmonic field (i.e. E = E 0 e -iωt ), the electron velocity will be of the form v(t) = v 0 e -iωt . Then, (2) leads to :

v 0 = - eE 0 m e (γ d -iω) . (3)
Combining the current density2 J = -n e ev with (3) and Ohm's law3 J = σ d E, one obtains :

σ d (ω) = n e e 2 m e 1 γ d -iω . ( 4
)
The relative permittivity is linked to its conductivity with the relation ε r

,d (ω) = ε ∞ + iσ d ωε0 .
Here, ε ∞ represents the core electrons contribution, which is equivalent to approximate the positive ions lattice as a continuous medium of permittivity ε ∞ . For the metals successfully described by the Drude model, ε ∞ is close to 1 (taking ε ∞ different from 1 is sometimes refered to as the Drude-Sommerfeld model). Eventually, one gets :

ε r,d (ω) = ε ∞ - ω 2 d ω 2 + iωγ d , (5)
where ω d = nee 2 meε0 is the plasma frequency of the electrons. The latter is linked to a typical relaxation time of slightly perturbed plasma electrons around their thermodynamic equilibrium state. By separating the real and imaginary parts of (5) :

ε r,d = ε 1 + iε 2 = ε ∞ - ω 2 d γ 2 d + ω 2 + i γ d ω 2 d ω (γ 2 d + ω 2 )
, and by considering ω ≪ γ d (which remains particularly accurate in the visible frequency field for the metals described by the Drude model), the above expressions can be simplified as :

ε 1 (ω) ≃ ε ∞ - ω 2 d ω 2 and ε 2 (ω) ≃ ω 2 d γ d ω 3
, which enlights the role of the parameter γ d in the dissipation induced by the Drude model4 .

The real and imaginary parts of the Drude permittivity function of the silver are plotted in figure 1, along with experimental curves from Johnson and Christy [START_REF] Johnson | Optical constants of the noble metals[END_REF]. The used parameters can be found in table 1. One should notice that, if the real part fits the Drude model prediction, the experimental imaginary parts show features that are not predicted by the model. Those latter root in interband phenomenons, i.e. core electrons contributions that jump to the conduction band. These features can be better fitted with the Drude-Lorentz model, which is discussed in the following section.

Drude-Lorentz model

The initial Drude model describes the behavior of the free electrons contained in the conduction band. Nevertheless, for certain metals (such as noble ones), electronic transitions between valence and conduction band occur in the visible frequency range, giving them their particular colors. The contributions of these electrons can be taken into account in the permittivity function with additional Lorentz terms. Contrary to the Drude-described electrons, the Lorentz ones are, in some way, "bound" to their ion cores. Therefore, it seems logical to reuse the expression (2), including an additional spring term to it :

∂ 2 x ∂t 2 + γ l ∂x ∂t + ω 2 l x = - e m e

E(t).

Following the same development as for the Drude model, one easily obtains the expression of the Lorentz permittivity :

ε r,l (ω) = - ∆εω 2 l ω 2 -ω 2 l + iωγ l .
The total permittivity can then be written by adding the Drude and Lorentz terms5 :

ε r,dl (ω) = ε ∞ - ω 2 d ω 2 + iωγ d - ∆εω 2 l ω 2 -ω 2 l + iωγ l . ( 6
)
As done previously, the previous permittivity can be splitted into a real and an imaginary parts :

ε 1 (ω) = ε ∞ - ω 2 d γ 2 d + ω 2 - ∆εω 2 l ω 2 -ω 2 l (ω 2 -ω 2 l ) 2 + γ 2 l ω 2 , and ε 2 (ω) = γ d ω 2 d ω (γ 2 d + ω 2 ) + ∆εω 2 l γ l ω (ω 2 -ω 2 l ) 2 + γ 2 l ω 2 .
For some metals such as gold, silver or copper, the addition of Lorentz terms brings a much better fit between experimental and theoretical values. In the case of silver, a Drude-Lorentz model with one Lorentz pole brings a much better adequation between experiment and theory, especially in the high-frequency range for the imaginary part, as can be seen on figure 3. The Drude-Lorentz model parameters used for this fitting are summed up in table 2. Nevertheless, it is to be understood that Drude and Drude-Lorentz parameters are fitted from experimental data over a particular range of frequencies. This has several consequences, such as (i) the variety of parameter sets that can be found in the literature for the different models (particular attention must be paid to the frequency range of interest); (ii) in regions of lower regularity of the permittivity function, additional adequatly-chosen poles can be considered to improve the fitting (the next section presents a way of fitting arbitrary permittivity functions).

Generalized dispersive model

Given an experimental permittivity function, a Padé approximant would be a convenient analytical coefficient-based function to approach the experimental data. Thus, one could write :

ε r,g (ω) = n i=1 α i (jω) i m k=1 α k (jω) k .
The fundamental theorem of algebra allows to write the previous expression as a decomposition of a constant, one zero-order pole (ZOP), a set of first-order generalized poles (FOGP), and a set of second-order generalized poles (SOGP) :

Inria ε r,g (ω) = ε ∞ - σ jω - l∈L1 a l jω -b l - l∈L2 c l -jωd l ω 2 -e l + jωf l , (7)
where ε ∞ , σ, a l , b l , c l , d l , e l , f l are real constants, and L 1 , L 2 are non-overlapping set of indices. This general writing allows an important flexibility for two reasons : (i) it unifies most of the common dispersion models in a single formulation (the right coefficient choices are presented in table 3 for the Drude, Drude-Lorentz and Debye media), and (ii) it permits to fit any experimental data set in a reasonable number of poles (and thus a reasonable number of coefficients), as will be presented later. One should notice that this approach is very similar to the Critical Points (CP) [START_REF] Vial | A new model of dispersion for metals leading to a more accurate modeling of plasmonic structures using the FDTD method[END_REF] and the Complex-Conjugate Pole-Residue Pairs (CCPRP) [START_REF] Han | Model dispersive media in FDTD method with complex-conjudate pole-residue pairs[END_REF]. In facts, these developments are the same in some mathematical sense, since the three of them exploit the fundamental theorem of algebra. The difference between them lies in the way they choose to use it : the CCPRP and the CP allow complex coefficients in their developments, and can therefore write the decomposition of the permittivity function in terms of single-order poles only, whereas choosing real coefficients leads to a collection of first-order and second-order poles. Model

a l b l c l d l e l f l Drude 0 0 ω 2 d 0 0 γ d Drude-Lorentz 0 0 ∆εω 2 l 0 ω 2 l γ l Debye ∆εγ de -γ de 0 0 0 0
The matter of fitting the coefficients of (7) to experimental data remains to be detailed. Various techniques can be used, through which the many existing versions of least square method. Vector fitting techniques [START_REF] Gustavsen | Rational approximation of frequency domain responses by vector fitting[END_REF] are also well developed for the CCPRP formulation. In our case, a Simulated Annealing (SA) algorithm has been used to determine the coefficients to use given an experimental data set. The original paper describing this method is [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF], but the reader can find extensive presentations and discussions by himself over the internet. Notice that the used algorithm has been written by William L. Goffe, and can be found at http://ideas.repec.org/c/wpa/wuwppr/9406001.html.

A set of coefficients has been calculated for the silver in the same range of frequencies as in the previous sections, using two SOGP. The values of the coefficients are presented in table 4, whereas a plot of the real and imaginary parts of the permittivity function is presented in figure 2. This plot is to be compared with the Drude-Lorentz one (see figure 3), which can also be seen as a two second-order poles fitting. Although the general expression (7) is not based on a physical model, it seems to display better fittings properties than classical poles such as Lorentz ones. As in the Lorentz case, additional poles can enhance the precision of the fitting.

Summary

The physical origins of the dispersion for propagating waves have been briefly outlined and common models for EMW dispersion in metals have been presented. A generalized disper-sion models has been proposed that seems to improve the fitting of permittivity function over arbitrary sets of experimental data.

The next part of this report focuses on the Discontinuous Galerkin (DG) method in the time domain, and presents its formulation in the case of the Maxwell's equations in vacuum. Then, the respective continuous equations and Discontinuous Galerkin Time Domain (DGTD) schemes are deduced for the Drude model, and for the generalized dispersive model.

Inria

DGTD method for non-dispersive media

This section focuses on the DGTD method for Maxwell equations. Let Ω ⊂ R 3 be a bounded convex domain, and n the outward normal to its boundary ∂Ω. The electric permittivity and permeability of free space are respectively denoted ε 0 and µ 0 . Then, in free space, the Maxwell equations can be written as follows :

     µ 0 ∂H ∂t = -∇ × E, ε 0 ∂E ∂t = ∇ × H -J, (8)
along with constitutive relations B = µ 0 H and D = ε 0 E and metallic boundary conditions.

DG methods have been originally introduced in 1973 by Reed and Hill [START_REF] Reed | Triangular mesh method for the neutron transport equation[END_REF], and has been widely used since in the computational fluid dynamics field. However, their application to the time-domain Maxwell equations are more recent. DG methods can be seen as classical finite element methods (FEM) for which the global continuity of the approximation has been lifted. This implies that the support of each basis function is restrained to a discretization cell, which leads to local formulations implying no large mass matrix inversion in the process. Afterward, connexion between the cells is restored by the use of a numerical flux to evaluate the boundary integrals. The choice of the numerical flux has a great influence on the mathematical properties of the DG discretization, as energy preservation, for example. Both centered [START_REF] Fezoui | Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equation on unstructured meshes[END_REF] or upwind [START_REF] Hesthaven | Nodal high-order methods on unstructured grids[END_REF] fluxes can be used.

The discontinuity of the approximation makes room for numerous methodologic improvements, such as local approximation orders [START_REF] Fahs | Development of a hp-like discontinuous Galerkin time-domain method on non-conforming simplical meshes for electromagnetic wave propagation[END_REF], high parallelization features ([BFLP06], [START_REF] Cabel | Multi-GPU acceleration of a DGTD method for modeling human exposure to electromagnetic waves[END_REF]) and the use of non-conformeous [START_REF] Fahs | A high-order non-conforming discontinuous Galerkin method for time-domain electromagnetics[END_REF] and hybrid meshes [START_REF] Durochat | High order nonconforming multi-element discontinuous Galerkin method for time-domain electromagnetics[END_REF], for example. Also, a wide choice of time-integration schemes can be used for the discretization of time derivatives, including Leap-Frog (LF) and Runge-Kutta (RK). As well as local approximation orders can be used, local time-stepping [START_REF] Piperno | Symplectic local time-stepping in non-dissipative DGTD methods applied to wave propagation problems[END_REF] as well as locally implicit formulations [START_REF] Moya | Temporal convergence of a locally implicit discontinuous Galerkin method for Maxwell's equations[END_REF] are some of the main features that have been studied during the last years in the DG framework. More recently, Space-Time Discontinuous Galerkin methods (STDG) have emerged. These formulations exploit a temporal discretization similar to the spatial one instead of exploiting classical advancing-in-time schemes like LF or RK.

The following section presents the non-dissipative DGTD method originally developped by Fezoui et al. in [FLLP05]. Let Ω h be a discretization of Ω, relying on a quasi-uniform triangulation T h verifying T h = N i=1 T i . The internal faces of the discretization are denoted a ik = T i T k , and n ik is defined as the unit normal vector to the face a ik , oriented from T i toward T k . For each cell T i , V i is the set of indices {k | T i T k = ∅}. Then, the quasi-uniform assumption implies that :

∀T i ∈ T h , ∀k ∈ V i , ∃δ, h k ≤ δ h h i . The semi-discrete fields are denoted (H h , E h , J h ), and on each cell T i the restrictions (H i , E i , J i ) = H h Ti , E h Ti , J h Ti are defined.
For each cell, a set of scalar basis functions (φ ij ) 1 j di is defined, where d i is the number of degrees of freedom (d.o.f.) per dimension (therefore, a 3D problem implies 3 d i d.o.f. per cell).

Weak formulation

It is now possible to write the weak formulation of the problem (8) in a cell T i . By taking the dot-product of each term with a vectorial test function ψ and then integrating over the cell, one obtains the following variational problem :

Find (E, H) ∈ H 0 (rot, Ω) × H (rot, Ω) such that ∀ ψ ∈ H (rot, Ω) ,        Ti µ 0 ∂H ∂t • ψ + Ti ∇ × E • ψ = 0, Ti ε 0 ∂E ∂t • ψ - Ti ∇ × H • ψ = - Ti J • ψ,
which can be rewritten as follows using classical vectorial calculus and Green formulae :

       Ti µ 0 ∂H ∂t • ψ + Ti E • ∇ × ψ = ∂Ti (ψ × E) • n i , Ti ε 0 ∂E ∂t • ψ - Ti H • ∇ × ψ = - Ti J • ψ - ∂Ti (ψ × H) • n i .
One immediatly notices that the previous equality only holds if the boundary terms exist.

Considering the properties of the mixed product, the latter becomes :

(ψ × E) • n i = (E × n i ) • ψ,
which implies that taking E in H 0 (rot, Ω) requires the normal and tangential trace of ψ on ∂T i to exist. This implies that one should take ψ in H 1 (Ω) instead of H (rot, Ω).

Remark :

A set of di scalar basis functions φij have previously been defined on the cell Ti. Nevertheless, the test functions are to be naturally chosen as vectorial elements. Therefore, one should define three vectorial basis functions (for the three space dimensions) for each scalar one. We now denote :

φ 1 ij =   φij 0 0   , φ 2 ij =   0 φij 0   and φ 3 ij =   0 0 φij   .

Space discretization

Volumic integrals

We now seek the approximations E h and H h of E and H in the following approximation space :

V h = v ∈ L 2 (Ω) 3 , v Ti ∈ (P p (T i )) 3 ∀ T i ∈ T h , (9)
where P p (T i ) is the space of polynomials of maximum degree p on T i . The contribution of each cell is therefore defined as

E i = E h Ti
, where E i is locally expanded as :

Inria E i = 3 d=1 di j E d ij φ d ij
We now choose the test functions ψ equal to the 3

d i vectors φ d ij . The spatial discretization of Ti ε 0 ∂E ∂t • ψ then leads to 3 d i terms of the following form : Ti ε 0 T φ d ij • φ d il ∂E d ij ∂t . (10)
The previous term can be cast under a matrix form M ε0 i E i , where the mass matrix is block

diagonal of size 3 d i × 3 d i : M ε0 i =    M ε0 i 0 0 0 M ε0 i 0 0 0 M ε0 i    , with M ε0 i jl = Ti ε 0 T φ d ij • φ d il
The vector E i introduced previously has 3 d i components, and is defined as follows :

E i =    E 1 ij 1 j di E 2 ij 1 j di E 3 ij 1 j di   
From now on, the mass matrices with an exponent should be understood as :

M x i jl = Ti x T φ d ij • φ d il ,
and their associated block matrices denoted M x i .

Surface integrals

Given that E and H are discontinuous at cell boundaries, the surface integrals need a special treatment. We choose a centered approximation, which permits to rewrite ∂Ti

(ψ × E) • n as : k∈Vi a ik φ d ij × E i + E k 2 • n ik ,
with a similar choice for the surface integral involving H.

Semi-discrete equations

Taking into account the previous definitions, one obtains the following 6 d i scalar equations :

           M µ0 i ∂ H i ∂t j + Ti E i • ∇ × φ d ij = k∈Vi a ik φ d ij × E i + E k 2 • n ik , M ε0 i ∂ E i ∂t j - Ti H i • ∇ × φ d ij = - Ti J i • φ d ij - k∈Vi a ik φ d ij × H i + H k 2 • n ik . (11)
Performing an integration by parts on the surface integrals leads to :

                     M µ0 i ∂ H i ∂t j = - 1 2 Ti E i • ∇ × φ d ij + ∇ × E i • φ d ij + 1 2 k∈Vi a ik φ d ij • (E k × n ik ) , M ε0 i ∂ E i ∂t j = 1 2 Ti H i • ∇ × φ d ij + ∇ × H i • φ d ij - Ti J i • φ d ij - 1 2 k∈Vi a ik φ d ij • (H k × n ik ) . ( 12
)
Reasoning similarly as for (10) permits to rewrite (12) as :

           M µ0 i ∂ H i ∂t = -K i E i + k∈Vi S ik E k , M ε0 i ∂ E i ∂t = K i H i - k∈Vi S ik H k -M i J i ,
where the 3 d i × 3 d i local stiffness matrix and the 3 d i × 3 d k surface matrix are respectively defined by :

K i jl = 1 2 Ti φ d ij • ∇ × φ d il + φ d il • ∇ × φ d ij and S ik jl = 1 2 Ti φ d ij • φ d kl × n ik .
As mentioned previously, the basis functions are chosen in the approximation space V h .

A common choice consists in using Lagrange polynomials, though other choices are possible [CFL10]. A set of p i + 1 interpolation nodes (x j ) 0 j pi is defined in the cell, and the basis functions are then the Lagrange function L p k (x) equal to 1 on the x k node, and to 0 on all the other x j , j = k.

Time discretization

For the time derivatives, a second-order leap frog scheme (LF2) is used, where the E i are evaluated at the time station t n = n∆t, whereas the H i , J i are evaluated at the time station when knowing those of

t n+ 1 2 = n + 1 2 ∆t.
E n i , H n+ 1 2 i and J n+ 1 2 i
with the following discretization :

Inria            M µ0 i ∆t H n+ 3 2 i -H n+ 1 2 i = -K i E n+1 i + k∈Vi S ik E n+1 k , M ε0 i ∆t E n+1 i -E n i = K i H n+ 1 2 i - k∈Vi S ik H n+ 1 2 k -J n+ 1 2 i . (13)
Remark : Different time discretization can be used, like fourth-order leap frog (LF4) scheme or second and fourth-order Runge-Kutta (RK2 and RK4) schemes, under their respective stability conditions.

DGTD formulation in dispersive media

The case of Maxwell's equations in vacuum has been presented in the previous section. After a brief review of existing works, the extension of the DGTD formulation to Drude-like media is presented. Then, the same treatment is done for the generalized dispersive model.

Previous works

Before they were exploited with DGTD methods, dispersion models have been extensively used in the context of FDTD, and a very large amount of references are therefore available on this topic. We only settle here for giving a very few of them, such as [START_REF] Okoniewski | Drude dispersion in ADE FDTD revisited[END_REF] that presents an ADE-FDTD (Auxiliary Differential Equation) algorithm with the Drude model, or [START_REF] Luebbers | FDTD calculation of scattering from frequency-dependent materials[END_REF] that uses a RC-FDTD (Recursive Convolution) method for Debye, Drude and Lorentz models.

Although the exploitation of dispersion models in the DGTD framework is not as rich as the FDTD one, an important amount of studies have already been conducted. In [START_REF] Lanteri | Convergence of a discontinuous Galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF], the authors study a DGTD-CF (Centered Fluxes) method for the Maxwell equations coupled with a Debye dispersion model. Stability and convergence are proved, and a bound on the error is given. In [START_REF] Ji | High-order DGDT methods for dispersive Maxwell's equations and modelling of silver nanowire coupling[END_REF], the authors study the 2D Maxwell equations in a Drude-like medium with a DGTD-CF of the fourth order in space, and with a RK4 scheme in time.

A significant number of contributions on the numerical analysis of the schemes for dispersive media models have been made by J. Li, such as [START_REF] Li | Mathematical and numerical study of wave propagation in negative-index material[END_REF] and [START_REF] Li | Numerical convergence and physical fidelity analysis for maxwell's equations in metamaterials[END_REF] to name a few. Numerous applications of the DG method in the area photonics have been issued : [START_REF] Busch | Discontinuous Galerkin methods in nanophotonics[END_REF] presents a DGTD and DGFD (Discontinuous Galerkin Frequency Domain) method as well as realistic cases, and discuss topics like sources treatment and boundary conditions, like absorbing ones (ABC). In [START_REF] Stannigel | Discontinuous Galerkin time-domain computations of metallic nanostructures[END_REF], the authors focus on the field enhancements observed in the vicinity of metallic V-shaped nanostructures described by Drude model. A nice overview of the DGTD method coupled with Drude and Drude-Lorentz models can be found in [START_REF] Timothy | Analysis of metallic nanostructures by a discontinuous Galerkin time-domain Maxwell solver on graphics processing units[END_REF], where the author exploit a GPU implementation to focuses on the analysis of various nanostructures features.

DGTD method in Drude-like media

Maxwell-Drude equations

We now consider the case of a frequency-dependent medium, under the hypothesis of a Drude-Sommerfeld model :

ε r (ω) = ε ∞ - ω 2 d ω 2 + iωγ .
Considering a constant permeability and a homogeneous and isotropic medium, one writes the general Maxwell equations as follows :

     ∇ × E = - ∂B ∂t , ∇ × H = ∂D ∂t , (14)
along with the following constitutive relations :

RR n°8298 D = ε 0 ε ∞ E + P, B = µ 0 H. (15)
Combining ( 14) and (15) leads to :

     ∇ × E = -µ 0 ∂H ∂t , ∇ × H = ε 0 ε ∞ ∂E ∂t + ∂P ∂t .
In the frequential domain, and under the aforementioned hypothesis, the polarization P is linked to the electric field through P = -

ε0ω 2 d ω 2 +iγ d ω Ê
, where • denotes the Fourier transform of the associated field in the temporal domain. An inverse Fourier transform gives :

∂ 2 P ∂t 2 + γ d ∂P ∂t = ε 0 ω 2 d E. ( 16
)
By defining the dipolar current vector J p = ∂P ∂t , ( 14)-( 16) can be rewritten as follows :

             µ 0 ∂H ∂t = -∇ × E, ε 0 ε ∞ ∂E ∂t = ∇ × H -J p , ∂J p ∂t + γ d J p = ε 0 ω 2 d E. (17)

Normalization

The previous set of equations is now normalized. A normalized variable is denoted X when the original variable is X. First, we define vacuum impedance and light velocity :

Z 0 = µ 0 ε 0 and c 0 = 1 √ ε 0 µ 0 .
Then, the following substitutions are applied :

H = Z 0 H, E = E, J p = Z 0 J p , t = c 0 t, γ d = γ d c 0 and ω 2 d = ω 2 d c 2 0 .
One can now rewrite (17) in the form :

                 µ 0 c 0 Z 0 ∂ H ∂ t = -∇ × E, ε 0 c 0 Z 0 ε ∞ ∂ E ∂ t = ∇ × H -J p , ∂ J p ∂ t + γ d J p = Z 0 ε 0 c 2 0 c 0 ω 2 d E.
Using the equalities

µ 0 c 0 Z 0 = 1 and ε 0 c 0 Z 0 = 1, one gets : Inria                  ∂ H ∂ t = -∇ × E, ∂ E ∂ t = 1 ε ∞ ∇ × H -J p , ∂ J p ∂ t = ω 2 d E -γ d J p . (18)
It could be noticed that the first two equations are expressed in V.m -2 , whereas the third one is in V.m -3 .

Remark : To simplify the writings, the X notation for the normalized variables will be omitted from now on.

DGTD formulation for the Maxwell-Drude equations

The extension of the DGTD method to the Maxwell-Drude equations is straightforward from the classical Maxwell equations discretization (see section 2). The last equation of (18) leads to a simple vectorial equation :

1 ∆t J n+ 3 2 i -J n+ 1 2 i = -γ d J n+ 1 2 i + ω 2 d E n+1 i .
The whole system can then be written as :

                   M i ∆t H n+ 3 2 i -H n+ 1 2 i = -K i E n+1 i + k∈νi S ik E n+1 k , M ε∞ i ∆t E n+1 i -E n i = K i H n+ 1 2 i - k∈νi S ik H n+ 1 2 k -M i J n+ 1 2 i , 1 ∆t J n+ 3 2 i -J n+ 1 2 i = ω 2 d E n+1 i - γ d 2 J n+ 3 2 i + J n+ 1 2 i . ( 19
)

DGTD method for the generalized dispersive model

A similar development as the one of the previous section can be made for the general permittivity function (7).

Continuous equations

As in the Drude case, polarizations and currents are introduced to account for the dispersive behavior in the temporal domain. The (normalized) continuous equations are therefore :

                                             ∂H ∂t = -∇ × E, ∂E ∂t = 1 ε ∞ ∇ × H -J 0 + l∈L1 J l + l∈L2 J l , J 0 = σE, J l = a l E -b l P l ∀l ∈ L 1 , ∂P l ∂t = J l ∀l ∈ L 1 , ∂J l ∂t = d l ∂E ∂t + c l E -f l J l -e l P l ∀l ∈ L 2 , ∂P l ∂t = J l ∀l ∈ L 2 .
(20)

DGTD formulation

In this formulation, the currents are evaluated at the even time stations, whereas the polarizations are evaluated at the odd ones. The resulting scheme in time is, as well as in the Drude case, of order two (see section 5.2). As one might notice, some fields have been substituted to avoid unnecessary calculations :

                                                                   M i ∆t H n+ 3 2 i -H n+ 1 2 i = -K i E n+1 i + k∈νi S ik E n+1 k , M i ∆t E n+1 i -E n i = 1 ε ∞ K i H n+ 1 2 i - k∈νi S ik H n+ 1 2 k -σ M i 2 E n+1 i + E n i + M i l∈L1 b l P n+ 1 2 l,i - M i 2 l∈L2 J n+1 l,i + J n l,i
,

1 ∆t P n+ 3 2 l,i -P n+ 1 2 l,i = a l E n+1 i - b l 2 P n+ 3 2 l,i + P n+ 1 2 l,i ∀l ∈ L 1 , 1 ∆t P n+ 3 2 l,i -P n+ 1 2 l,i = J n+1 l,i ∀l ∈ L 2 , 1 ∆t J n+1 l,i -J n l,i = d l ∆t E n+1 i -E n i + c l 2 E n+1 i + E n i , - f l 2 J n+1 l,i + J n l,i -e l P n+ 1 2 l,i ∀l ∈ L 2 . ( 21
)
where σ = σ + l∈L1 a l . The resulting updating scheme is presented in figure 4.

Inria E n H n+ 1 2 J n L2 P n+ 1 2 L2 P n+ 1 2 L1 E n+1

RR n°8298 4 Theoretical study of the Maxwell-Drude equations

A stability and convergence study of the Maxwell-Drude equations is now presented. This part is highly inspired from [START_REF] Lanteri | Convergence of a discontinuous Galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF].

Stability study of the Maxwell-Drude equations

Continuous equations

First, the energy associated to the differential system ( 18) is defined at a given time t :

ξ(t) = 1 2 H(t) 2 L 2 (Ω) + ε ∞ E(t) 2 L 2 (Ω) + 1 ω 2 d J p (t) 2 L 2 (Ω) . (22)
Assuming that (E, H) have a sufficient regularity (i.e. (E, H) ∈ H 0 (rot, Ω) × H (rot, Ω) in space, and C 1 in time, for example), one would like to prove that this energy decreases with time. The L 2 scalar product on Ω of each equation leads to :

∂ξ ∂t = - Ω (∇ × E) • H + Ω (∇ × H) • E - γ d ω 2 d Ω J p • J p .
The Green formula applied to the first term implies :

Ω (∇ × E) • H = Ω (∇ × H) • E - ∂Ω (H × E) • n.
Using the property of the mixed product, and given the metallic boundary conditions, the boundary integral is equal to zero. Therefore :

∂ξ ∂t = - γ d ω 2 d J p (t) 2 L 2 (Ω) ≤ 0 (23)

Semi-discrete case

We now seek the approximated fields in the V h space that has been previously defined in equation (9). Moreover, ∀ W h ∈ V h , the following notations are used :

• The restriction W i = W h Ti ,
• The mean of W h across an interface between two cells T i and T k :

{W h } ik = W i a ik + W k a ik 2 ,
• The tangential jump of W h through an interface :

W h ik = W k a ik -W i a ik × n ik .
We focus on the differential system (18). Taking the L 2 scalar product on T i of each equation with test functions ϕ, ψ, φ sufficiently regular leads to :

Lanteri, Scheid & Viquerat                  Ti ∂H ∂t • ϕ = - Ti E • (∇ × ϕ) + ∂Ti ϕ • (E × n) , ε ∞ Ti ∂E ∂t • ψ = Ti H • (∇ × ψ) - ∂Ti (ψ × H) • n - Ti J p • ψ, 1 ω 2 d Ti ∂J p ∂t • φ = Ti E • φ - γ d ω 2 d Ti J p • φ.
One defines the semi-discrete fields (H h , E h , J h ) as solutions of the following weak formulation :

∀ (ϕ h , ψ h , φ h ) ∈ V 3 h , ∀t ∈ [0, T ], ∀i ∈ [0, N T ],                    Ti ∂H h ∂t • ϕ h = - Ti E h • (∇ × ϕ h ) + k∈Vi a ik ϕ h • ({E h } ik × n) , ε ∞ Ti ∂E h ∂t • ψ h = Ti H h • (∇ × ψ h ) - k∈Vi a ik (ψ h × {H h } ik ) • n - Ti J h • ψ h , 1 ω 2 d Ti ∂J h ∂t • φ h = Ti E h • φ h - γ d ω 2 d Ti J h • φ h . ( 24
)
and the associated semi-discrete energy ξ h is therefore defined, according to the definition of ξ, as follows :

ξ h = 1 2 H h 2 L 2 (Ω) + ε ∞ E h 2 L 2 (Ω) + 1 ω 2 d J h 2 L 2 (Ω) .
Taking (ϕ h , ψ h , φ h ) = (H h , E h , J h ), and then summing the three equations of (24) together and over all cells, one gets :

∂ξ h ∂t = Ω (H h • (∇ × E h ) -E h • (∇ × H h )) - γ d ω 2 d Ω J h • J h + Fint {H h } E h - Fint {E h } H h + ∂Ω E h • (H h × n) .
A final integration by parts 6 leads to :

∂ξ h ∂t = - γ d ω 2 d J h 2 L 2 (Ω) . (25)
Hence, ξ h (t) ≤ ξ h (0), and the stability of the semi-discrete scheme is ensured. For the rest of this study, we need to write the semi-discrete formulation over the whole domain.

Therefore, let U = (X, Y, Z) and U ′ = (X ′ , Y ′ , Z ′ ), and let the following bilinear forms :

m U, U ′ = Ω X • X ′ + ε ∞ Ω Y • Y ′ + 1 ω 2 d Ω Z • Z ′ a U, U ′ = Ω (X • (∇ × Y ′ ) -Y • (∇ × X ′ )) - Ω Z • Y ′ + Ω Y • Z ′ - γ d ω 2 d Ω Z • Z ′ b U, U ′ = Fint {X} Y ′ - Fint {Y } X ′ + ∂Ω Y ′ • (X × n) .
6 Realized on each cell T i and then summed over the whole domain.

Inria

Let U h ∈ V h be the solution of the semi-discrete problem. Given thee latter definitions, the semi-discrete formulation can be rewritten on the entire domain as :

m ∂U h ∂t , U ′ h = a U h , U ′ h + b U h , U ′ h , ∀U h ∈ V 6 h . (26)
Moreover, the solution U of the continuous equations verifies :

m ∂U ∂t , U ′ h = a U, U ′ h + b U, U ′ h , ∀U h ∈ V 6 h . (27)

Fully discrete scheme

A second-order leap-frog (LF2) time discretization is chosen. The electric field on one hand, and the magnetic field and polarization current on the other hand are evaluated respectively at the even and odd time stations. The fully discrete scheme can then be written as follows :

                                 Ti H n+ 3 2 i -H n+ 1 2 i ∆t • ϕ h = - Ti E n+1 i • (∇ × ϕ h ) + k∈Vi a ik ϕ h • E n+1 h ik × n ik , ε ∞ Ti E n+1 i -E n i ∆t • ψ h = Ti H n+ 1 2 i • (∇ × ψ h ) - k∈Vi a ik ψ h × H n+ 1 2 i ik • n ik - Ti J n+ 1 2 i • ψ h , 1 ω 2 d Ti J n+ 3 2 i -J n+ 1 2 i ∆t • φ h = Ti E n+1 i • φ h - γ d ω 2 d Ti J n+ 3 2 i + J n+ 1 2 i 2 • φ h , (28)
and we define the associated energy on the cell T i by :

ξ n i = 1 2 Ti H n+ 1 2 i • H n-1 2 i + ε ∞ Ti E n i • E n i + 1 ω 2 d Ti J n+ 1 2 i • J n-1 2 i ( 29
)
The total energy at a given time t n = n∆t is calculated as :

ξ n = N T i=0 ξ n i
Since they will be useful later, we remind to the reader the following classical inequalities [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF].

Lemma 1. let T i be a cell of the mesh T h . Then, ∃C > 0, such that ∀ u ∈ P p (T i ) :

∇ × u Ti ≤ C h u Ti , (30)
and

u 2 a ik ≤ C h u 2 Ti , (31)
RR n°8298

where • Ti and • a ik represent the L 2 norms defined respectively on the cell T i and on the face of the cell a ik .

Proposition 1 (Stability). The formulation ( 28) is stable under the following condition :

∆t < min h C , 2 ω d + γ d , 4ε ∞ C h -ω d (32)
Proof. One would like to know under which conditions the energy ξ n can be written as a definite positive form of the variables H n-1 2 , E n and J n-1 2 . To do so, we seek a lower bound of ξ n by using the equations of the system (28) to replace the occurences of H n+ 1 2 and H n-1 2 with the help of well-chosen test functions. More accurately, the following substitutions are used :

• The first equation of ( 28) is used at time t n with ϕ h = H n-1 2 i ;

• The second one at t n with ψ h = E n i ;

• The third one at t n with φ h = J n-

1 2 i .
These lead to :

ξ n i = 1 2 H n-1 2 i 2 Ti + ε ∞ E n i 2 Ti + α βω 2 d J n-1 2 i 2 Ti -∆t Ti ∇ × H n-1 2 i • E n i + ∆t β Ti E n i • J n-1 2 i + ∆t k∈Vi a ik H n-1 2 i • ({E n h } ik × n ik ) ,
where α = 1 -∆tγ d 2 and β = 1 + ∆tγ d 2 . Splitting the curl term in two parts, and integrating by parts one out of the two, one gets (the power indices are omitted up to the end of this proof) :

ξ n i = 1 2 H i 2 Ti + ε ∞ E i 2 Ti + α βω 2 d J i 2 Ti - ∆t 2 Ti (E i • (∇ × H i ) + H i • (∇ × E i )) + ∆t β Ti E i • J i + ∆t k∈Vi a ik H i • ({E h } ik × n ik ) - ∆t 2 k∈Vi a ik (H i × E i ) • n ik .
The surface integrals then rearrange as follows, given the metallic boundary conditions at the boundaries of the domain :

∆t k∈Vi a ik H i • ({E h } ik × n ik ) - ∆t 2 k∈Vi a ik (H i × E i ) • n ik = ∆t 2 k∈Fint Vi a ik H i • (E k × n ik ) - ∆t 2 k∈∂Ω Vi a ik H i • (E i × n ik ) .
Inria DGTD method in dispersive media

35

The different integrals are then bounded as follows using (30) and (31) :

Ti (E i • (∇ × H i ) + H i • (∇ × E i )) ≤ C h E i Ti H i Ti , a ik H i • (E i × n ik ) ≤ C h E i Ti H i Ti , Ti E i • J i ≤ E i Ti J i Ti .
Using the previous relations as well as the classical inequality ab ≤ 1 2 a 2 + b 2 , one gets :

ξ n i ≥ 1 2 H i 2 Ti + ε ∞ E i 2 Ti + α βω 2 d J i 2 Ti - ∆t 8
C h

E i 2 Ti + H i 2 Ti - ∆t 4 ω d β E i 2 Ti + 1 ω 2 d J i 2 Ti (33) - ∆t 8 C h a ik ∈Fint H i 2 Ti + E k 2 T k + a ik ∈∂Ω H i 2 Ti + E i 2 Ti
.

In order to make the CFL conditions more readable, it is assumed that :

∆t ≤ 2 γ d , (34)
which implies :

1 1 β 1 2 .
Therefore, adjusting C, (33) can be rewritten as :

ξ n i ≥ 1 2 1 - C∆t 4h H i 2 Ti + 1 2 ε ∞ - C∆t 4h - ∆tω d 4 E i 2 Ti - C∆t 8h k∈Vi E k 2 T k + 1 4ω 2 d 1 - γ d ∆t 2 - ω d ∆t 2 J i 2 Ti .
Then, summing over all the cells and adjusting C, one obtains an inequality involving ξ n :

ξ n ≥ 1 2 1 - C∆t h H 2 Ω + 1 2 ε ∞ - ∆t 4 C h -ω d E 2 Ω + 1 4ω 2 d 1 - ∆t 2 (γ d + ω d ) J 2 Ω .
Each one of the three induced conditions is now considered separately :

∆t < h C , ∆t < 4ε ∞ C h -ω d , ∆t < 2 ω d + γ d .
One should eventually notice that the condition (34) is contained in the last condition abovewritten, which leads to the desired result.

Proposition 2 (Bound of the discrete energy). Under CFL condition, the discrete energy ξ n can be bounded in the following way :

ξ n ≤ ξ 0 1-θ 1+θ n ∀n ∈ N *
with θ ≥ 0.

Proof. In order to bound the discrete energy, one would first be interested in knowing the sign of ξ n+1ξ n , ∀ n ≥ 0. To do so, the equations of ( 28) are used at different times and for different test functions :

ξ n+1 i -ξ n i = ∆t Ti H n+ 1 2 i • ∇ × E [n+ 1 2 ] i -E [n+ 1 2 ] i • ∇ × H n+ 1 2 i + ∆t k∈Vi a ik H n+ 1 2 i • E [n+ 1 2 ] i × n ik -E [n+ 1 2 ] i • H n+ 1 2 i × n ik + ∆t Ti J n+ 1 2 i • 1 2β E n+1 i + 1 2α E n i -E [n+ 1 2 ] i - α 2 -β 2 2αβω 2 d Ti J n+ 1 2 i • J n+ 1 2 i .
Integrating by parts the curl terms twice yields :

ξ n+1 i -ξ n i = - ∆t 2 k∈Vi a ik E [n+ 1 2 ] i × H n+ 1 2 k + E [n+ 1 2 ] k × H n+ 1 2 i • n ik + ∆t 2 γ d 4 Ti J n+ 1 2 i
• -

E n+1 i 1 + ∆tγ d 2 + E n i 1 -∆tγ d 2 η + 1 ω 2 d ∆tγ d ∆t 2 γ 2 d 4 -1 θ0 J n+ 1 2 i 2 Ti . If ∆t ≤ 2 γ d , then θ 0 < 0 and |θ| < |θ 1 | = 2 ∆t 2 γ 2 d 4 -1
. One can obtain the following bound for η :

|η| ≤ ∆t 4 |θ 1 | β J n+ 1 2 i Ti E n i Ti + α J n+ 1 2 i Ti E n+1 i Ti ≤ ∆tω d 2 |θ 1 | β J n+ 1 2 i Ti ω d max E n i Ti , E n+1 i Ti ρi ≤ ∆tω d 2 |θ 1 |    J n+ 1 2 i 2 Ti ω 2 d + ρ 2 i    . Inria where α = 1 + ∆tγ d 2 and α = 1 - ∆tγ d 2 . As specified earlier, if ∆t ≤ 2 γ d
, then one has :

β ≤ 2 and α ≤ 1.

Then, by summing the contributions over all the mesh cells, and by considering metallic boundaries, it is possible to write :

ξ n+1 -ξ n ≤ N T i=0      ∆tω d |θ 1 | 2    J n+ 1 2 i 2 Ti ω 2 d + ρ 2 i    + θ 0 ω 2 d J n+ 1 2 i 2 Ti ≤0      ≤ N T i=0 |θ 1 | ∆tω d 2    J n+ 1 2 i 2 Ti ω 2 d + ρ 2 i    .
Let C be a generic constant. Under CFL condition, the first term involving

N T i=1 J n+ 1 2 i 2
Ti is bounded by Cξ n+1 , whereas the

N T i=1
ρ 2 i one is bounded by C ξ n + ξ n+1 . This yields :

ξ n+1 ≤ |θ 1 | ∆tω d ξ n+1 + 1 + |θ 1 | ∆tω d 2 ξ n .
The latter leads to :

ξ n+1 ≤ 1 + |θ1|∆tω d 2 1 -|θ 1 | ∆tω d ξ n ≤ 1 + |θ 1 | ∆tω d 1 -|θ 1 | ∆tω d ξ n = 1 + θ 1 -θ ξ n ,
with θ ≥ 0. It is also required that θ ≤ 1, which can be shown to be equivalent to :

∆t ≤ γ 2 d 2 -2ω d + 4ω 2 d + γ 2 d .
The latter condition can be proved to be less restrictive than the ∆t ≤ 2 γ d

. Therefore, we have obtained the expected result.

Convergence of the fully discrete Maxwell-Drude DG formulation

In this section, the convergence of the Maxwell-Drude DG scheme is proven, after a preliminary convergence result about the semi-discrete formulation. We start with the following lemma. Lemma 2. Let T i be a cell of the mesh T h , and π h (•) a linear continuous projector from H s+1 (T i ) onto P p (T i ) with s ≥ 0 and k ≥ 1. Then, for u ∈ H s+1 (T i ) and m = 0, 1 : (36)

Convergence of the semi-discrete formulation

Theorem 1 (Convergence of the semi-discrete formulation). Let H, E, J p be the solution of ( 17) and

(H h , E h , J h ) ∈ C 1 [0, T ] , V 3 h the semi-discrete solution of (24). If H, E, J p ∈ C 0 [0, T ] , H s+1 (Ω)
9 for s ≥ 0, then there exists C ≥ 0 independent of h such that :

max t∈[0,T ] γ(t) 1 2 ≤ Ch min(s,p) H, E, J p C 0 ([0,T ],H s+1 (Ω) 9 ) ,
with :

γ(t) = π h (H) -H h 2 L 2 (Ω) + π h (E) -E h 2 L 2 (Ω) + π h J p -J h 2 L 2 (Ω) ∀t ∈ [0, T ] .
Proof. The orthogonal L 2 projection on the V 3 h space is defined as π h (U) = (π h (H) , π h (E) , π h J p , and the consistency error as ε

(t) = 1 2 m (π h (U) -U h , π h (U) -U h ).
Then, it is easy to bound ξ in the following way :

ε(t) ≥ 1 2 min 1, ε ∞ , 1 ω 2 d π h (H -H h ) 2 L 2 (Ω) + π h (E -E h ) 2 L 2 (Ω) + π h J p -J h 2 L 2 (Ω)
In addition, given that π h (U) -U h ∈ V h , one could easily prove that :

m π h ∂U ∂t - ∂U ∂t , π h (U) -U h = 0, (37)
a π h ∂U ∂t - ∂U ∂t , π h (U) -U h = 0. (38)
Since m, a and b are bilinear forms, and given the equalities (26 -27 -37 -38), it is possible to write (the proof is exactly the same as in [START_REF] Lanteri | Convergence of a discontinuous Galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF], pp 12 -13) :

m ∂π h (U) ∂t - ∂U h ∂t , π h (U) -U h = a (π h (U) -U h , π h (U) -U h ) + b (π h (U) -U h , π h (U) -U h ) + b (U -π h (U) , π h (U) -U h ) .
Given what has already been proven in the stability study of the semi-discrete formulation, it is straightforward to prove that we have :

a (π h (U) -U h , π h (U) -U h ) + b (π h (U) -U h , π h (U) -U h ) = - γ d ω 2 d π h J p -J h 2 Ω .
Then, the term b (Uπ h (U) , π h (U) -U h ) remains to be bounded. Let us write :

Inria b (U -π h (U) , π h (U) -U h ) = Fint {H -π h (H)} π h (E) -E h ζ - Fint {E -π h (E)} π h (H) -H h µ + ∂Ω (π h (E) -E h ) • ((H -π h (H)) × n) κ .
First, one considers the ζ term. Applying the Cauchy-Schwarz inequality leads to :

|ζ| ≤ a ik ∈Fint a ik {H -π h (H)} ik π h (E) -E h ik ≤ a ik ∈Fint a ik {H -π h (H)} ik π h (E) -E h ik ≤ a ik ∈Fint a ik |{H -π h (H)} ik | 2 1 2 a ik | π h (E) -E h ik | 2 1 2
.

The first sum can be bounded as follows :

a ik |{H -π h (H)} ik | 2 ≤ C H i -π h (H i ) 2 L 2 (a ik ) + H k -π h (H k ) 2 L 2 (a ik ) ≤ C H i -π h (H i ) 2 L 2 (∂Ti) + H k -π h (H k ) 2 L 2 (∂T k ) ≤ C h min(s,p)+ 1 2 Ti H s+1,Ti 2 + h min(s,p)+ 1 2 T k H s+1,T k 2 ≤ Cδ 2 min(s,p)+1 h 2 min(s,p)+1 Ti H 2 s+1,Ti + H 2 s+1,T k ,
the last inequality resulting from the quasi-uniform assumption made on the triangulation of the domain. Bounding the second integral can be done almost the same way, and requires the inverse inequalities presented earlier :

a ik | π h (E) -E h ik | 2 ≤ Cδ -1 h -1 Ti π h (E) -E h 2 0,Ti + π h (E) -E h 2 0,T k .
Then, by adjusting the constant C, one obtains :

|ζ| ≤ Cδ min(s,p) (i,k)|a ik ∈Fint h min(s,p) Ti H 2 s+1,Ti + H 2 s+1,T k 1 2 π h (E) -E h 2 0,Ti + π h (E) -E h 2 0,T k 1 2 ≤ Cδ min(s,p) h min(s,p) (i,k)|a ik ∈Fint H 2 s+1,Ti + H 2 s+1,T k 1 2 π h (E) -E h 2 0,Ti + π h (E) -E h 2 0,T k 1 2 , RR n°8298 40

Lanteri, Scheid & Viquerat

where h = max Ti h Ti . Eventually, one must consider the sum with respect to all the internal faces, which can be rewritten as follows by adjusting once again the generic constant C :

|ζ| ≤ Cδ min(s,p) h min(s,p) Ti H s+1,Ti π h (E) -E h 0,Ti .

Bounding µ can be done the exact same way :

|µ| ≤ Cδ min(s,p) h min(s,p) Ti E s+1,Ti π h (H) -H h 0,Ti .

Moreover, one could easily prove that κ = 0, given the metallic boundary condition on ∂Ω.

Therefore :

|ζ| + |µ| ≤ Ti Cδ min(s,p) h min(s,p) Ti H s+1,Ti π h (E) -E h 0,Ti + E s+1,Ti π h (H) -H h 0,Ti ≤ C max Ti δ min(s,p) h min(s,p) Ti   Ti H 2 s+1,Ti 1 2 Ti π h (E) -E h 2 0,Ti 1 2 + Ti E 2 s+1,Ti 1 2 Ti π h (H) -H h 2 0,Ti 1 2   ≤ Cδ min(s,p) h min(s,p) H s+1,Ω π h (E) -E h 0,Ω + E s+1,Ω π h (H) -H h 0,Ω ≤ Cδ min(s,p) h min(s,p) (H, E) s+1,Ω π h (E) -E h 2 0,Ω + π h (H) -H h 2 0,Ω 1 2 . Since ε(t) = 1 2 t 0 m ∂ (π h (U) -U h )
∂s , π h (U) -U h ds, and assuming that ε(0) = 0, one can write :

ε(t) = 1 2 t 0 - γ d ω 2 d J h 2 Ω ds + 1 2 t 0 b (U -π h (U) , π h (U) -U h ) ds ≤ 1 2 t 0 b (U -π h (U) , π h (U) -U h ) ds.
In regard of the previous results, one obtains :

ε(t) ≤ Ch min(s,p) t 0 (H, E) s+1,Ω π h (E) -E h 2 0,Ω + π h (H) -H h 2 0,Ω 1 2 ds ≤ Ch min(s,p) t 0 H, E, J p s+1,Ω γ(s) 1 2 ds, therefore : γ(t) ≤ Ch min(s,p) υ H, E, J p C 0 ([0,T ],H s+1 (Ω) 3 ) , where υ = max [0,T ] γ(t) 1 2 .
Taking the maximum value of the left hand side over [0, T ] permits to conclude. Inria

Convergence of the fully discrete formulation

We consider the fully discrete formulation presented previously in equation ( 28). We choose a constant time step ∆t verifying the CFL condition (32) and such that N ∆t = T , where N is the total number of timesteps. Two classical results about Taylor-Lagrange expansions will be required :

Lemma 3. Let be U ∈ C 3 ([t n , t n+1 ]), then ∃ (c n , c n+1 ) ∈]t n , t n+ 1 2 [×]t n+ 1 2 , t n+1 [ and (d n , d n+1 ) ∈ ]t n , t n+ 1 2 [×]t n+ 1 2 , t n+1 [ such that : U h (t n+1 ) -U h (t n ) = ∆t ∂U h ∂t t n+ 1 2 + ∆t 3 28 ∂ 3 U h ∂t 3 (c n+1 ) + ∂ 3 U h ∂t 3 (c n ) , and
1 2 (U h (t n+1 ) + U h (t n )) = U h t n+ 1 2 + ∆t 2 16 ∂ 2 U h ∂t 2 (d n+1 ) + ∂ 2 U h ∂t 2 (d n ) .
Theorem 2 (Convergence of the fully discrete formulation). Let be :

H, E, J p ∈ C 3 [0, T ] , L 2 (Ω) 9 C 0 [0, T ] , H s+1 (Ω) 9 .
Under a CFL condition as in (32), the following error estimate holds :

max n∈[0,N ] H t n+ 1 2 -H n+ 1 2 h 2 L 2 (Ω) 3 + E (t n ) -E n h 2 L 2 (Ω) 3 + J p t n+ 1 2 -J n+ 1 2 h 2 L 2 (Ω) 3 1 2 ≤ C ∆t 2 + h min(s,k) H, E, J p C 3 ([0,T ],L 2 (Ω) 9 ) + H, E, J p C 0 ([0,T ],H s+1 (Ω) 9 )
Proof. We define the consistency error as follows :

ε n+1 h = E h (t n+1 ) -E n+1 h 2 L 2 (Ω) + H h t n+ 3 2 -H n+ 3 2 h 2 L 2 (Ω) + J h t n+ 3 2 -J n+ 3 2 h 2 L 2 (Ω) 1 2
,

where E n+1 h , H n+ 3 2 h and J n+ 3 2 h
are defined as :

                                     Ti H n+ 3 2 h -H h t n+ 1 2 ∆t • ϕ h = - Ti E h (t n+1 ) • (∇ × ϕ h ) + k∈Vi a ik ϕ h • ({E h (t n+1 )} ik × n ik ) , ε ∞ Ti E n+1 h -E h (t n ) ∆t • ψ h = Ti H h t n+ 1 2 • (∇ × ψ h ) - Ti J h t n+ 1 2 • ψ h - k∈Vi a ik ψ h × H h t n+ 1 2 ik • n ik , 1 ω 2 d Ti J n+ 3 2 h -J h t n+ 1 2 ∆t • φ h = Ti E h (t n+1 ) • φ h - γ d ω 2 d Ti J h t n+ 3 2 + J h t n+ 1 2 2 • φ h .

RR n°8298

The following developments mainly focus on the J equation, since the treatment of the equation for the update of E and H is very similar to the one done in [START_REF] Lanteri | Convergence of a discontinuous Galerkin scheme for the mixed time domain Maxwell's equations in dispersive media[END_REF]. One could rewrite the last equation as follows :

1 ω 2 d ∆t Ti J n+ 3 2 h -J h t n+ 3 2 + J h t n+ 3 2 -J h t n+ 1 2 • φ h = Ti E h (t n+1 ) • φ h - γ d ω 2 d Ti J h t n+ 3 2 + J h t n+ 1 2 2 • φ h .
The semi-discrete equation involving J h used at time t = t n+1 reads as :

1 ω 2 d Ti ∂J h ∂t (t n+1 ) • φ h = Ti E h (t n+1 ) • φ h - γ d ω 2 d Ti J h (t n+1 ) • φ h .
Then, substituting the E h term in the fully discrete equation leads to :

1 ω 2 d ∆t Ti J n+ 3 2 h -J h t n+ 3 2 + J h t n+ 3 2 -J h t n+ 1 2 • φ h = 1 ω 2 d Ti ∂J h ∂t (t n+1 ) • φ h + γ d ω 2 d Ti J h (t n+1 ) • φ h - γ d ω 2 d Ti J h t n+ 3 2 + J h t n+ 1 2 2 • φ h .
Exploiting the Taylor-Lagrange equalities and notations from lemma 3 gives :

1 ∆t Ti J n+ 3 2 i -J i t n+ 3 2 • φ h + ∆t 2 28 Ti ∂ 3 J h ∂t 3 (c n+1 ) + ∂ 3 J h ∂t 3 (c n ) + ∆t 2 16 Ti ∂ 2 J h ∂t 2 (d n+1 ) + ∂ 2 J h ∂t 2 (d n ) = 0, therefore : J n+ 3 2 h -J h t n+ 3 2 L 2 (Ω) ≤ C∆t 3 ∂ 2 J h ∂t 2 (t) L 2 (Ω) + ∂ 3 J h ∂t 3 (t) L 2 (Ω)
.

Following the same idea for the other two fields allows to write :

|ε n h | ≤ C∆t 3 U C 3 ([0,T ],L 2 (Ω) 9 ) . ( 39
) One could now define ĵn i (φ h ) = 1 ω 2 d ∆t Ti J h t n+ 3 2 -J n+ 3 2 h
• φ h as well as ên i and ĥn i associated to their related fields. One has :

1 ω 2 d ∆t Ti J h t n+ 3 2 -J h t n+ 1 2 • φ h + γ d ω 2 d Ti J h t n+ 3 2 + J h t n+ 1 2 2 • φ h = Ti E h (t n+1 ) • φ h + ĵn i (φ h ). Inria We define ĵn h (φ h ) = N T i=0 ĵn i (φ h ) = 1 ω 2 d ∆t Ω J h t n+ 3 2 -J n+ 3 2 h
• φ h , which can be easily bounded as follows :

ĵn h (φ h ) ≤ C ∆t J n+ 3 2 h -J h t n+ 3 2 L 2 (Ω) φ h L 2 (Ω) .
Following the same development for ên h and ĥn h , one eventually proves that, • being the operator norm associated to the considered linear forms on L 2 (Ω) 3 :

ĵn h + ĥn h + ên h ≤ C∆t 2 U C 3 ([0,T ],L 2 (Ω) 9 ) . ( 40
)
One now focuses on the error energy. One defines J

n+ 1 2 h = J h t n+ 1 2 -J n+ 1 2 h
, as well as their

E and H counterparts H n+ 1 2 h
and E n h . The error energy then writes as :

εn i = 1 2 ε ∞ E n i 2 L 2 (Ω) + H n-1 2 i , H n+ 1 2 i Ti + 1 ω 2 d J n-1 2 i , J n+ 1 2 i Ti , where J n+ 1 2 i verifies : 1 ω 2 d ∆t Ti J n+ 3 2 i -J n+ 1 2 i • φ h = Ti E n+1 i • φ h - γ d ω 2 d Ti J n+ 3 2 i + J n+ 1 2 i 2 • φ h + ĵn i .
Combinations similar to the ones in the stability study prove that, under a CFL-like condition :

εn+1 -εn = i∈[0,N T ] ∆t 2 γ d 4 Ti J n+ 1 2 i • - E n+1 i 1 + ∆tγ d 2 + E n i 1 -∆tγ d 2 + 1 ω 2 d ∆tγ d ∆t 2 γ 2 d 4 -1 J n+ 1 2 i 2 Ti + ên h E [n+ 1 2 ] h + 1 2 ĥn h H n+ 1 2 h + ĥn+1 h H n+ 1 2 h + 1 2 ĵn h J n+ 1 2 h + ĵn+1 h J n+ 1 2 h .
Following the same ideas, one obtains :

H n h 2 L 2 (Ω) + E n h 2 L 2 (Ω) + J n h 2 L 2 (Ω) 1 2 ≤ C∆t 2 U C 3 ([0,T ],L 2 (Ω) 9 ) . (41)
It only remains to put together the theorem (3) along with (47) and the lemma 2 :

U (t n ) -U n h L 2 (Ω) 9 ≤ Ch min(s,k) U C 0 ([0,T ],L 2 (Ω) 9 ) + Ch min(s,k) U C 0 ([0,T ],H s+1 (Ω) 9 ) + C∆t 2 U C 3 ([0,T ],L 2 (Ω) 9 ) .
Taking the maximum of the latter over all n ∈ [0, N ] then leads to the desired result.

RR n°8298

5 Theoretical study of the generalized dispersive model

Stability study of the generalized dispersive formulation

Continuous equations

Let N be the total number of fields involved :

N = 2 + card (L 1 ) + 2card (L 2 ) .
We define the energy associated to the system (20) at a given time t as follows :

ξ(t) = 1 2 H(t) 2 L 2 (Ω) + ε ∞ E(t) 2 L 2 (Ω) + l∈L1 b l a l P l (t) 2 L 2 (Ω) + l∈L2 e l c l + d l f l P l (t) 2 L 2 (Ω) + l∈L2 1 c l + d l f l J l (t) -d l E(t) 2 L 2 (Ω) .
Assuming a sufficient regularity for each of the fields, using equation ( 20) leads to :

∂ξ ∂t = -σ + l∈L1 a l + l∈L2 c l d l c l + d l f l E 2 L 2 (Ω) - l∈L2 f l c l + d l f l J l 2 L 2 (Ω) + l∈L2 e l d l c l + d l f l E • P l - l∈L1 b 2 l a l P 2 L 2 (Ω) + l∈L1 2b l E • P l .
The L 1 terms can be rewritten the following way :

-a l E 2 L 2 (Ω) - b 2 l a l P l 2 L 2 (Ω) + 2b l E • P l = - b 2 l a l a l b l E -P l 2 L 2 (Ω)
.

Then, every term is negative, except for the E • P l one. The latter can be bounded in the following way :

l∈L2 e l d l c l + d l f l E • P l ≤ αξ, where α = 1 √ ε ∞ max L2 d l e l c l + d l f l ξ.
One might notice that d l = 0 or e l = 0 implies a decreasing energy. Although it is not provided in the current study, a physical interpretation of the d l parameter might be of interest for future fitting improvements, since it appears to be responsible of the non-decreasing energy. Then : ∂ξ ∂t αξ which directly implies that, ∀t ∈ [0, T ] : ξ(t) ξ(0) exp(αt).

Semi-discrete case

From now on, subscript indices l will be omitted in the sums over L 1 or L 2 sets, and a vector field written A Lu,h must be understood as A l∈Lu,h . We define the semi-discrete fields H h , E h , P L1,h , P L2,h , J L2,h as solutions of the following weak formulation :

∀ (ϕ h , ψ h , φ h , κ h , Π h ) ∈ V 3 h , ∀t ∈ [0, T ], ∀i ∈ [0, N T ],                                                            Ti ∂H h ∂t • ϕ h = - Ti E h • (∇ × ϕ h ) + k∈Vi a ik ϕ h • ({E h } ik × n) , ε ∞ Ti ∂E h ∂t • ψ h = Ti H h • (∇ × ψ h ) - k∈Vi a ik (ψ h × {H h } ik ) • n - Ti σE h • ψ h - L1 Ti (a l E h -P l,h ) • ψ h - L2 Ti J l,h • ψ h , Ti ∂P l,h ∂t • φ h = Ti a l E h • φ h - Ti b l P l,h • φ h l ∈ L 1 , Ti ∂P l,h ∂t • κ h = Ti J l,h • κ h l ∈ L 2 , Ti ∂J l,h ∂t -d l ∂E h ∂t • Π h = Ti c l E h • Π h - Ti f l J l,h • Π h - Ti e l P l,h • Π h l ∈ L 2 . (42)
Following what has been done in the Drude case for the treatment of the boundaries, as well as the ideas of the continuous equations, one easily obtains : ξ h (t) ξ h (0) exp(αt).

Fully discrete scheme

The time discretization of the considered system can be written as follows :

Inria                                                                              Ti H n+ 3 2 i -H n+ 1 2 i ∆t • ϕ h = - Ti E n+1 i • (∇ × ϕ h ) + k∈Vi a ik ϕ h • E n+1 h ik × n ik , ε ∞ Ti E n+1 i -E n i ∆t • ψ h = Ti H n+ 1 2 i • (∇ × ψ h ) - k∈Vi a ik ψ h × H n+ 1 2 i ik • n ik -σ + l∈L1 a l Ti E [n+ 1 2 ] i • ψ h + l∈L1 Ti P n+ 1 2 l,i • ψ h - l∈L2 Ti J [n+ 1 2 ] l,i • ψ h , Ti P n+ 3 2 l,i -P n+ 1 2 l,i ∆t • φ h = Ti E n+1 i • φ h -b l Ti P n+ 3 2 l,i + P n+ 1 2 l,i 2 • φ h l ∈ L 1 , Ti P n+ 3 2 l,i -P n+ 1 2 l,i ∆t • κ h = Ti J n+1 l,i • κ h l ∈ L 2 , Ti J l,i -d l E i n+1 -J l,i -d l E i n ∆t • Π h = Ti E [n+ 1 2 ] i • Π h -f l Ti J [n+ 1 2 ] l,i • Π h -e l Ti P n+ 1 2 l,i • Π h l ∈ L 2 , (43)
and its associated energy for the cell T i is :

ξ n i = 1 2 Ti H n+ 1 2 i • H n-1 2 i + ε ∞ Ti E n i • E n i + l∈L1 a l b l Ti P n+ 1 2 l,i • P n-1 2 l,i + l∈L2 1 c l + d l f l Ti J l,i -d l E i n • J l,i -d l E i n + l∈L2 e l c l + d l f l Ti P n+ 1 2 l,i • P n-1 2 l,i
. Proposition 3 (Stability). The formulation (43) is stable under the following condition :

∆t ≤ min   h C , min   2 e l , 1 d l , 2 a l + b l , ε ∞ C h + 1 2 L1 a 2 l 2b l + 1 2 L2 e l d l c l +d l f l     . (44)
Proof. An adequate choice of test functions for the scheme (43) similar to what was done in the Drude case allows us to write :

Ti P n+ 1 2 l,i • P n-1 2 l,i = 1 1 + l (1 -l ) P n-1 2 l,i 2 Ti + a l ∆t Ti E n i • P n-1 2 l,i
, for l ∈ L 1 , with the definition l = b l ∆t 2 . For l ∈ L 2 , one has :

Ti P n+ 1 2 l,i • P n-1 2 l,i = P n-1 2 l,i 2 Ti + ∆t Ti J n l,i • P n-1 2 l,i .

RR n°8298

The treatment of the H term is exactly identical to what has been donne in the Drude section. Therefore :

ξ n i = 1 2 H i 2 Ti + ε ∞ E i 2 Ti + l∈L1 a l (1 -l ) b l (1 + l ) P l,i 2 Ti + l∈L2 e l c l + d l f l P l,i 2 Ti + l∈L2 1 c l + d l f l J l,i -d l E i 2 Ti + ∆t k∈Vi a ik H i • ({E h } ik × n ik ) - ∆t 2 k∈Vi a ik (H i × E i ) • n ik + l∈L2 e l c l + d l f l ∆t Ti J l,i • P l,i + l∈L1 a l b l (1 + l ) a l ∆t Ti E i • P l,i . If one assumes that ∆t ≤ 2 b l ∀l ∈ L 1 , then 1 1+ l ≥ 1 2 . The minoration of ξ n
i is similar to the Drude case, and one eventually obtains :

ξ n ≥ 1 2 1 - C∆t h H 2 Ω + 1 2 ε ∞ - C∆t h - 1 2 L1 a 2 l ∆t 2b l - 1 2 L2 e l d l ∆t c l + d l f l E 2 Ω + 1 2 L2 1 c l + d l f l 1 - e l ∆t 2 2 J l -d l E 2 Ω + 1 2 L2 e l c l + d l f l (1 -∆td l ) P l 2 Ω + 1 2 L1 a l 2b l 1 -l - a l ∆t 2 P l 2 Ω .
which leads to the desired result.

Proposition 4 (Bound on the discrete energy). Under CFL condition, the discrete energy ξ n can be bounded in the following way :

ξ n ≤ ξ 0 1-θ 1+2θ n ,
with θ ≥ 0.

Proof. Following the same ideas as for the Drude model, one obtains :

Inria ξ n+1 i -ξ n i = - ∆t 2 k∈Vi a ik E [n+ 1 2 ] i • H n+ 1 2 k + E [n+ 1 2 ] k • H n+ 1 2 i -∆t σ + L1 a l E [n+ 1 2 ] 2 Ti + ∆t L1 b l + a 2 l b l (1 + b l ∆t) Ti P n+ 1 2 l,i • E [n+ 1 2 ] i -∆t L2 Ti J [n+ 1 2 ] l,i • E [n+ 1 2 ] i + ∆t L2 e l c l + d l f l Ti J [n+ 1 2 ] l,i • P n+ 1 2 l,i + ∆t L2 1 c l + d l f l Ti c l E [n+ 1 2 ] i -f l J [n+ 1 2 ] l,i -e l P n+ 1 2 l,i • J l,i -d l E i [n+ 1 2 ] , which simplifies in ξ n+1 i -ξ n i = - ∆t 2 k∈Vi a ik E [n+ 1 2 ] i • H n+ 1 2 k + E [n+ 1 2 ] k • H n+ 1 2 i -∆t σ + L1 a l + L2 c l d l c l + d l f l E [n+ 1 2 ] i 2 Ti -∆t L2 f l c l + d l f l J [n+ 1 2 ] l,i 2 Ti + ∆t L1 b l + a 2 l b l (1 + b l ∆t) Ti P n+ 1 2 l,i • E [n+ 1 2 ] i + ∆t L2 d l e l c l + d l f l Ti P n+ 1 2 l,i • E [n+ 1 2 ] i .
Summing over all the cells and then bounding the scalar products the same way as before gives (remind that 1 1+b l ∆t ≤ 1), under CFL condition :

ξ n+1 i -ξ n i ≤ ∆t 2 L1 b l 1 + a 2 l b 2 l P n+ 1 2 l 2 Ω + E [n+ 1 2 ] 2 Ω + ∆t 2 L2 d l e l c l + d l f l P n+ 1 2 l 2 Ω + E [n+ 1 2 ] 2 Ω ≤ θ 2ξ n + ξ n+1 , with θ = ∆t 2 L1 b l 1 + a 2 l b 2 l + L2 d l e l c l + d l f l
. The last inequality then leads to the result.

Convergence of the fully discrete generalized dispersive DG formulation

Convergence of the semi-discrete formulation

As in the Drude case, we start by defining the following bilinear forms. Let U = (V, W, X, Y, Z)

and

U ′ = (V ′ , W ′ , X ′ , Y ′ , Z ′ ), then : m U, U ′ = Ω V • V ′ + ε ∞ Ω W • W ′ + L1 b l a l Ω X • X ′ + L2 1 c l + d l f l Ω (Y -d l W ) • (Y ′ -d l W ′ ) + L2 1 c l + d l f l Ω Z • Z ′ , a U, U ′ = Ω (V • (∇ × W ′ ) -W • (∇ × V ′ )) - Ω σW • W ′ - L1 Ω (a l W • W ′ -b l X • W ′ ) + L1 Ω (a l W • X ′ -b l X • X ′ ) - L2 Ω Y • W ′ + L2 Ω Y • Z ′ + L2 Ω (c l W -f l Y -e l Z) • (Y ′ -d l W ′ ) , b U, U ′ = Fint {V } W ′ - Fint {W } V ′ + ∂Ω W ′ • (V × n) .
The semi-discrete scheme over the entire domain can be written as :

m ∂U h ∂t , U ′ h = a U h , U ′ h + b U h , U ′ h .
Moreover, the solution U of the continuous equations being solution of the semi-discrete scheme, one has :

m ∂U ∂t , U ′ h = a U, U ′ h + b U, U ′ h .
Theorem 3 (Convergence of the semi-discrete formulation). Let :

H, E, (P l ) l∈L1 , (P l ) l∈L2 , (J l -d l E) l∈L2
be the solution of (20), and :

H h , E h , (P h,l ) l∈L1 , (P h,l ) l∈L2 , J h,l -d l E l∈L2 ∈ C 1 [0, T ] , V N h
the semi-discrete solution of the associated semi-discrete formulation. If :

H, E, (P l ) l∈L1 , (P l ) l∈L2 , (J l -d l E) l∈L2 ∈ C 0 [0, T ] , H s+1 (Ω) 3N
for s ≥ 0, then there exists C ≥ 0 independent of h such that :

γ(t) 1 2 ≤ Ch min(s,p) H, E, (P l ) l∈L1 , (P l ) l∈L2 , (J l -d l E) l∈L2 C 0 ([0,T ],H s+1 (Ω) 3N ) e αt ,

Inria

where α has been defined in section 5.1.1, and :

γ(t) = π h (H) -H h 2 L 2 (Ω) + π h (E) -E h 2 L 2 (Ω) + l∈L1 π h (P l ) l∈L1 -(P h,l ) l∈L1 2 L 2 (Ω) + l∈L2 π h (P l ) l∈L2 -(P h,l ) l∈L2 2 L 2 (Ω) + l∈L2 π h (J l -d l E) l∈L2 -J h,l -d l E h l∈L2 2 L 2 (Ω)
.

Proof. Like previously, we have :

ε(t) = 1 2 t 0 m ∂ (π h (U) -U h ) ∂s , π h (U) -U h ds,
where U is the solution of the continuous problem. Following what has been done in the Drude case, we have :

m ∂π h (U) ∂t - ∂U h ∂t , π h (U) -U h = a (π h (U) -U h , π h (U) -U h ) + b (π h (U) -U h , π h (U) -U h ) + b (U -π h (U) , π h (U) -U h ) .
As highlighted in section 5.1.1, the d l parameter prevents the energy to be decreasing.

Indeed, unlike the Drude case, one has :

a (π h (U) -U h , π h (U) -U h ) + b (π h (U) -U h , π h (U) -U h ) = L2 Ω e l d l (π h (Z) -Z h ) • (π h (W ) -W h ) ,
which can be bounded by a Cauchy-Schwarz inequality :

a (π h (U) -U h , π h (U) -U h ) + b (π h (U) -U h , π h (U) -U h ) ≤ αε,
where α has been defined in section 5.1.1. The bounding of the term b (U -

π h (U) , π h (U) -U h )
is strictly identical to what has been done for the Drude convergence proof. Therefore :

ε ′ (t) ≤ αε(t) + Ch min(s,p) (H, E) s+1,Ω π h (E) -E h 2 0,Ω + π h (H) -H h 2 0,Ω 1 2
≤ αε(t) + Ch min(s,p) H, E, (P l ) l∈L1 , (P l ) l∈L2 , (J ld l E) l∈L2 s+1,Ω γ(t)

1 2
≤ αε(t) + Ch min(s,p) H, E, (P l ) l∈L1 , (P l ) l∈L2 , (J l -

d l E) l∈L2 s+1,Ω ε(t) 1 2 .
Dividing by ε(t) 1 2 , one obtains :

ε ′ (t) ε(t) 1 2
≤ αε(t)

1 2 + Ch min(s,p) H, E, (P l ) l∈L1 , (P l ) l∈L2 , (J ld l E) l∈L2 s+1,Ω .

RR n°8298

Then, a Grönwall inequality gives, assuming that ε(0) = 0 : ε(t) 1 2 ≤ 1 α Ch min(s,p) H, E, (P l ) l∈L1 , (P l ) l∈L2 , (J ld l E) l∈L2 s+1,Ω e αt -1

Therefore :

γ(t) 1 2 ≤ Ch min(s,p) δ H, E, (P l ) l∈L1 , (P l ) l∈L2 , (J l -d l E) l∈L2 C 0 ([0,T ],H s+1 (Ω) 3N ) e αt ,
which is the expected result.

Convergence of the fully discrete formulation

Theorem 4 (Convergence of the fully discrete formulation). Let be :

H, E, (P l ) l∈L1 , (P l ) l∈L2 , (J l -d l E) l∈L2 ∈ C 3 [0, T ] , L 2 (Ω) 3N C 0 [0, T ] , H s+1 (Ω) 3N .
Under the CFL condition (44), the following error estimate holds :

max n∈[0,N ] H t n+ 1 2 -H n+ 1 2 h 2 L 2 (Ω) 3 + E (t n ) -E n h 2 L 2 (Ω) 3 + L1 P l∈L1 t n+ 1 2 -P n+ 1 2 h,l∈L1 2
L 2 (Ω) 3 + L2 P l∈L2 t n+ 1 2 -P n+ 1 2 h,l∈L2 2
L 2 (Ω) 3 + L2 (J l -d l E) l∈L2 (t n ) -J h,l -d l E h n l∈L2 2 L 2 (Ω) 3 1 2 ≤ C ∆t 2 + h min(s,p)
H, E, (P l ) l∈L1 , (P l ) l∈L2 , (J ld l E) l∈L2 C 3 ([0,T ],L 2 (Ω) 3N ) + H, E, (P l ) l∈L1 , (P l ) l∈L2 , (J ld l E) l∈L2 C 0 ([0,T ],H s+1 (Ω) 3N ) .

Proof. The consistency error is defined as follows :

ε n+1 h = E h (t n+1 ) -E n+1 h 2 L 2 (Ω) + H h t n+ 3 2 -H n+ 3 2 h 2 L 2 (Ω) + L1 P h,l∈L1 t n+ 3 2 -P n+ 3 2 h,l∈L1 2
L 2 (Ω) + L2 P h,l∈L2 t n+ 3 2 -P n+ 3 2 h,l∈L2 2
L 2 (Ω) + L2 J h,l -d l E h l∈L2 (t n+1 ) -J h,l -d l E h n+1 l∈L2 2 L 2 (Ω) 1 2
,

where E n+1 h , H n+ 3 2 h , P n+ 3 2 h,l∈L1 , P n+ 3 2 h,l∈L2 and J h,l -d l E h n+1 l∈L2
are defined as :

Inria                                                                                                Ti H n+ 3 2 i -H i t n+ 1 2 ∆t • ϕ h = - Ti E i (t n+1 ) • (∇ × ϕ h ) + k∈Vi a ik ϕ h • ({E h (t n+1 )} ik × n ik ) , ε ∞ Ti E n+1 i -E i (t n ) ∆t • ψ h = Ti H i t n+ 1 2 • (∇ × ψ h ) - k∈Vi a ik ψ h × H i t n+ 1 2 ik • n ik -σ + L1 a l Ti E i (t n+1 ) + E i (t n ) 2 • ψ h + L1 Ti P l,i t n+ 1 2 • ψ h - L2 Ti J l,i (t n+1 ) + J l,i (t n ) 2 • ψ h , Ti P n+ 3 2 l,i -P l,i t n+ 1 2 ∆t • φ h = Ti E i (t n+1 ) • φ h -b l Ti P l,i t n+ 3 2 + P l,i t n+ 1 2 2 • φ h l ∈ L 1 , Ti P n+ 3 2 l,i -P l,i t n+ 1 2 ∆t • κ h = Ti J l,i (t n+1 ) • κ h l ∈ L 2 , Ti J l,i -d l E i n+1 -J l,i -d l E i (t n ) ∆t • Π h = Ti E i (t n+1 ) + E i (t n ) 2 • Π h -f l Ti J l,i (t n+1 ) + J l,i (t n ) 2 • Π h -e l Ti P l,i t n+ 1 2 l ∈ L 2 .
Following the same procedure as for the Drude case, one can prove that :

|ε n h | ≤ C∆t 3 U C 3 ([0,T ],L 2 (Ω) 15 ) , (45)
which implies :

ĥn h + ên h + pn h,l∈L1 + pn h,l∈L2 + ĵn h,l∈L2 -d l ên h ≤ C∆t 2 U C 3 ([0,T ],L 2 (Ω) 15 ) . ( 46
)
where the natural definitions of pn h,l∈L1 , pn h,l∈L2 and ĵn h,l∈L2d l ên h are coherent with the one given in the Drude case. The error energy is defined as :

εn i = 1 2 ε ∞ E n i 2 L 2 (Ω) + H n-1 2 i , H n+ 1 2 i Ti + L1 a l b l P n-1 2 l,i , P n+ 1 2 l,i Ti + L2 e l c l + d l f l P n-1 2 l,i , P n+ 1 2 l,i Ti + L2 1 c l + d l f l J n l,i -d l E n i 2 L 2 (Ω) ,
where the notations U n i has been defined to the Drude case. Under CFL condition, one readily gets :

εn+1 -εn = i∈[0,N T ] ∆t L1 b l + a 2 l b l (1 + b l ∆t) Ti P n+ 1 2 l,i • E [n+ 1 2 ] i + L2 d l e l c l + d l f l Ti P n+ 1 2 l,i • E [n+ 1 2 ] i -σ + L1 a l + L2 c l d l c l + d l f l E [n+ 1 2 ] i 2 Ti - L2 f l c l + d l f l J [n+ 1 2 ] l,i 2 Ti + ên h E [n+ 1 2 ] h + 1 2 ĥn h H n+ 1 2 h + ĥn+1 h H n+ 1 2 h + 1 2 L1 pn h,l P n+ 1 2 h,l + pn+1 h,l P n+ 1 2 h,l + 1 2 L2 pn h,l P n+ 1 2 h,l + pn+1 h,l P n+ 1 2 h,l + 1 2 L2 ĵn h,l -d l ên h (J l,i -d l E i ) [n+ 1 2 ] .
Considering inequality (46) as well as classical bounding methods, the following result is obtained :

H n h 2 L 2 (Ω) + E n h 2 L 2 (Ω) + P n h,l∈L1 2
L 2 (Ω) + P n h,l∈L2 2
L 2 (Ω) + (J h,l -d l E h ) n l∈L2 2 L 2 (Ω) 1 2 ≤ C∆t 2 U C 3 ([0,T ],L 2 (Ω) 3N ) . ( 47
)
The end of the proof is then similar to the Drude case.

Inria

Numerical results for the Drude model

This section presents some simulation results for different situations : first, an artificiallybuilt case is considered to validate the implementation of the method. Then, a physical situation is considered, and comparisons are made against theoretical results.

Artificial validation case

Presentation

A cubic cavity of length L with metallic walls is considered. The goal is to simulate a standing wave in the cavity, given a volumic source. To do so, the analytical solution is built from the classical cavity modes. For the E and H fields, the following expressions are looked for

E(x, t) =   -f (x) 0 g(x)   cos (x k t) , H(x, t) =   h(x) i(x) j(x)   sin (x k t) ,
where x m = n k π L , n k representing the number of the mode. We used : Given that the previous modes cannot be solutions of (18), a volumic source is added in the second equation :

               f (x,
∂E ∂t = 1 ε ∞ ∇ × H -J p + J s , (48)
where the source current J s is defined as :

J s (x, t) =   (β sin(x k t) -γ d cos(x k t)) f (x) 0 (-β sin(x k t) + γ d cos(x k t)) g(x)   , (49)
assuming that J p (x, y, z, 0) is taken equal to zero. The parameters used in the previous expression are defined as follows :

α = ω 2 d x 2 k + γ 2 d , β = ε ∞ x k -3 x 2 m x k -αx k .
Being given these expressions, the analytical solution for the polarization current is :

J p (x, t) =   -α (x k sin(x k t) + γ d cos(x k t)) f (x) 0 α (x k sin(x k t) + γ d cos(x k t)) g(x)   . (50)
The chosen Drude parameters for this case are the ones given in section 1.2.1, table 1. The (normalized) maximal time t max , (here 2.0 10 -6 ), and the cavity length are chosen in order to simulate two temporal periods. The spatial period is related to the parameter n k , that describes the spatial shape of the mode. Eventually, the frequency f is taken equal to 3.0 10 5 GHz.

Results

We denote E, H, J p the exact solution of (48), and (E h , H h , J h ) the calculated solution. The results obtained from the DG calculation are coherent and close to the analytical ones : one can refer to figure 5 for some exact and calculated field plots. At a given time t n , the total L 2 error is :

e n h = H t n+ 1 2 -H n+ 1 2 h 2 L 2 (Ω) + ε ∞ E (t n ) -E n h 2 L 2 (Ω) (51) + 1 ω 2 d J p t n+ 1 2 -J n+ 1 2 h 2 L 2 (Ω) 1 2 (52)
For a given time step ∆t verifying the CFL condition, the convergence order is calculated with DG-P 1 and DG-P 2 methods. The results are summarized in tables 5(a) and 5(b) respectively. One should notice that, for higher approximation orders, the convergence order should remain bounded to 2, given second-order accuracy of the LF2 time approximation. A visual representation of the time evolution L 2 error for different orders of approximation is displayed on figure 6. This case permits to validate the dispersive media DG code for the Drude model. Nevertheless, it does not represent any physical situation. The following sections will therefore present more realistic cases. 6.2 Near-field enhancement of a gold nanosphere

Presentation

A gold sphere of radius R and centered at (0, 0, 0) immersed in vacuum is considered, within a spherical domain. Silver-Muller absorbing boundary conditions (ABC) are used, and the sphere is illuminated with a sinusoidal plane wave propagating in the ẑ direction, which amplitude is modulated in time with a gaussian profile. Its central frequency is denoted as f c , and τ represents its initial phase delay :

E inc (t) = sin (2πf c (t -4τ )) e -( t-4τ τ ) 2 e x (53)
The permittivity of gold is considered to follow a Drude model, which parameters are summed up in table 6. The situation is represented on figure 7. In this case, analytical solutions exist in the frequency domain. Their expressions result from the Mie theory, which is extensively described in [vdH81]. The nearfield is here considered, and the sphere radius is taken equal to 20 nm.

Results

We focus on the amplitude of the total field in the vicinity of the gold nanosphere. In this case, the Mie theory predicts an enhanced field at the poles of the sphere. The reference Inria z y x solution used is given by a Matlab script7 which also exploits a Drude model. The reference solution being a frequency domain one, the DFT of the DG solution at the central frequency is calculated at each point of the domain.

As can be seen on figure 8, there is a very good agreement between the reference and numerical solutions, regarding both the maximum field value and the spatial extension of the resonance lobes. A more precise analysis of the results can be made by looking at the evolution of the modulus of the electric field along the contour line of equation {y = 0, z = 0}, which crosses the nanosphere at its center along the x+ direction. Such plots have been made for various mesh refinements and orders of interpolation, and are presented in figure 9. The L 1 error levels for these plots can be found on table 7, along with the meshes caracteristics. One could see that increasing refinement and interpolation order lead to a sharper field jump at the sphere interface. It also seems that neither of them manage to obtain a significantly better fit of the decreasing field outside of the sphere. This might be a consequence of the first-order geometrical approximation of the sphere, but that point remains to be verified. -As expected, increasing the interpolation order from 1 to 2 implies a roughly comparable relative increase in the computational cost (approximately 300%). On the M 1 mesh, this leads to a 25% drop in the L 1 error, whereas on the M 2 mesh, it only lowers of 12%. This is a logical consequence of the first-order approximation of the geometry of the sphere by straight-edged tetrahedrons, and the use of curvilinear elements will be considered in a future work. Another improvement could consist in a p-adaptative formulation, which would get the best from both high interpolation orders and mesh refinement, using large cells with high p values away from the interface, and small cells with low p values in the vicinity of the interface.

Remark : It has been noticed that a good numerical approximation of the field jump across the sphere interface also relies on a good homogeneity of the surface mesh of the sphere. The presented results were made using a particularly homogenous surface mesh for the sphere, as can be seen on figure 10(a). An example of a poor quality surface mesh is also shown.

Conclusion

On a physical point of view, it should be noted that the resonance phenomenon presented hereabove is a consequence of the dispersive behaviour of the metal only, and therefore does not appear when using a PEC sphere. These kinds of resonances can be combined in numerous ways (in [START_REF] Teixeira | Time-domain finite-difference and finite-element methods for Maxwell equations in complex media[END_REF], the reader can find the description of a L-junction device made of a linear assembly of metallic spheres; an extensive description of split-ring resonators is presented in [START_REF] Timothy | Analysis of metallic nanostructures by a discontinuous Galerkin time-domain Maxwell solver on graphics processing units[END_REF]; although, the latter references only represent a fraction of the litterature dealing with these kind of devices, and the reader can find many more without difficulties) with various geometries to achieve useful nanoscale devices. This textbook case, limited to a simple model and geometry, opens the way to more complicated geometries, models and methodology developments. While the following section is devoted to a generalized dispersive model, curvilinear elements and p-adaptative formulations will be at the center of future efforts.

Inria 7 Reflection coefficient of a silver slab described by various dispersion models

Presentation

This section presents a validation case for the DGTD formulation presented in section 3.3.2. As shown on figure 11, a silver slab is illuminated by a plane wave, described by equation (53). The central frequency f c and initial phase delay τ are respectively chosen equal to 900 THz and 4 10 -16 s in order to obtain an acceptable spectral density over the whole chosen frequency range, which is [300, 1500] THz. The slab is parallelepipedic, its side length being equal to 150 nm, and its thickness l to 10 nm. The incident field on the slab is noted E i , and the reflected field E r . The reflection coefficient is then defined as :

R = |E r | |E i | .
For a monochromatic wave of given frequency f m and amplitude 1, the reflection coefficient can be calculated analytically. Let ε s (f c )andµ s (f c ) be the permittivity and permeability of the slab at frequency f c , and ε 0 , µ 0 those of vacuum. Their respective impedances are therefore expressed as :

Z s (f c ) = µ s ε s (f c ), Z 0 = µ 0 ε 0 .
A few lines of hand calculation to write the matching conditions of the fields at the two interfaces between the slab and the vacuum leads to : R(f c ) = g 0s + g s0 e -2ikl 1 + g 0s g s0 e -2ikl , where :

g 0s = Z s -Z 0 Z s + Z 0 = -g s0 ,
and k is the exctinction coefficient of the material, e.g. k(f c ) = ε 0 ε 2 (f c ).

During the simulation, the time evolution of the diffracted field is recorded at a given probe point. Then, its normalized Fourier transform is calculated and compared to the analytical model. The dispersive behavior of the silver slab is modelized successively by a Drude model, a Drude-Lorentz model, a 2-SOGP model and a 4-SOGP model. The coefficients for these models and several others, calculated by the SA algorithm, can be found in the appendix A. For each of these models, the analytical reflection spectrum is calculated by the classical impedance formulas described above. Therefore, each numerical solution is compared to its exact solution, and its quality assessed through the calculation of the L 1 error over the reflection spectrum. Moreover, several of these results are also compared with the exact solution of the problem obtained using the raw Johnson & Christy data [START_REF] Johnson | Optical constants of the noble metals[END_REF].

Results

For all the following results, the used mesh consists of 58,826 vertices and 318,318 tetrahedrons, and a P 2 approximation has been used. We first assess the method by comparing of each model with its own exact solution. The calculated errors can be found in table 8, and plots for the Drude-Lorentz and the 2-SOGP models are presented on figure 12. One immediatly notices a very good agreement between the numerical and the exact solution over the whole frequency range, and a steady level of error regarding the number of second-order poles used. The impact of the number of poles on the computational time and on the amount of memory allocated was also assessed, and a visual representation of its influence is presented on figure 13. The time and memory values are normalized respectively by a reference time, t r , and a reference memory value, m r , which corresponds to the zero-SOGP case, e.g. all the computational domain is made of vacuum. The displayed computational times correspond to a 64 CPU cores parallel case, whereas the allocated memory values were recorded in a sequential case. For the above-described mesh and for a physical calculation time of 2 10 -1 s, one obtains t r = 3267s and m r = 958 Mo. Several things are to be noticed about this plot : first, the additional computational time when switching from non-dispersive to dispersive behaviour is quite high (approximately 30%).

On the contrary, the time cost of each additional pole is almost negligible. This might seem illogical at first sight, but can be easily explained : when a dispersive material is considered, the fields update requires a loop over all the tetrahedrons of the domain. For each of these, a test must be done to determine wether or not the current tetrahedron is in a dispersive domain, and if yes in which one. Considering the high cost of an IF statement done for each tetrahedron at each timestep, the additional cost observed on figure 13 is now understandable. The memory allocation rises regularly with the addition of poles, and the cost of a single SOGP can be evaluated to roughly 15 % of the reference memory occupation m r . An improved implementation is currently under study, that should particularly reduce the extra computational time required when dispersive materials are considered. A last comparison is made, showing the main interest of the generalized dispersive model : on figure 14, the plots of silver imaginary part permittivity and the predicted reflectance spectra are presented aside for J&C, 2-SOGP and 4-SOGP data. The positive impact of a good fitting of the material properties is particularly visible here at the resonance frequency, since the relative error on the amplitude of the resonance drops from 90% for 2-SOGP fitting to roughly 33% for 4-SOGP, for an extra computational time of 5% only.

Conclusion

The interest of using a generalized dispersive model has been shown in the particularly simple situation of a reflective slab. It has been shown that very good approximations of the real dispersive behaviour of metals such as silver can be obtained with a limited number of poles, and the cost of such improvements have been assessed. As an aside, improved implementations of the model might be worth considering in order to reduce the latter. The interest of mixed ZOGP/FOGP/SOGP fitting of permittivity functions remain to be evaluated.

Inria 8 Conclusion

This document presents the theoretical and numerical aspects of the Maxwell equations coupled to (i) a Drude model and (ii) a generalized dispersive model, discretized by a DGTD method. The physical dispersion basics are first presented, followed by a description of the DGTD method and the resulting formulation in cases (i) and (ii). Then, a priori stability and convergence results are proved for both models. Eventually, three numerical experiments are led, demonstrating the interest of correctly describing the dispersive behavior of metals in nanoscale devices.

Computing the interaction of electromagnetic waves in the THz range with nanoscale devices have proved to be a greedy kind of problems : considering a 50 nm-sized device and a 600 THz incident plane wave (typical values of these kinds of problems), one readily finds that the classical criterion of meshing with a precision of roughly h ≃ λ 10 gives h ≃ 50 nm, which is the size of the considered scatterer. This draws several conclusions : considering the latter statement, one would be tempted to mesh the computational domain with a high discrepancy in the cell sizes between the diffracting object and the vacuum (or other medium) that surrounds it : this is to be taken with care, for this could be the cause of spurious reflections and of a very small timestep. In many situations, a good way of overcoming such problems is to use a higher geometrical approximation of the considered object : this would allow to loosen the mesh and to put higher orders of approximation to good use. Furthermore, such problems leading to very high computational efforts, high performance computing (HPC) options are to be considered with even more interest. These strategies will be at the center of future efforts.

  Real and imaginary parts of the silver relative permittivity predicted by the Drude model compared to experimental data . . . . . . . . . . . . . . . . . . . . 2 Real and imaginary parts of the silver relative permittivity predicted by the Drude-Lorentz model compared to experimental data . . . . . . . . . . . . . . . 3 Real and imaginary parts of the silver relative permittivity predicted by the 2-SOGP model compared to experimental data . . . . . . . . . . . . . . . . . . . 4 Update scheme in the genrally dispersive case . . . . . . . . . . . . . . . . . . . 5 Exact and calculated E x fields with a P 1 approximation for the artificial case . . 6 L 2 error for different orders of approximation, for the artificial case . . . . . . . 7 Physical situation : gold nanosphere illuminated with a plane wave . . . . . . . 8 Comparison of DG and Mie solutions for the near-field enhancement of a gold nanosphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Mie and DG 1D plot of the electric field modulus across the dispersive gold nanosphere . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Comparison of two surfacic meshes for the gold nanosphere . . . . . . . . . . . 11 Physical situation : silver slab illuminated with a plane wave . . . . . . . . . . . 12 Reflection coefficient of a silver slab described by several dispersion models . . 13 Computational time and memory allocation for the generalized dispersive formulation for various models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Comparison of the fitting of imaginary part of the silver permittivity by 2-SOGP and 4-SOGP and its effect on the prediction of the reflectance spectrum . . . . 15 Real and imaginary parts of the silver and gold relative permittivity predicted by the 4-SOGP model compared to experimental data . . . . . . . . . . . . . . . List of Tables 1 Drude parameters for silver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Drude-Lorentz parameters for silver . . . . . . . . . . . . . . . . . . . . . . . . . 3 Parameters of the generally dispersive permittivity function for the Drude, Drude-Lorentz and Debye models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2-SOGP parameters for silver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Convergence orders for the artificial case with P 1 and P 2 approximations . . . . 6 Parameters set for the gold nanosphere case . . . . . . . . . . . . . . . . . . . . 7 L 1 errors for various meshes and orders of approximation in the case of the gold nanosphere nearfield enhancement . . . . . . . . . . . . . . . . . . . . . . . 8 L 1 errors, computational times and allocated memory for various dispersive model in the silver slab case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Coefficients of various dispersive models for silver . . . . . . . . . . . . . . . . . 10 Coefficients of various dispersive models for gold . . . . . . . . . . . . . . . . .

Figure 1 |

 1 Figure 1 | Real and imaginary parts of the silver relative permittivity predicted by the Drude model compared to experimental data from Johnson & Christy.

Figure 2 |

 2 Figure 2 | Real and imaginary parts of the silver relative permittivity predicted by the Drude-Lorentz model compared to experimental data from Johnson & Christy.

Figure 3 |

 3 Figure 3 | Real and imaginary parts of the silver relative permittivity predicted by the 2-SOGP model compared to experimental data from Johnson & Christy.

Figure 4 |

 4 Figure 4 | Update scheme in the generally dispersive case.

  |uπ h (u)| m,Ti ≤ Ch min(s,p)+1-m Ti u s+1,Ti , uπ h (u) 0,∂Ti ≤ Ch min(s,p)+ 1 2 Ti u s+1,Ti .

  y, z) = cos(x m x) sin(x m y) sin(x m z), g(x, y, z) = sin(x m x) sin(x m y) cos(x m z), h(x, y, z) = sin(x m x) cos(x m y) cos(x m z), i(x, y, z) = cos(x m x) sin(x m y) cos(x m z), j(x, y, z) = cos(x m x) cos(x m y) sin(x m z).

Figure 5 |

 5 Figure 5 | Exact and calculatedEx fields with a P 1 approximation for the artificial case.

Figure 7 |

 7 Figure 7 | Physical situation : gold nanosphere illuminated with a plane wave.

Figure 8 |Figure 9 |

 89 Figure 8 | Comparison of DG and Mie solutions for the near-field enhancement of a gold nanosphere.

Figure 10 |

 10 Figure 10 | Comparison of two surfacic meshes for the gold nanosphere.

Figure 11 |

 11 Figure 11 | Physical situation : silver slab illuminated with a plane wave.

Figure 12 |

 12 Figure 12 | Reflection coefficient of a silver slab described by 12(a) a Drude-Lorentz model and 12(b) a 2-SOGP model.

Figure 13 |Figure 14 |

 1314 Figure 13 | Computational time and memory allocation for the generalized dispersive formulation for various models.

Table 1 |

 1 Drude parameters for silver in the [300, 1500] THz frequency range.

	ε∞	ω d	γ d
	-	GHz	GHz
	3.7362	1.3871 10 7	4.5154 10 4

Table 2 | Drude-Lorentz parameters for silver in

 2 the [300, 1500] THz frequency range.

	ε∞	ω d	γ d	∆ε	ω l	γ l
	-	GHz	GHz	-	GHz	GHz
	2.7311	1.4084 10 7	6.6786 10 3	1.6336	8.1286 10 6	3.6448 10 6

Table 3 | Parameters of the generally dispersive permittivity function for

 3 the Drude, Drude-Lorentz and Debye models.

Table 4 | 2-SOGP parameters for silver in

 4 the [300, 1500] THz frequency range.

	ε∞	c 1	d 1	e 1	f 1	c 2	d 2	e 2	f 2
	-	GHz 2	GHz	GHz 2	GHz	GHz 2	GHz	GHz 2	GHz
	1.2944	1.8909 10 14	2.6584 10 6	0.0	0.0	5.6165 10 13	1.2005 10 7	4.3932 10 13	3.1709 10 6

  This leads to seek the values of E n+1

	i	and H	n+ 3 2 i

•

  The first one at time t n and time t n+1 with ϕ h = H The last one at time t n and time t n+1 with φ h = J

	n+ 1 2 i	;
	• The second one at time t n with ψ h = E [n+ 1 2 ]
	n+ 1 2 i	.
	Combining these different expressions yields :

i ;

•

Table 5 | Convergence orders with

 5

	(a) Convergence rate with P 1	(b) Convergence rate with P 2
	Refinement	Convergence rate	Refinement	Convergence rate
	1 50	-	1 25	-
	1 75	1.2575	1 50	2.2004
	1 100	1.1197	1 75	2.0826
	1 125	1.1000	1 100	2.0366
	1 150	1.0614	1 125	2.0432
	RR n°8298

P 1 and P 2 approximations. The total L 2 error is calculated with the discrete equivalent of expression (51).

error for different orders of approximation.

Table 6 | Parameters set for the gold nanosphere case.

 6

	ε∞	ω 0	γ d	fc	τ
	-	GHz	GHz	GHz	s
	1	1.19 10 7	1.41 10 5	4.5 10 5	2 10 -15

Table 7 | L 1 errors for various meshes and orders of approximation in the case of the gold nanosphere nearfield enhancement.

 7

	Mesh	ns	nt	r	L 1 P 1	L 1 P 2
	M1	26128	153517	17.66	1.1081 10 -8	8.3051 10 -9
	M2	146818	881154	19.77	7.6520 10 -9	6.6986 10 -9
	M3	389955	2338433	24.77	6.6025 10 -9

Table 8 |

 8 L 1 errors,

computational times and allocated memory for various dispersive model in the silver slab case.

	Dispersion model	L 1 error	t(s)	t tr	m mr
	Vacuum	-	3267	1	1
	Drude	0.0804	4300	1.316	1.151
	Drude-Lorentz	0.0765	4366	1.336	1.305
	2-SOGP	0.0820	4340	1.328	1.305
	4-SOGP	0.0941	4571	1.399	1.608

RR n°8298

τ f is the electron mean free path.

ne stands for the electronic density.

σ d represents the conductivity of the considered electrons.RR n°8298

Indeed, γ d represents the friction experienced by the electrons.Inria

It is of course possible to add more Lorentz poles, in order to describe more electronic transition resonances.RR n°8298

Code developed by Guangran Kevin Zhu, available at http://www.mathworks.com/matlabcentral/ fileexchange/31119-sphere-scattering RR n°8298

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

A Coefficients for the generalized dispersive formulation

This appendix provides dispersion coefficients for the generalized dispersive formulation for silver and gold over the [300,1500] THz frequency range. The experimental data, from [START_REF] Johnson | Optical constants of the noble metals[END_REF], is fitted to different models with a simulated annealing algorithm [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF]. For each of the two metals, Drude, Drude-Lorentz, 2-SOGP and 4-SOGP coefficients are given. To ease the reading, special units are used : ω's, γ's, d's and f 's are given in P Hz, whereas c's and e's are given in P Hz 2 . Plots of the 4-SOGP permittivities are displayed in figure 15.

RR n°8298