Marc Duruflé

Alexandre Gras

Philippe Lalanne

Non-uniqueness of the Quasinormal Mode Expansion of Electromagnetic Lorentz Dispersive Materials

Keywords: electromagnetic resonance, quasinormal mode, microcavity, nanoresonator, modal expansion

Any optical structure possesses resonance modes and its response to an excitation can be decomposed onto the quasinormal and numerical modes of discretized Maxwell's operator. In this report, we consider a dielectric permittivity that is a N-pole meromorphic function of the pulsation ω (Lorentz model). We propose a common formalism and obtain different formulas for the modal expansion. The non-uniqueness of the excitation coefficient is due to a choice of the linearization of Maxwell's equation with respect to ω and of the form of the source term. We make the link between the numerical discrete modal expansion and analytical formulas that can be found in the literature. We detail the formulation of dispersive Perfectly Matched Layers (PML) in order to keep a linear eigenvalue problem. We also give an algorithm to regain an orthogonal basis for degenerate modes. Numerical results validate the different formulas and compare their accuracy.

Non-unicité de la décomposition en modes quasi-normaux d'un champ électromagnétique pour des matériaux dispersifs régis par le modèle de Lorentz

Résumé : Un dispositif optique possède des modes de résonances et sa réponse à une excitation peut être décomposée sur les modes quasi-normaux et les modes numériques de l'opérateur de Maxwell discrétisé. Dans ce papier, on considère une permittivité diélectrique qui est une fonction méromorphe de ω avec N pôles (modèle de Lorentz). Nous proposons un formalisme géneral pour obtenir différentes formules pour la décomposition selon les modes. La non-unicité du coefficient d'excitation est dû au choix de la linéarisation des équations de Maxwell par rapport à ω et de la forme du terme source. Nous établissons le lien entre la décomposition modale discrète et des formules analytiques qui peuvent être trouveées dans la littérature. Nous détaillons la formulation avec des PML (Perfectly Matched Layers) dispersives afin de conserver un problème aux valeurs propres linéaire. Nous donnons aussi un algorithme pour retrouver une base orthogonale pour les modes dégénérés. Des résultats numériques valident les différentes formules et comparent leur précision.

Mots-clés : résonance électromagnétique, mode quasi-normale, microcavité, nano-résonateur, décomposition modale B Computation of biorthogonal vector for 2-D PML C Computation of biorthogonal vector for 3-D PML 1 Introduction

Quasinormal Modes

Optical micro and nanoresonators, be they plasmonic, photonic or hybrid, enhance and localize the electromagnetic energy at wavelength or subwavelength scales and are key components in many photonic applications. Their optical response is characterized by one of a few resonant features resulting from the excitation of one or a few dominant modes, the natural resonance modes of the resonators. These modes conveniently labelled by the integer m = 1, 2... are characterized by their electric and magnetic field vectors distributions, E m (r) and H m (r). These vectors are solutions of the following eigenvalue boundary problem [START_REF] Lalanne | Light interaction with photonic and plasmonic resonances[END_REF])

       -i ωm ε(ω m )E m -∇ × H m = 0, -i ωm µ(ω m )H m + ∇ × E m = 0, + Boundary conditions (1)
where ε(ω) and µ(ω) are respectively the dielectric permittivity and magnetic permeability and depend on the position r and pulsation ω. The exp(-iωt) convention for time harmonic fields is assumed. The fields E m (r) have continuous tangent traces across interfaces between subdomains and satisfy the outgoing-wave conditions at infinity. They are often called quasinormal modes (QNMs) to emphasize that their harmonic evolution is characterized by an exponential damping in time (they are the eigenstates of a non-Hermitian operator), so to say their pulsation ωm is complex with Im(ω m ) < 0. Micro and nanoresonators play a leading role in many areas in nanophotonics, from quantum information processing to ultrasensitive biosensing, nonlinear optics, and various optical metasurfaces. This pushes a strong pressure on the development of QNM theory and QNM numerical methods that explicitly consider QNMs in the analysis, providing important clues towards the interpretation of the resonator response.

Quasinormal Mode expansion of the scattered field

The scattered field [E S (r, ω), H S (r, ω)] is solution of time-harmonic Maxwell's equations As a result, the Maxwell's equations that will considered in the sequel are given as

       -i ω ε(ω)E S -∇ × H S = iω(ε(ω) -ε b ) E inc ,
       -i ω ε(ω)E S -∇ × H S = J, -i ω µ(ω)H S + ∇ × E S = 0, + Sommerfeld condition .
(2)

Efficiently computing this scattered field for a large number of pulsations consists in expanding the solution into the QNM basis :

E S (r, ω) = m α m (ω)E m (r, ω),
where the α m 's are the complex modal excitation coefficients, which measure how much the QNMs are excited by the driving field illuminating the resonator with a real frequency ω.

Inria

There are still some complicated mathematical issues in relation with the actual physical problem for which the open space is infinite and Maxwell's operator are continuous. For instance, the conditions under which the completeness of the QNM expansions of Eq. ( 1) is guaranteed are not still fully understood [START_REF] Colom | Modal expansion of the scattered field: Causality, nondivergence, and nonresonant contribution[END_REF][START_REF] Abdelrahman | Completeness and divergence-free behavior of the quasi-normal modes using causality principle[END_REF]. There are also several known and correct expressions for the α m 's [START_REF] Lalanne | Light interaction with photonic and plasmonic resonances[END_REF]), but we do not know which offer the best performance, e.g. the fastest convergence rate towards the actual solution as the number of QNMs retained in the expansion increases.

However, for practical geometries of interest in nanophotonics, the QNMs are computed numerically and it would be unrealistic to expect computing many QNMs over a broad spectral range, ideally in the entire lower half-plane of the complex plane (Im(ω m ) < 0). Rather we have to consider a discretized version of the initial Maxwell's equations and the physical domain is bounded by perfectly-matched layers (PMLs). The discretized operator is a matrix of finite dimension, and its spectrum is composed of a finite number of QNMs (often the relevant ones involved in the resonator dynamics in the spectral range of interest) completed by a large number of PML modes, which have much less physical significance but warrant completeness [START_REF] Vial | Quasimodal expansion of electromagnetic fields in open two-dimensional structures[END_REF][START_REF] Yan | Rigorous modal analysis of plasmonic nanoresonators[END_REF][START_REF] Lalanne | Light interaction with photonic and plasmonic resonances[END_REF].

Efficient QNMs solvers exist for computing and normalizing QNMs and PML modes for various geometries, such as plasmonic crystals, metal gratings and plasmonic nanoantennas [START_REF] Lalanne | Quasinormal mode solvers for resonators with dispersive materials[END_REF]; even freeware [START_REF] Bai | Efficient and intuitive method for the analysis of light scattering by a resonant nanostructure[END_REF] or improved commercial software packages [START_REF] Yan | Rigorous modal analysis of plasmonic nanoresonators[END_REF]) can be used. Thus the important remaining step is the reconstruction problem, i.e. the computation of the modal coefficients α m 's and the reconstruction of the scattered field.

While the concept of the QNM expansion as described by equation ( 1) is a common feature among papers on modal expansions, the expression of the excitation coefficient α m changes depending on the method used to project the scattered field onto the QNM basis. The contour integral method described in [START_REF] Binkowski | An auxiliary field approach for computing optical resonances in dispersive media[END_REF]Zschiedrich et al., 2018) necessitates an integration path in the complex frequency plane which encloses the eigenmodes of interests. The scattered field at the real frequency can thus be expressed as a sum of two contributions : a resonant contribution made up of contour integrals around the eigenfrequencies in the complex plane, and a non-resonant contribution, a setup that's echoed in [START_REF] Colom | Modal expansion of the scattered field: Causality, nondivergence, and nonresonant contribution[END_REF].

There also exists the resonant state expansion method described in [START_REF] Muljarov | Resonant-state expansion of dispersive open optical systems: Creating gold from sand[END_REF] where the modal expansion is obtained from the decomposition of the Green's function in QNMs. For dispersive materials whose permittivity can be described with a Drude model, the auxiliary field method and biorthogonal relations [START_REF] Yan | Rigorous modal analysis of plasmonic nanoresonators[END_REF] can be used to derive an analytical expression of α m . An improved formula for the field inside the resonator has been obtained in [START_REF] Wu | Intrinsic multipolar contents of nanoresonators for tailored scattering[END_REF]. For the same material permittivity, however, another formula for α m was found in [START_REF] Zolla | Photonics in highly dispersive media: the exact modal expansion[END_REF], without the use of auxiliary fields. Finally, we can cite also [START_REF] Zimmerling | A Lanczos model-order reduction technique to efficiently simulate electromagnetic wave propagation in dispersive media[END_REF], where the eigenmodes are solved simultaneously with the reconstruction, such that only significant eigenmodes are kept. This report will focus on the reconstruction problem of resonators composed of Drude-Lorentz materials. For these cases, we wish to make explicit possible sources of the non-uniqueness of the expressions of the excitation coefficients α m by using a common formalism. This common formalism is based on a discrete version of Maxwell's equations, which we will use to derive the different formulas. It was important for us to obtain these formulas from a discrete point of view in order to ensure that the numerical modal expansion converges indeed towards the expected solution. It will be shown that different linearizations of the nonlinear eigenvalue problem with the use of auxiliary fields yield different formulas of α m . With this approach, already known formulas are found. Furthermore, we show that an infinity of formulas can be derived easily by choosing appropriate source terms. Numerical results will show that the tested formulas converge towards the same reference solution. It will also be explained how degenerate eigenvalues (i.e. multiple eigenvalues) can be treated correctly with a simple Gram-Schmidt orthogonalization procedure.

In section 2, the Lorentz model used for the dielectric permittivity is detailed and the common formalism used to obtain different formulas for α m is introduced. In the sections 3 and 4.3, two different linearizations of Maxwell's equation are detailed providing different formulas of α m . The computational domain has to be truncated, e.g. with Perfectly Matched Layers. In order to keep a real eigenvalue problem, we have chosen to implement dispersive PMLs, that are detailed in subsection 4.6. Finally, numerical results are presented in section 5 in order to compare the accuracy of the different formulas that will have been rederived.

General setting

Drude-Lorentz model for the permittivity

The dielectric permittivity depends on the pulsation ω and can be approximated efficiently with Padé approximant, which usually necessitates a lower degree than a simple polynomial approximation. Assuming that all poles of the Padé approximation are simple (see (Sehmi et al., 2017b)), we can approximate the permittivity as

ε(ω) = L k=1 iσ k ω -Ω k
where Ω k are poles in the half complex plane (Im(Ω k ) < 0), σ k are complex coefficients and L is the number of poles. However, this first approximation is not physically relevant since it does not respect the causality principle that requires that:

ε(ω) = ε(-ω)
where ω denotes the conjugate of ω. Therefore, a physically relevant approximation is the following one (see for instance (Sehmi et al., 2017b))

ε(ω) = ε ∞ + D d=1 iη k ω + iγ k + L k=1 iσ k ω -Ω k + iσ k ω + Ωk (3)
where σ k and Ω k are complex coefficients and γk , η k are real coefficients. We assume that Im(ω k ) < 0 and γk > 0 such that all the poles of ε(ω) are located in the half complex plane Im(z) < 0. For the last sum of (3), we can reduce to the same denominator to obtain

iσ k ω -Ω k + iσ k ω + Ωk = - c k -iωσ k ω 2 -ω 2 0,k + iγ k ω
with the relations

ω 0,k = |Ω k |, γ k = -2Im(Ω k ), σk = 2Re(σ k ), c k = 2Im σ k Ωk .
The constants ω 0,k , γ k , σk and c k are real. For the first sum of (3), if we consider an even number D, we can regroup terms two by two:

iη k ω + iγ k + iη k+1 ω + iγ k+1 = - c k -iωσ k ω 2 -ω 2 0,k + iγ k ω Inria with the relations ω 0,k = γk γk+1 , γ k = γk + γk+1 , σk = η k + η k+1 , c k = η k γk+1 + η k+1 γk
Again, the constants ω 0,k , γ k , σk and c k are real. For metallic materials, good approximations are obtained by choosing D = 2 (see Sehmi et al. (2017b)). For other materials, we can take D = 0 (see (Sehmi et al., 2017a)). As a result, we will consider the following model of permittivity

ε(ω) = ε ∞ - L k=1 c k -iωσ k ω 2 -ω 2 0,k + iγ k ω , (4)
with real coefficients c k , σ k , ω 0,k , γ k . This model can also be obtained by considering the motion of electrons (see [START_REF] Wooten | Optical Properties of Solids[END_REF]) for which we have σ k = 0.

Discrete modal expansion

For a Drude-Lorentz model (as described in the previous sub-section) of ε(ω), auxiliary unknowns can be introduced in order to obtain a linear eigenvalue problem. After this linearization procedure and after discretization (e.g. with Finite Element Method), the time-harmonic Maxwell's Equations can be written

-iωM h U h + K h U h = F h , (5)
where M h is the mass matrix, K h is the stiffness matrix, and F h is the source term (h denotes the mesh size). U h is the main unknown that will contain components of E and other auxiliary unknowns introduced to obtain a linear eigenvalue problem. The matrices M h and K h are real and independent of ω, an example of matrices will be given in section 3. The right eigenvectors x m and left eigenvectors x ⊥ m solve the following linear eigenvalue problems which are transposed from each other

K h x m = iω m M h x m , K T h x ⊥ m = iω m M T h x ⊥ m .
From a discrete point of view, once the discrete linear system (5) is set, the biorthogonal projection of the unknown U h provides an unique formula for α m :

α m = 1 i(ω m -ω) F h , x ⊥ m , (6)
This biorthogonal projection is obtained by considering the relation (5) and taking the scalar product with the left eigenvector x ⊥ m . In this document, the convention x, y = x i y i is used (without conjugate). Because Maxwell's equations usually provide a complex symmetric linear system for the unknown E h (when other unknowns are substituted), the left eigenvector can be computed directly from the right eigenvector. The expressions will be given for the different linearizations introduced. The formula (6) holds if the left and right eigenvectors

x m , x ⊥ m are normalized such that M h x m , x ⊥ m = 1. (7)
As detailed in section 3, the discrete relation ( 6) can be rewritten with integrals in order to obtain analytic formulas for α m . With this procedure, we have obtained the following formulas (already known in the literature)

• The formula 5.11 in [START_REF] Lalanne | Light interaction with photonic and plasmonic resonances[END_REF]:

α m = 1 i(ω m -ω) Ωres J(r) • E m (r)dr (8)
which is obtained by a direct linearization of the system (2) (detailed in section 3).

• The formula proposed in [START_REF] Yan | Rigorous modal analysis of plasmonic nanoresonators[END_REF] (equivalent to formula 5.6 in [START_REF] Lalanne | Light interaction with photonic and plasmonic resonances[END_REF]):

α m = Ωres (ε b -ε ∞ )E inc • E m dΩ + ωm ωm -ω Ωres (ε(ω m ) -ε b )E inc • E m dΩ (9)
which is obtained by choosing a different source F h (detailed in section 4.1).

• The formula proposed in [START_REF] Zolla | Photonics in highly dispersive media: the exact modal expansion[END_REF] (equivalent to formula 5.10 in [START_REF] Lalanne | Light interaction with photonic and plasmonic resonances[END_REF]):

α m = ω i ωm (ω m -ω) Ωres J(r) • E m (r) dr (10)
which is obtained by starting from the second-order formulation of Maxwell's equations with curl-curl operator (detailed in section 4.3).

In these formulas, Ω res denotes the domain of the resonator for which ε(ω) is different from ε b (hence it is the support of the source term J). All these formulas hold if the modes E m are normalized as follows

Ω ∂(ω m ε(ω m )) ∂ ωm E m • E m - ∂ (ω m µ(ω m )) ∂ ωm H m • H m dΩ = 1. ( 11
)
where Ω is the computational domain. This is the usual normalization [START_REF] Muljarov | Resonant-state expansion of dispersive open optical systems: Creating gold from sand[END_REF]; [START_REF] Sauvan | Theory of the spontaneous optical emission of nanosize photonic and plasmon resonators[END_REF]; [START_REF] Bai | Efficient and intuitive method for the analysis of light scattering by a resonant nanostructure[END_REF]. An infinite set of formulas can be found by splitting the source on the different fields. For the linearization given in section 3, by writing the generalized source term as

F = [f 1 , f 2 , f 3,k , f 4,k ],
we can find the following generalization of the modal excitation coefficient:

α m = 1 i(ω m -ω) Ωres f 1 • E m -f 2 • H m +   L k=1 ω 2 0,k f 3,k + iω m f 4,k -ω 2 m + ω 2 0,k -iγ k ωm • E m   dΩ res (12)
provided that

f 1 + L k=1 iωf 4,k + ω 2 0,k f 3,k -ω 2 -iωγ k + ω 2 0,k + i ω ∇ × 1 µ b f 2 = J. ( 13
)
The derivation is detailed in section 4.2. With a particular choice of f 3,k , f 4,k (given in subsection 4.4), we obtain the following formula:

α m = 1 i(ω m -ω) Ωres (ε(r, ωm ) -ε ∞ (r)) (ε(r, ω) -ε ∞ (r)) J(r) • E m (r)dr (14)
This formula can also be obtained by expanding an auxiliary field into its modal components and computing E S from this auxiliary field. This process is detailed in subsection 4.4 and linked with alternative expressions of Green's kernel given in [START_REF] Muljarov | Resonant-state expansion of dispersive open optical systems: Creating gold from sand[END_REF] and another formula proposed by [START_REF] Wu | Intrinsic multipolar contents of nanoresonators for tailored scattering[END_REF].

The modal solution is given as

E modal S = N m=1 α m E m ( 15
)
where N is the number of modes conserved. The formulas (8), (10), ( 9), ( 12) and ( 14) for coefficients α m will provide a field E modal S that will converge to the scattered field E S when N tends to the size of matrix M h . Their convergence rate, however, may differ.

Inria 3 Eigenmode expansion for first-order formulation of Maxwell's equations

In this section, we note E, H the solutions of Maxwell's system (2). The permittivity ε(ω) is modeled by the relation (4) and the permeability µ = µ b is non-dispersive.

Discrete expansion

In order to linearize Maxwell's equations with respect to ω, auxiliary fields P k (known as the polarization) and Q k (known as the current) are introduced :

P k = c k -iωσ k ω 2 -ω 2 0,k + iγ k ω E, Q k = -iωP k
These auxiliary fields exist only inside the resonator. With elementary algebraic manipulations, we can reformulate Maxwell's system (2) as the following source problem

                     -i ω ε ∞ E + L k=1 Q k -∇ × H = J -i ω µ b H + ∇ × E = 0 -iωP k -Q k = 0 -iωQ k + γ k Q k + ω 2 0,k P k -c k E + iωσ k E = 0 (16)
inside the resonator. Outside the resonator, we consider homogeneous Maxwell's equations

-i ω ε b E -∇ × H = 0 -i ωµ b H + ∇ × E = 0 (17)
We can write this system using the linear operators K and

M KU -iωMU = F with K =     0 -∇× 0 1 ∇× 0 0 0 0 0 0 -1 -c k 0 ω 2 0,k γ k     , M =     ε e 0 0 0 0 µ b 0 0 0 0 1 0 -σ k 0 0 1     , U =     E H P k Q k     , F =     J 0 0 0     ,
where

ε e = ε ∞ in Ω res ε b , elsewhere.
The unknowns P k and Q k have degrees of freedom only inside the resonator. After discretization, the Maxwell's system is then given as

-iωM h U h + K h U h = F h ( 18
)
where

U h = (E h , H h , P h , Q h )
, and E h , H h , P h , Q h contain the components of E, H, P k , Q k on basis functions. The source term F h is given as

(F h ) i = Ωres J(r) • ϕ i (r) dr RR n°9348
where ϕ i are basis functions for unknown E h . Matrices M h and K h and basis functions ϕ i are given in appendix A.1. The right eigenvectors x m solve the eigenproblem

K h x m = λ m M h x m ( 19
)
where the eigenvalue λ m is linked with ωm by

λ m = iω m .
x ⊥ m is the left eigenvector of generalized eigenproblem (19) and solve the following eigenproblem

K T h x ⊥ m = λ m M T h x ⊥ m ( 20
)
Remark 1. The eigenvectors x ⊥ m and x m can be chosen such that they satisfy the following biorthogonality relation

M h x m , x ⊥ n = δ m,n where δ m,n is the Kronecker symbol. Proof. Assuming that M -1 h K h is diagonalizable, we have M -1 h K h = VDV -1
where D is a diagonal matrix with eigenvalues λ m on the diagonal and V the matrix whose columns are formed with right eigenvectors x m . The left eigenvectors of M -1 h K h denoted w m are the rows of matrix V -1 . Since VV -1 = I, vectors x m and w m are biorthogonal

x m , w n = δ m,n
The left eigenvectors w m can also be found by searching right eigenvectors of the transpose of

M -1 h K h : K T h M -T h w m = λ m w m By introducing x ⊥ m = M -T h w m , we obtain K T h x ⊥ m = λ m M h x ⊥
m , which corresponds to the eigenvalue problem (20). The biorthogonality relation is obtained directly by stating:

x m , w n = x m , M T h x ⊥ m = M h x m , x ⊥ m = δ m,n
The solution U h is expanded with right eigenvectors x m (they form a basis if the matrix

M -1 h K h is diagonalizable): U h = m α m x m
By injecting this expansion in (18) and using ( 19)), we obtain

m α m (-i ω + i ωm )M h x m = F h
The modal coefficient α m is directly obtained by taking the scalar product , with the left

eigenvector x ⊥ m α m (-i ω + i ωm ) = F h , x ⊥ m
As a result we obtain the relation

α m = 1 i(ω m -ω) F h , x ⊥ m (21)
which is the announced result in the sub-section 2.2.

Inria

Link with continuous expansion

The left eigenvector x ⊥ m of the system (18) can be obtained directly from the right eigenvector x m with the following theorem:

Theorem 1. The left eigenvector x ⊥ m = E ⊥ m , H ⊥ m , P ⊥ m , Q ⊥ m of the eigenproblem (19) is given as                      E ⊥ m = E m H ⊥ m = -H m P ⊥ m,k = ω 2 0,k c k -λ m σ k P m,k Q ⊥ m,k = - 1 c k -λ m σ k Q m,k
where

x m = (E m , H m , P m , Q m ) is the right eigenvector of (19).
Proof. The proof is written at a continuous level without replacing the curl operator by the associated finite element matrix. The transposition of this proof at a discrete level is actually immediate since the last three equations ( 18) are expressed directly on quadrature points as detailed in appendix A. The left eigenvector

x ⊥ m = E ⊥ m , H ⊥ m , P ⊥ m , Q ⊥ m solves the eigenproblem (20):                    λ m ε e E ⊥ m -λ m L k=1 σ k Q ⊥ m,k = ∇ × H ⊥ m - L k=1 c k Q ⊥ m,k λ m µ b H ⊥ m = -∇ × E ⊥ m λ m P ⊥ m,k = ω 2 0,k Q ⊥ m,k λ m Q ⊥ m,k = E ⊥ m -P ⊥ m,k + γ k Q ⊥ m,k (22)
We express P ⊥ m,k and P ⊥ m,k with only E ⊥ m :

P ⊥ m,k = ω 2 0,k λ 2 m + ω 2 0,k -γ k λ m E ⊥ m , Q ⊥ m,k = λ m λ 2 m + ω 2 0,k -γ k λ m E ⊥ m
The unknown H ⊥ m is also substituted to obtain a system in E m only:

λ 2 m ε e + c k -λ m σ k λ 2 m + ω 2 0,k -γ k λ m 1 R E ⊥ m + ∇ × 1 µ b ∇ × E ⊥ m = 0
where 1 R is the indicator function of Ω res . This is the same eigenvalue problem solved by E m . Therefore, with an appropriate normalization, we can set E ⊥ m = E m . By using the second equation of ( 22), we obtain H ⊥ m = -H m . Since, we have

P m,k = c k -λ m σ k λ 2 m + ω 2 0,k -γ k λ m E m , Q m,k = -λ m P m,k
we obtain

P ⊥ m,k = ω 2 0,k c k -λ m σ k P m,k , Q ⊥ m,k = - 1 c k -λ m σ k Q m,k .
The following theorem makes the link between the discrete formula ( 21) and the continuous formula (8):

Theorem 2. The formula (21) is the discrete equivalent of (8).

Proof. We use the fact that the source is non-null only for the equation in E:

F h , x ⊥ m = i E ⊥ m,i Ωres J(r) • ϕ i (r)dr where E ⊥ m,i is the i-th component of E ⊥ m .
By swapping the sum and the integral, we obtain

F h , x ⊥ m = Ωres J(r) • E ⊥ m (r)dr
Finally, we use that

E ⊥ m = E m (theorem 1) to obtain F h , x ⊥ m = Ωres J(r) • E m (r)dr.
We have a similar theorem for the normalizations:

Theorem 3. The normalization (7) with the left eigenvector given in theorem 1 is the discrete equivalent of (11).

Proof. Similarly to the previous proof, we have the following equality

M h x m , x ⊥ m = Ω ε e E m • E ⊥ m + µ b H m • H ⊥ m dΩ + Ωres L k=1 P m,k • P ⊥ m,k + (Q m,k -σ k E m,k ) • Q ⊥ m dΩ res
The left eigenvector is substituted by using theorem 1 and the relation λ m = iω m to obtain

M h x m , x m = Ω ε e E m • E m -µ b H m • H m dΩ + Ωres L k=1 ω 2 0,k c k -iω m σ k P m,k • P m,k - 1 c k -iω m σ k (Q m,k -σ k E m,k ) • Q m dΩ res P m and Q m are replaced to obtain M h x m , x m = Ω ε e E m • E m -µ b H m • H m dΩ + Ωres L k=1 (ω 2 0,k + ω2 m )c k + iσ k ωm (-2ω 2 0,k + iγ k ωm ) (-ω 2 m + ω 2 0,k -iγ k ωm ) 2 E m • E m dΩ res
A simple computation of the derivative of ωε(ω) with respect to ω provide the same expressions.

As a result, we have proven that

M h x m , x m = Ω ∂ (ω m ε(ω m )) ∂ ωm E m • E m -µ b H m cdotH m dΩ.
This relation proves that the normalization ( 7) is the discrete equivalent of (11).

Inria

Remark 2. The normalization can be written with only unknown E m :

M h x m , x m = Ω ∂ (ω m ε(ω m )) ∂ ωm E m • E m + ε(ω m )E m • E m dΩ.
The result is obtained directly by using the relation -

H m = 1 iωµ b ∇ × E m
ω 2 m Ω ε(ω m )E m • E m dΩ + Ω 1 µ b ∇ × E m • ∇ × E m dΩ = 0.
4 Derivation of other formulas and issues

Derivation with an alternative source

To obtain the formula (8), first, we have written Maxwell's equations directly for the scattered field E S (r, ω), H S (r, ω) and then introduced the auxiliary fields P and Q. In the aforementioned paper Yan et al. ( 2018), Maxwell's equations are first written for the total field, and the auxiliary unknowns P and Q are introduced at this step. Hence the unknowns E, H, P, Q solve the system ( 16) with J = 0. As a second step, we subtract the equations solved by the incident field (homogeneous Maxwell's equation with indices ε b and µ b ), and use the relations

[E(r, ω), H(r, ω)] = E S (r, ω) + E inc (r, ω), H S (r, ω) + H inc (r, ω)
to obtain the system solved by (E S , H S , P T,k , Q T,k ) inside the resonator:

                     -i ω ε ∞ E S + L k=1 Q T,k -∇ × H S = iω(ε ∞ -ε b )E inc -i ω µ b H S + ∇ × E S = 0 -iωP T,k -Q T,k = 0 -iωQ T,k + γ k Q T,k + ω 2 0 P T,k -(c k -iωσ k )E S = (c k -iωσ k )E inc (23)
The auxiliary fields P T,k , Q T,k here represent the total polarization and current whereas in the previous section we had introduced the scattered polarization P k and current Q k . Unlike the equations considered in section 3, we can see that the source term on the right hand side of the equations is no longer confined to the first equation. The coefficient α m becomes (by using formula (6) and biorthogonal vectors given in theorem (1)):

α m = 1 i(ω m -ω) Ωres iω(ε ∞ -ε b )E m • E inc + (c k -iωσ k )iω m (-ω 2 m + ω 2 0,k -iγ ωm ) E m • E inc dΩ res
that can be written as

α m = Ωres (ε b -ε ∞ )E inc • E m dΩ + ωm ωm -ω Ωres (ε(ω m ) -ε b )E inc • E m dΩ,
which is the announced formula (9). It is important to notice that the systems ( 23) and ( 16) provide exactly the same numerical solution E S (and also the same H S ). Only the auxiliary fields P and Q differ, that's why the source F h is different between the two approaches and two different formulas are obtained for α m . Other formulas for α m can be found by choosing a different distribution of the source over the four equations. This is the object of the next sub-section.

Remark 3. The formula (9) has been obtained through discretized Maxwell's equations and with biorthogonal projection. It has also been derived in a different way by using the divergence theorem and the continuous operator, not the discretized one, in [START_REF] Yan | Rigorous modal analysis of plasmonic nanoresonators[END_REF].

Generalized Sources

Let us split the source term J into a set of artificial sources denoted

f 1 , f 2 , f 3,k , f 4,k inside the resonator.                      -i ω ε ∞ E + L k=1 Q k -∇ × H = f 1 -i ω µ b H + ∇ × E = f 2 -iωP k -Q k = f 3,k -iωQ k + γ k Q k + ω 2 0,k P k -(c k -iωσ k )E = f 4,k (24)
By abuse of notation, we kept the same names for auxiliary fields H, P k , Q k , but these auxiliary fields are different from the fields introduced in section 3, only the field E will be the same. We can now obtain an infinity of formulas for α m :

Theorem 4. There exists an infinity of formulas for α m , e.g.

α m = 1 i(ω -ωm ) Ωres f 1 • E m -f 2 • H m +   L k=1 ω 2 0,k f 3,k + iω m f 4,k -ω 2 m + ω 2 0,k -iγ k ωm • E m   dΩ res

as soon as

Ωres

f 1 + L k=1 iωf 4,k + ω 2 0,k f 3,k -ω 2 -iωγ k + ω 2 0,k • ϕ i + i ωµ b f 2 • ∇ × ϕ i dΩ res = Ωres J • ϕ i dΩ res (25)
for all test-functions ϕ i . If the function f 2 is chosen continuous in Ω, we obtain the condition (13).

Proof. The fourth equation of ( 24) is multiplied by -iω and P is substituted to obtain:

Q k = -iωf 4,k -ω 2 0,k f 3,k -ω 2 -iωγ k + ω 2 0,k -iω (c k -iωσ k ) -ω 2 -iωγ k + ω 2 0,k

E

The second equation of ( 24) provides

-iωH = f 2 µ b - 1 µ b ∇ × E
We inject these expressions of H and Q k in the variational formulation of the first equation of ( 24) to obtain the variational formulation satisfied by E:

-ω 2 Ω ε(ω)E•ϕ i + 1 µ b (∇×E)•(∇×ϕ i ) dΩ = -iω Ω   f 1 + iωf 4,k + ω 2 0,k f 3,k -ω 2 -iωγ k + ω 2 0,k + i ωµ b f 2 • ∇ × ϕ i   dΩ

Inria

By equalling the right-hand side with -iω Ω J • ϕ i dΩ, we obtain the condition (25). This condition implies the continuous condition (13) if f 2 is continuous in Ω. The formula for α m is given by

α m = 1 i(ω -ωm ) Ω f 1 • E ⊥ m + f 2 • H ⊥ m + L k=1 f 3,k • P ⊥ m,k + f 4,k • Q ⊥ m,k dΩ.
By using theorem (1), we obtain the claimed expression for α m .

By choosing different splittings of the source (i.e. different functions f 1 , f 2 , f 3 , f 4 that satisfy the relationship above), we will obtain different formulas for α m . The modal solution obtained with these different formulas (see equation ( 15)) will converge towards the same electric field E S when the number of modes is increased.

Derivation with second-order formulation of Maxwell's equations

In this section we propose a different linearization of the problem by starting from the second order formulation. With this alternative linearization, we obtain the formula (10) for the coefficients α m . Let us start from the second-order formulation of Maxwell's equations

-ω 2 ε(ω)E + ∇ × 1 µ b ∇ × E = -iωJ.
In order to linearize this equation, let us introduce the field E = -iωE and the auxiliary field

P = (ε(ω) -ε ∞ ) E and Q = -iωP .
We obtain the following system of linear equations:

             -i ωE -E = 0 -i ωε ∞ E + L k=1 Q k + ∇ × 1 µ b ∇ × E = -i ωJ -i ωP k -Q k = 0 -iωQ k + γ k Q k + ω 2 0,k P k -(c k -iωσ k )E = 0 , ( 26
)
inside the resonator which gives the following stiffness and mass operators K and M for the vector U = [E, E , P , Q ]:

K =      0 -1 0 0 1 µ b ∇ × ∇× 0 0 1 0 0 0 -1 0 -c k ω 2 0,k γ k      , M =     1 0 0 0 0 ε e 0 0 0 0 1 0 0 -σ k 0 1     .
As a result, Maxwell's equations are rewritten as :

(-i ωM + K)U = F, where F = [0, -iωJ, 0, 0]
is the source term. After discretization, we have the following discrete system

(-iωM h + K h )U h = F h .
The matrices M h , K h are given in appendix A.2. For sake of simplicity, the following developments will be made with the continuous operators K and M, but the extension to the discrete operators M h and K h is straightforward. We obtain the following system of equations for the biorthogonal eigenvectors

(x ⊥ m = [E ⊥ , E ⊥ , P ⊥ , Q ⊥ ]):                    -i ωm E ⊥ + ∇ × 1 µ b ∇ × E ⊥ = 0 -i ωm ε ∞ E ⊥ -E ⊥ - L k=1 (c k -iω m σ k )Q ⊥,k = 0 -i ωm P ⊥,k + ω 2 0,k Q ⊥,k = 0 -i ωm Q ⊥,k + γ k Q ⊥,k + E ⊥,k -P ⊥,k = 0,
By eliminating the other variables, we can show that E ⊥ verifies

-ω 2 m ε(ω m )E ⊥ + ∇ × 1 µ 0 ∇ × E ⊥ = 0,
and subsequentially :

                     E ⊥ = E m E ⊥ = -iω m ε(ω m )E m P ⊥,k = ω 2 0,k ω 2 0,k -iγ k ωm -ω2 m E m Q ⊥,k = i ωm ω 2 0,k -iγ k ωm -ω2 m E m ,
where E m is the E-component of the the left eigenvector x m . We can now obtain the excitation coefficient :

α m = 1 i (ω m -ω) F h , x ⊥ m M h x m , x ⊥ m = -iω Ωres J(r) • E m (r)dr i(ω m -ω) N m ,
where the coefficient N m appears since we choose the normalization (11) of the first order formulation. N m is given as

N m = M h x m , x ⊥ m = Ω E m •E ⊥ +ε ∞ E m •E ⊥ dΩ+ Ωres L k=1 P m,k • P ⊥,k + Q m,k -σ k E m • Q ⊥,k dΩ res .
By substituting E m , P k , Q k , E ⊥ , E ⊥ , P ⊥ , Q ⊥ by their expression, we obtain

N m = -iω m Ω ε(ω m )E m • E m + ε e E m • E m dΩ + Ωres c k (ω 2 0,k + ω2 m ) + iσ k ωm (iγ k ωm -2ω 2 0,k ) (ω 2 0 -iγ k ωm -ω2 m ) 2 E m • E m dΩ res .

Inria

We recognize the normalization used by the first order formulation multiplied by -iω m . As a result, if E m is normalized by ( 11), we obtain that

N m = -iω m ,
which gives us this expression for the excitation coefficient:

α m = ω iω m (ω m -ω) Ωres J • E m dΩ.
We recognize the formula (10).

Obtention of other formulas by using auxiliary fields

With the introduced linearization of section 3, the polarization field P = k P k can be reconstructed inside the resonator with the modal expansion:

P modal = N m=1 α m P m = N m=1 α m (ε(ω m ) -ε ∞ )E m
and the scattered field E S computed from P modal :

E S ≈ P modal (ε(ω) -ε ∞ ) = N m=1 α m (ε(ω m ) -ε ∞ ) (ε(ω) -ε ∞ ) E m (27)
At least this modal expansion will converge inside the resonator towards the scattered field E S when the number of modes N increases. We can reinterpret the equation ( 27) as a new formula for α m :

α m = 1 i(ω m -ω) Ωres (ε(r, ωm ) -ε ∞ (r)) (ε(r, ω) -ε ∞ (r)) J(r) • E m (r)dr
which is the formula (14) announced in section 2.2. We can also choose the auxiliary field Q = k Q k and obtain the following formula

α m = 1 i(ω m -ω) Ωres ωm (ε(r, ωm ) -ε ∞ (r)) ω(ε(r, ω) -ε ∞ (r)) J(r) • E m (r)dr (28)
It turns out that the two formulas ( 14) and ( 28) are also valid outside the resonator. For the formula ( 14), it can be proven by choosing appropriately the sources f 3,k , f 4,k :

Theorem 5. The formula (14) provides a modal expansion valid outside the resonator.

Proof. We choose

f 1 = f 2 = 0, f 3,k = c k ω 2 0,k J ε(ω) -ε ∞ , f 4,k = -σ k J ε(ω) -ε ∞
in the system (24). We have

L k=1 iωf 4,k + ω 2 0,k f 3,k -ω 2 -iωγ k + ω 2 0,k = L k=1 c k -iωσ k -ω 2 -iωγ k + ω 2 0,k J ε(ω) -ε ∞ = J
As a result the condition ( 25) is satisfied and we can apply the theorem 4 to obtain a valid formula for α m in all the computational domain:

α m = 1 i(ω -ωm ) Ωres   L k=1 ω 2 0,k f 3,k + iω m f 4,k -ω 2 m + ω 2 0,k -iγ k ωm • E m   dΩ res which gives α m = 1 i(ω -ωm ) Ωres L k=1 c k -iω m σ k -ω 2 m + ω 2 0,k -iγ k ωm J • E m ε(ω) -ε ∞ dΩ res
The formula ( 14) is obtained.

Remark 4. The ratio

(ε(r, ωm ) -ε ∞ (r)) (ε(r, ω) -ε ∞ (r))
in formula ( 14) is located inside the integral in order to handle correctly the case of several resonators with different coefficients ε ∞ , c k , γ k , σ k , ω 0,k .

Case of an alternative source

In [START_REF] Wu | Intrinsic multipolar contents of nanoresonators for tailored scattering[END_REF], this process is applied with the alternative source of subsection 4.1. As a result, the following formula is obtained for the scattered field E S :

E S ≈ -E inc + N m=1 α W m (ε res (ω m ) -ε ∞ ) (ε res (ω) -ε ∞ ) (29)
where α W m is given by the equation ( 9):

α W m = Ωres (ε b -ε ∞ )E inc • E m dΩ + ωm ωm -ω Ωres (ε(ω m ) -ε b )E inc • E m dΩ
and ε res (ω) is the permittivity inside the resonator. The formula ( 29) is valid only inside the resonator and cannot be reinterpreted as a new formula for α m because of the presence of -E inc .

For a metallic material (e.g. D > 0, γ 0 = 0), the vector (0, 0, P, 0) is a static eigenvector of (19) and is excited with the right hand side chosen in (23). As a result, the formula (29) will not be used for metallic materials. We have observed numerically that it did not converge to the correct solution in this case.

Link with formulas with poles of permittivity

In [START_REF] Muljarov | Resonant-state expansion of dispersive open optical systems: Creating gold from sand[END_REF] (formula A9 of appendix A), the following alternative formula is proposed

α m = 1 i(ω m -ω) Ωres (ω -Ω k ) (ω m -Ω k ) J(r) • E m (r)dr (30)
where Ω k is a pole of the refractive permittivity ε(ω). This formula is proved at a continuous level, but we have also observed that it remains valid at a discrete level for any pole of ε(ω):

Theorem 6. The formula (30) provides a modal expansion valid outside the resonator.

Inria

Proof. The proof is similar to the proof of theorem 5. It starts from the following linearization (we consider the model (3) with D = 0) of Maxwell's equations

                 -iωε ∞ E - k P k -∇ × H = 0 -iωµ b H + ∇ × E = 0 -iω P k + iΩ k P k + iωσ k E = f 3,k -iω Q k -i Ωk Q k + iωσ k E = 0 (31)
This system has the same solution E as the system ( 16) by choosing

f 3,k = (-iω + iΩ k )J
The left eigenvector associated with the system (31) is given as

P ⊥ m,k = E m -iω m + iΩ k
As a result, we obtain the following formula for α m

α m = 1 i(ω m -ω) Ωres f 3,k • P ⊥ m,k dΩ res = 1 i(ω m -ω) Ωres (-ω + Ω k ) (-ω m + Ω k ) J • E m dΩ res
which is the formula (30). As a result this formula provides a valid modal expansion in all the computational domain.

The formula ( 14) can be seen as a combination of formulas (30) since all the poles of ε(ω) will be involved.

Treatment of degenerate eigenvalues

A set of degenerate modes {x k } m1≤k≤m2 , are solutions of the eigenvalue problem at the same eigenfrequency ωm1 . Degenerate eigenvectors do not necessarily form an orthogonal sub-basis with respects to M h . However, using Gram-Schmidt orthogonalization process, an orthogonal sub-basis with respects to M h can be constructed from the set of degenerate modes by algorithm 1. By applying this procedure, the formula (6) holds for degenerate eigenvalues with normaliza-

Algorithm 1 Algorithm to apply Gram-Schmidt orthogonalization to vectors x m for m=m 1 to m 2 do Initialize y = x m for j = m 1 to m -1 do Compute α = M h x m , x ⊥ j
Substitute y by y -αx j end for Compute left eigenvector y ⊥ from right eigenvector y by using theorem 1 Substitute x m by y/ M h y, y ⊥ Store x ⊥ j = y ⊥ / M h y, y ⊥ end for tion (7). This process can also be done with continuous eigenmodes by replacing

M h x m , x ⊥ j by Ω ∂(ω m ε(ω m )) ∂ ωm E m • E j - ∂ (ω m µ(ω m )) ∂ ωm H m • H j dΩ.
Here µ depends on ω inside the PML layers, which are detailed in the next sub-section. For elements inside the PML, the right eigenvector will be computed with formula (33) instead of theorem 1. In practice, we have observed degenerate eigenvalues due to PML layers or accumulation points. Accumulation points occur for poles and zeros of the permittivity ε(ω). In the 2-D transverse electric case (i.e. E = ue z ), we did not observe accumulation points for zeros of ε(ω).

PML

In this section, we describe how dispersive PMLs are handled. The damping coefficients σ x , σ y and σ z inside a PML where x > x 0 , y > y 0 or z > z 0 are parabolic:

σ 1 = σ x = 3 log(1000) 2a 3 (x -x 0 ) 2 v max σ σ 2 = σ y = 3 log(1000) 2a 3 (y -y 0 ) 2 v max σ σ 3 = σ z = 3 log(1000) 2a 3 (z -z 0 ) 2 v max σ.
The coefficient σ serves to adjust the reflection coefficient of the PML. v max is the speed of the wave inside the PML and a is the thickness of PML layer. In this section, we describe the formulation used for dispersive PMLs. The matrices M h , K h are no longer symmetric. We provide relations between the left eigenvector x ⊥ m and right eigenvector x m . As a result we do not need to compute the eigenvectors of the adjoint problem, since we can compute x ⊥ m directly from the right eigenvector x m .

2-D case

In Transverse Electric case, we have

E = u e z , H = v x e x + v y e y .
We use a split formulation of the PMLs where u = u 1 + u 2 inside the PML. The unknowns u 1 , u 2 , and v = (v x , v y ) are solutions of:

                     -iω ε b u 1 + ε b σ x u 1 - ∂v x ∂x = 0 -iω ε b u 2 + ε b σ y u 2 - ∂v y ∂y = 0 -iωµ b v + µ b σ x 0 0 σ y v -∇(u 1 + u 2 ) = 0
u = 0 at the border of the PML.

We consider the unknowns:

u = u 1 + u 2 u * = u 1 -u 2 .
Inria u, u * , v, are solutions of the following system,

             -iω ε b u + ε b σ x + σ y 2 u + ε b σ x -σ y 2 u * -div v = 0 -iω ε b u * + ε b σ x + σ y 2 u * + ε b σ x -σ y 2 u - ∂v x ∂x - ∂v y ∂y = 0 -iωµ b v + µ b σ v -∇ u = 0. (32)
The unknown u * exists only in the PML domain. In the physical domain, only unknowns u and v are present, and we solve

-iω ε(ω) u -div v = -iωj -iωµ b v + µ b σ v -∇ u = 0,
where j is the source term. Of course, additional unknowns p and q are added in Ω res to linearize the system in ω. After discretization, we will obtain :

-iωM h U h + K h U h = F h .
The matrix M h is symmetric, while K h is not. The left eigenvector x ⊥ m and the right eigenvector x m are written as:

x ⊥ m =   u ⊥ m u * ,⊥ m v ⊥ m   , x m =   u m u * m v m   .
We have obtained the following relations (λ m = iω m is the eigenvalue associated with x m and x ⊥ m ):

u ⊥ m = 1 - σ x + σ y 2λ m u m u * ,⊥ m = σ x -σ y 2λ m u m and v ⊥ m =     1 µ b (-λ m + σ x ) ∂u ⊥ m ∂x + ∂u * ,⊥ m ∂x 1 µ b (-λ m + σ y ) ∂u ⊥ m ∂y - ∂u * ,⊥ m ∂y     .
The proof is given in appendix B.

3-D case

In the PMLs we have:

                           -iωε b E + εT 2,3,1 E -∇ × H * = 0 -iωµ b H + µT 2,3,1 H + ∇ × E * = 0 -iωE * + T 3,1,2 E * + iωE -T 1,2,3 E = 0 -iωH * + T 3,1,2 H * + iωH -T 1,2,3 H = 0 E × n = 0 at the border of the PML, RR n°9348 with T i,j,k =   σ i 0 0 0 σ j 0 0 0 σ k   .
The unknowns E * and H * exist only in the PML domain. In the physical domain, there are only unknowns E and H (supplemented by unknowns P and Q in Ω res ) that solve (16). After discretization we will obtain:

-iωM h U h + K h U h = F h .
The matrices M h and K h are not symmetric (see appendix (C)). If we note x m = (E m , H m , E * m , H * m ) the right eigenvector, the left eigenvector x ⊥ m is given as:

x ⊥ m =            E * m -H * m 1 + T 2,3,1 -T 3,1,2 -λ m + T 3,1,2 ε b E m -1 + T 2,3,1 -T 3,1,2 -λ m + T 3,1,2 µ b H m            . ( 33
)
The proof is given in appendix C. Straightforward computations give that

M h x m , x ⊥ m = Ω ∂(ω m ε(ω m )) ∂ ωm E m • E m - ∂ ( ωm µ(ω m )) ∂ ωm H m • H m dΩ with ε(ω) = ε b (-iω + T 2,3,1 ) (-iω + T 3,1,2 ) -iω (-iω + T 1,2,3 ) , µ(ω) = µ b (-iω + T 2,3,1 ) (-iω + T 3,1,2 ) -iω (-iω + T 1,2,3 ) , (34)
inside the PML. We find the announced normalization (11) in the introduction.

Static modes

In the 2-D transverse electric case (E = u e z ), the static modes satisfy ∇u m = 0.

Moreover we impose a null Dirichlet condition on PML boundaries. As a result, the static modes satisfy u m = 0 on all the domain. This is also true at a discrete level because u is discretized with continuous finite elements. As a result, static modes have non-null components mainly for unknown v. The formulas (8), ( 9), ( 10), ( 14) and ( 28) can be used to reconstruct the unknown u if static modes are dropped.

For the 3-D case, the static modes satisfy

∇ × E m = 0.
They are usually constituted of gradients of functions. The associated eigenspace is usually large. In the numerical experiments, we have observed that 50% of the eigenmodes are static eigenmodes. Clearly, in large 3-D cases, it is not possible to compute and store all these static modes, that's why we have made the choice to compare formulas without these modes. Besides, we have observed that the matrix M -1 h K h is often not diagonalizable because of these static modes. Because of this lack of diagonalization, some eigenvectors are not correctly computed or the Gram-Schmidt process fails.

Inria 5 Numerical results

The numerical results have been obtained with the software Montjoie Duruflé (2018) for the computation of finite element matrices M h and K h given in section 3. The fields will be displayed inside the physical domain.

In this section, all the eigenvalues of the matrix M -1 h K h are computed with Lapack. We represent adimensionalized pulsations ω m defined as

ω m = ωm ω adim where ω adim = c 0 L 0 , L 0 = 10 -7 .
c 0 is the speed of light and L 0 the characteristical length (here 100nm). All of the eigenvalues such that |ω m | < 10 -6 are dropped in order to remove static modes. Since the eigenvalues are complex conjugate, only eigenvalues (and associated eigenvectors) such that Re(ω m ) ≥ 0 are stored. The eigenvalues such that λ m = σ i (σ i is the damping function in PMLs) are also excluded, since the auxiliary fields H, E * , H * cannot be eliminated (division by zero) for these eigenvalues. In practice, we have observed that the associated eigenvectors have null components (at machine precision) for the unknown E m and do not contribute to the field E S . Finally, if two pulsations ω i , ω j are close enough (i.e. |ω i -ω j | < 10 -6 ) they are considered degenerate.

In this section, we will compare the following formulas

• Formula (8) denoted as Usual, since it can be obtained from a direct linearization of Maxwell's equations.

• Formula (9) denoted as Alternative Source, since it can be obtained by choosing a different right hand side of Maxwell's equations.

• Formula (10) denoted as Order2, since it can be obtained from second-order formulation of Maxwell's equations.

• Formula ( 14) denoted as P-unknown, since it can be obtained by reconstructing the field E from the modal expansion of the unknown P.

• Formula (28) denoted as Q-unknown, since it can be obtained with the unknown Q.

• Formula (29) denoted as Alternative + P, since it can be obtained with the unknown P by choosing the alternative source. We observed that this formula does not work with metallic materials (an explanation is given in sub-section 4.4.1). As a result, this formula will be used for non-metallic materials only.

These formulas are rewritten in the table 1.

Table 1: Formulas used to computed the scattered field E S With the coefficients α m and β m defined as

β m = ω (ω m -ω) Ωres (ε(ω) -ε b )E inc • E m dΩ α m = Ωres (ε b -ε ∞ )E inc • E m dΩ + ωm ωm -ω Ωres (ε(ω m ) -ε b )E inc • E m dΩ

Formula Expression

Formula (8) (Usual)

E S = β m E m Formula (9) (Alternative Source) E S = α m E m Formula (10) (Order2) E S = ω ωm β m E m
Formula (14) (P-unknown)

E S = (ε(ω m ) -ε ∞ ) (ε(ω) -ε ∞ ) β m E m Formula (28) (Q-unknown) E S = ωm ω (ε(ω m ) -ε ∞ ) (ε(ω) -ε ∞ ) β m E m
Formula (29) (Alternative + P)

E S = -E inc + (ε(ω m ) -ε ∞ ) (ε(ω) -ε ∞ ) α m E m Inria

2-D golden disk

In the 2-D case, we consider a scalar electric field (corresponding to the z-component of E S ) and a vectorial magnetic field (corresponding to x and y-components of H S ). We first look at the case of the field diffracted by a metallic disk with a radius of 100 nm, where the metal is gold, which can be can modeled accurately with Drude-Lorentz model ( 3) with D = 2, L = 2. We have chosen the values computed in (Sehmi et al., 2017b) (table II column with L = 2), they are repeated in table 2. These constants are given in electron-volts, they must be multiplicated

ε ∞ η 1 γ1
η 2 γ2 2.6585 1056.9 0 -1056.9 0.07247

σ 1 Ω 1 σ 2 Ω 2 0.57604 + 0.18443i 2.5509-0.27427i 4.1891+4.2426i 2.8685-1.2195i Table 2: Constants ε ∞ , η k , γk , σ k , Ω k for gold. by 2πe
h where e = 1.602176634 • 10 -19 is the elementary charge and h = 6.62607015 • 10 -34 is the Planck constant. The physical computation domain is 400 nm long and 200 nm wide (see figure 1). PML layers are added to the mesh of figure 1. The thickness of PML is equal to 100nm with two cells in direction of PMLs. The damping of PMLs σ is taken equal to 3. The disk is made of gold while the surrounding medium is the vacuum. We represent in figure 2 Figure 1: Mesh used for the scattering of a disk the numerical adimensionalized pulsations ω m obtained with the mesh displayed in figure 1 and fourth-order finite elements. We can compare these pulsations with analytical QNMs for the disk (computed with Bessel functions). The comparison is displayed in figure 3. We see that QNMs are correctly computed, and we also observe the presence of other modes that we call PML modes. The matrices M h and K h have 6440 rows. Among the 2364 eigenvectors stored, 488 are associated with a degenerate eigenvalue. Degenerate eigenvalues are represented with red squares in figure 2. These degenerate eigenvalues are due to accumulation points (here due to poles of the dielectric permittivity) or due to PML layers. The field driving the system is a plane wave propagating along the x-axis, at the real frequency ω :

E inc = e ikx e z
where k = ω c 0 is the wave number. The scalar field u = E z is discretized with continuous finite elements (here Q 4 with the mesh of figure 1). The solution is plotted for two wavelengthes

Hankel functions).

We compute the field diffracted by the disk for 201 angular frequencies ω evenly spaced in the interval [ω 1 , ω 2 ] where

ω 1 = 2πc 0 λ 1 , ω 2 = 2πc 0 λ 2 .
The two wavelengthes λ 1 , λ 2 are equal to 1 µm and 250 nm. In figure 5, we display the relative error between the modal solution

E modal S = m α m Ẽm
and the direct FEM solution

E F EM S = (-iωM h + K h ) -1 F h
as a function of the spectrum width. For a given spectral width, the relative error is computed for the 201 angular frequencies and the maximum value of this error is retained and plotted. For a given spectral width L, only the modes whose eigenfrequencies ωm verify

Re(ω m ) ∈ [-L ω adim , L ω adim ] and Im(ω m ) ∈ [-ω adim L/2, 0]
are included in the expansion. The relative error is computed on the whole physical domain Ω p (PMLs are not included) by the formula

Relative Error = Ωp E modal S -E F EM S 2 dΩ p Ωp E F EM S 2 dΩ p
For this case, all the formulas converge nicely to the expected solution. The formula ( 14) is the most accurate while the formula (8) is the least accurate. Rela ive error Usual Al erna ive Source Order2 P unknown Q unknown

2-D silicon square

Figure 5: Relative error between the scattered field computed with the modal expansion and with a direct FEM solver as a function of the spectral width. Case of the golden disk.

ε ∞ σ 1 Ω 1 σ 2 Ω 2
σ 3 0.81568 1.6934 + 2.084i 3.3736 -0.11402i 5.2573 + 8.0106i 3.6519 -0.52378i -1.7164 + 5.9939i

Ω 3 σ 4 Ω 4 σ 5
Ω 5 4.287 -0.21116i -0.00528 + 0.32911i 5.3188 -0.18434i -3.8438 + 6.9298i 5.5064 -1.7892i Inria thickness of PML is equal to 100nm with one cell in direction of PMLs. The damping of PMLs σ is taken equal to 3. Eighth-order finite elements are used for this simulation. In figure 6, the eigenpulsations ω m are plotted in the complex plane. We see the accumulation points (with their multiple eigenvalues) that are located on poles of ε(ω). The matrices M h , K h contain 39640 rows. Among the 16933 eigenvectors stored, 4974 are associated with a degenerate eigenvalue.

The field driving the system is a plane wave propagating along the x-axis, at the real frequency ω :

E inc = e ikx e z
The total field is displayed in 7 for different wavelengthes. We see the change of behavior of silicon. This material has almost no damping for λ > 500nm and behaves almost like a metal for short wavelengthes.

In figure 8, we chose to represent the error inside the silicon resonator in order to compare the formula (29) with other formulas. The relative error is computed for 201 angular frequencies and the maximum value of this error is retained. Pulsations are evenly spaced in the interval [ω 1 , ω 2 ] where

ω 1 = 2πc 0 λ 1 , ω 2 = 2πc 0 λ 2 .
The two wavelengthes λ 1 , λ 2 are equal to 800 mm and 200 nm. We observed that the formula (29) does not work outside the resonator, but provides accurate results inside the resonator. For this case, the formula (8) is the least accurate. The large number of modes N is here due to the presence of accumulation points in the region of interest (five poles). Finally, in figure 9, the relative L 2 error inside the square is plotted versus the angular frequency ω when N = 8908 modes are used in the modal expansion. [500 nm, 1.5 µm] by using the function least_squares of Python package scipy.optimize. The obtained coefficients are given in table 4. The constants σ 1 and Ω 1 are given in electron-volts. In

ε ∞ σ 1 Ω 1 1.12648273 2.1759525332+20.775848892i 3.95095351-0.19089288i Table 4: Constants ε ∞ , σ k , Ω k for silicon.
figure 10, the permittivity has been represented for this reduced model (L = 1) and the original model given by table 3 (L = 5). We check that there is a good agreement until ω ≈ 2.5eV which corresponds to wavelengthes greater than 500 nm. The physical computation domain is the parallelepiped box [0, 150nm]×[0, 150nm]× [-150nm, 150nm] with a quarter of the dielectric ball (see figure 11). PML layers are added to the mesh of figure 11. The thickness of PML is equal to 100nm with only one cell in direction of PMLs. The damping of PMLs σ is taken equal to 1. The source is an incident plane wave oriented in z-direction and polarized in x-direction

E inc = e ikz e x
We impose a Perfectly conducting condition on plane x = 0 (i.e. E × n = 0) and a Neumann condition on plane y = 0 (i.e. H × n = 0) in order to have the same solution as for the whole sphere. Fourth order edge elements are used for the unknown E and the mesh of figure 11. We compute the field diffracted by the sphere for 201 angular frequencies ω evenly spaced in the interval [ω 1 , ω 2 ] where

ω 1 = 2πc 0 λ 1 , ω 2 = 2πc 0 λ 2 .
The two wavelengthes λ 1 , λ 2 are equal to 1.5 µm and 500 nm. The numerical solution obtained for 500nm is plotted in figure 12. The error between this solution and the analytical solution in 2-D, the formulas provide a modal solution that converges towards the direct FEM solution. For this case, the formula ( 14) is the most accurate while the formula ( 9) is the least accurate. Similarly to the 2-D case, only modes such that Re(ω m ) ∈ [-L ω adim , L ω adim ] and Im(ω m ) ∈ [-ω adim L/2, 0] are kept, where L is the spectral width. When a reduced spectrum is selected, the formula ( 14) is the most accurate. In figure 16, the relative error is computed on curl(E) but only inside the resonator. We can see that the formula (29) is quite accurate for small values of L, this formula does not converge because the incident field is projected on degrees of freedom by interpolation (instead of solving a linear system with the mass matrix). As a result, it stagnates at 1 % which is the level of accuracy of the interpolation for this case.

If the electric field is desired, a nice approach consists in discretizing H with edge elements (instead of E), reconstructing H with the modal expansion:

H modal = α m Hm .
and of computing E S by using Maxwell's equations

E S = 1 -iωε(ω) J + ∇ × H modal (35)
In figure 17, the relative error on the electric field has been computed by using this method. The formulas (8), ( 9), ( 14) or ( 28) can be used to obtain H modal with the coefficients α m (because they are obtained with f 2 = 0, see for instance theorem 5). The coefficients α m given by the formula (10) can be used only to reconstruct E modal (with equation ( 15)). The reason is that this formula has been established by introducing the unknowns E, E , P , Q (see section (4.3)). Thus, only these four unknowns can be reconstructed with this formula and not H. In figure 17, we observe that the reconstructed field E with this method converges correctly to the numerical eletrical field. For this case, the formula ( 14) is the most accurate.

Silver rod

We consider the case of a silver rod. The resonator is a cylinder of radius 240 nm and of height 400 nm. Only a quarter of the cylinder is meshed with Dirichlet/Neumann conditions to obtain the same solution as for the full cylinder. PML layers are added for x, y, z > 400nm and a Silver-Müller condition is set for z = -400nm. The silver material is modeled by a Drude-Lorentz model with D = 2, L = 2. In table 5, coefficients of this model are given, they have obtained with the function least_square as in the previous function. The measured values of ε have been taken from http://sspectra.com/sopra.html. In figure 18 the Silver-Müller condition by a PML for z < -400nm with the same parameters as for the PML z > 400nm (i.e. same damping σ and same thickness), degenerate eigenvalues are observed outside of poles and zeros of ε(ω).

The source is an incident plane wave oriented in z-direction and polarized in x-direction

E inc = e ikz e x
We compute the field diffracted by the sphere for 201 angular frequencies ω evenly spaced in the interval [ω 1 , ω 2 ] where

ω 1 = 2πc 0 λ 1 , ω 2 = 2πc 0 λ 2 .
The two wavelengthes λ 1 , λ 2 are equal to 900 nm and 400 nm. With this choice of wavelengthes, the real part of the permittivity does not belong to the critical interval -3, -1 3 . The mesh is very coarse, such that the numerical solution (see figure 20) is not accurate. The numerical error is estimated to be around 20 %. In figure 21, the different formulas are compared in the physical domain. We observe that the formula ( 14) is the most accurate for this case.

Conclusion

In this report, we have considered the solution of Maxwell's equations for dispersive materials modeled by a Drude-Lorentz model. A classical approach to compute the solution on a frequency range consists of computing the eigenmodes of the system (QNMs) and projecting the solution on these modes. We have shown how different linearizations of Maxwell's equations with respect to the pulsation ω can lead to several formulas for the projection coefficients α m . With this approach, we show how formulas classically found in the literature (e.g. ( 8), ( 9), ( 10)) can be proven at a discrete level. An infinity of formulas exist for the excitation coefficients α m , a general formula ( 12) is obtained by splitting the source into the different equations. The formulas ( 14), ( 30) and ( 28)) are obtained with a particular choice of splitting. We have established the link between the biorthogonal projection (that is usually the method of choice in numerics to project Inria the solution on eigenvectors) and analytic formulas that have been obtained in the continuous space. Because of the symmetry of Maxwell's equations, the left eigenvector can be computed from the right eigenvector without solving the transpose eigenvalue problem. The expressions of the left eigenvectors are provided by the theorem (1) for Drude-Lorentz material and by the equation ( 33) for dispersive PMLs. We provided in the algorithm 1 a simple way to handle correctly degenerate eigenvalues. These eigenvalues usually appear because of accumulation points (located at poles or zeros of the permittivity function) or because of a symmetry in the PML layout. Numerical experiments confirm that the formulas ( 8), ( 9), ( 10), ( 14), ( 30) and ( 28) can be used to reconstruct the electric field. If all the eigenmodes are included in the modal expansion, all these formulas provide the same solution, i.e. the finite element solution. In the tested cases, the formula ( 14) is the most accurate. In the case of non-metallic resonators (ω 0,k = 0), the formula ( 29) is of interest to obtain the field inside the resonator. In 3-D computations, since static modes are dropped, the electric field cannot be computed directly (only its curl is correctly computed). We have proposed to use formulas ( 8), ( 9), ( 14) or (28) to compute the magnetic field and recover the electric field by using Maxwell's equations (equation ( 35)). Zschiedrich, L., Binkowski, F., Nikolay, N., Benson, O., Kewes, G., and Burger, S. (2018). Rieszprojection-based theory of light-matter interaction in dispersive nanoresonators. Phys. Rev. A, 98:043806.

A Finite element matrices

A.1 First-order formulation

The computational domain Ω is made of hexahedra (in 3-D) :

Ω = K i
where K i is a single hexahedron. Edge elements based on Nédélec's first family are used for the unknowns E and discontinuous finite elements for unknowns H, P and Q :

E ∈ V h = u ∈ H(curl, Ω) such that ∀K i ∈ Ω, DF T i u • F i ∈ Q r-1,r,r × Q r,r-1,r × Q r,r,r-1 H, P, Q ∈ W h = u ∈ (L 2 (Ω)) 3 such that ∀K i ∈ Ω, u ∈ Q 3 r,r,r
where F i is the transformation from the unit cube [0, 1] 3 to the hexahedron K i , DF i its jacobian matrix and Q r1,r2,r3 is the following polynomial space

Q r1,r2,r3 = p(x) = r1 =0 r2 m=0 r3 m=0 a ,m,n x 1 x m 2 x n 3 a ,m,n ∈ C Inria
For the basis functions ϕ i generating the space V h , we use the basis functions detailed in [START_REF] Cohen | Non spurious spectral-like element methods for Maxwell's equations[END_REF]. For the basis functions ψ i generating the space W h , we use Lagrange interpolation polynomial

ψi1,i2,i3, (x) = ϕ GL i1 (x 1 )ϕ GL i2 (x 2 )ϕ GL i3 (x 3 ) where ϕ GL i (x) = j =i x -ξ j j =i ξ i -ξ j
with ξ j the Gauss-Lobatto points on the interval [0, 1]. The basis functions on the hexahedron K i are then given as

ψ i1,i2,i3 • F i = ψi1,i2,i3
The stiffness and mass matrices K h and M h associated with the system (16) are given by:

M h =     D E h 0 0 0 0 µ b D H h 0 0 0 0 D h 0 -σ k C T h 0 0 D h     K h =     S h -R h 0 C h R T h 0 0 0 0 0 0 -D h -c k C T h 0 ω 2 0,k D h γ k D h    
where h denotes the mesh size and

(D E h ) i,j = Ω ε e ϕ i (x) • ϕ j (x)dx (D h ) i,j = Ωres ϕ i (x) • ϕ j (x)dx (D H h ) i,j = Ω ψ i (x) • ψ j (x)dx (C h ) i,j = Ωres ϕ i (x) • ψ j (x)dx (R h ) i,j = Ω ψ j (x) • ∇ × ϕ i (x)dx (S h ) i,j = Γ S ε b µ b (ϕ i (x) × n) • (ϕ j (x) × n)dx
where n is the outgoing normale on the boundary Γ S where a Silver-Müller condition can be set.

Here, we consider ε e = ε ∞ in Ω res ε b , elsewhere. Degrees of freedom for P and Q are restricted to the domain Ω res . All the integrals are evaluated with Gauss-Lobatto points such that matrices D h and D H h are diagonal. For the 2-D case, the unknown E is scalar (we consider the Transverse Electric case) and discretized with nodal continuous elements, unknowns P and Q are also scalar and discretized with basis functions of E. The unknown H is vectorial and discretized with discontinuous elements. Basis functions are also based on Gauss-Lobatto to achieve mass lumping.

A.2 Second-order formulation

We use the same variational spaces and same basis functions as the previous section. The stiffness and mass matrices K h and M h associated with the system (26) are given by:

M h =     D I h 0 0 0 0 D E h 0 0 0 0 D h 0 0 -σ k C T h 0 D h     K h =     0 -D I h 0 0 K E h S h 0 C h 0 0 0 -D h 0 -c k C T h ω 2 0,k D h γ k D h    
with the additional matrices

(D I h ) i,j = Ω ϕ i (x) • ϕ j (x)dx (K E h ) i,j = Ω 1 µ b (∇ × ϕ i (x)) • ∇ × ϕ j (x) dx

B Computation of biorthogonal vector for 2-D PML

The proof of the relations given in subsection 4.6.1 is done with continuous operators M and K.

Its extension to discrete operators (i.e. matrices M h and K h ) is straightforward thanks to mass lumping. We have the eigenvalue problem:

0 0 0 -µ b   , K =          ε b σ x + σ y 2 ε b σ x -σ y 2 -div ε b σ x -σ y 2 ε b σ x + σ y 2 -div ⊥ ∇ 0 -µ b σ         
where div ⊥ v = ∂v where λ m = iω m is the eigenvalue with

M =     ε b 0 0 0 0 µ b 0 0 -1 0 1 0 0 -1 0 1     , K =     ε b T 2,3,1 0 0 -∇× 0 µ b T 2,3,1 ∇× 0 -T 1,2,3 0 T 3,1,2 0 0 -T 1,2,3 0 T 3,1,2     , U =     E H E * H *    

Inria

In order to find the left eigenvector of this system, we consider the adjoint eigenvalue problem to this system.

K T U ⊥ = λM T U ⊥ ,
with U ⊥ = E ⊥ , H ⊥ , E * ,⊥ , H * ,⊥ , which grants us the following system of equations: We are now going to try to identify the different components of U ⊥ . First off we can show E ⊥ = E * . The third equation and second equation of (36) give:

           ε b (-λ m + T 2,
E * ,⊥ = -∇ × H ⊥ -λ m + T 3,1,2 = -1 -λ m + T 3,1,2 ∇ × -λ m + T 1,2,3 µ b (-λ m + T 2,3,1 ) H * ,⊥
The first equation and fourth equation of ( 36), provide

E * ,⊥ = ε b (-λ m + T 2,3,1 ) -λ + T 1,2,3 E ⊥ , H * ,⊥ = ∇ × E ⊥ -λ m + T 3,1,2
By substituting these expressions in the previous equation, we obtain an equation in E ⊥ only:

-

  i ω µ(ω)H S + ∇ × E S = iω(µ(ω) -µ b )H inc , + Sommerfeld condition, where E inc , H inc is the incident field, and ε b , µ b the background indices. The incident fields E inc , H inc satisfy homogeneous Maxwell's equations with indices ε b , µ b . Let us introduce J = iω(ε(ω) -ε b ) E inc and we consider only non-ferromagnetic media such that µ(ω) = µ b = µ 0 in the physical domain.

  and the variational formulation satisfied by E m (with only Dirichlet or Neumann boundary conditions):

Figure 2 :

 2 Figure 2: All numerical adimensionaliazed eigenvalues for the golden disk, degenerate eigenvalues are located within red squares.

Figure 3 :

 3 Figure 3: Numerical eigenvalues located in the box [0, 2] × [-1, 0] for the golden disk. In red, analytical eigenvalues of Quasi Normal Modes.

  We consider the scattering by a silicon square[-400nm, 400nm] 2 surrounded by vacuum. The computational domain is a rectangle [-500nm, 500nm] × [0, 500nm]. A Neumann condition is imposed on the axis of symmetry y = 0. The silicon is modeled by a Lorentz model with five poles as proposed inSehmi et al. (2017a), the values for Ω k and σ k are repeated in table 3. The

Figure 6 :

 6 Figure 6: All numerical adimensionaliazed eigenvalues for the silicon square, degenerate eigenvalues are located within red squares. On the right, eigenvalues contained in the box [0, 5] × [-1, 0].

Figure 7 :

 7 Figure 7: Real part of the total field for wavelengthes λ = 200nm, 400nm, 500nm, 600nm for the silicon square (left to right and top to bottom).

Figure 8 :

 8 Figure 8: Relative L 2 error for the interior of the silicon square versus the number of modes N .

Figure 9 :Figure 10 :

 910 Figure 9: Relative L 2 error versus the pulsation ω for the interior of the silicon square

Figure 11 :Figure 13 :

 1113 Figure 11: Mesh used for the scattering of a sphere

Figure 14 :Figure 15 :

 1415 Figure 14: Numerical adimensionalized pulsations for the silicon sphere. Numerical eigenvalues are in blue, analytical QNMs in red.

Figure 16 :

 16 Figure16: Relative error on curl of E S versus the spectral width. Case of the silicon sphere, only interior of the sphere is taken into account.

Figure 17 :

 17 Figure 17: Relative error on electric field E S (as computed in (35)) versus the spectral width. Case of the silicon sphere.

Figure 19 :

 19 Figure19: Numerical eigenpulsations ω m in the complex plane for the silver rod.

Figure 20 :Figure 21 :

 2021 Figure 20: Real part of the total field (component E x of electric field) for three planes Oxy, Oxz and Oyz. Case of a silver rod (λ = 400 nm).

C

  σ x + σ y ) = 0 and we denote y asy = -div µ -1 b (-λ m + σ) -1 ∇u ⊥ m + ∇ ⊥ (σ x -σ y )u ⊥ m 2λ m -(σ x + σ y )The part of the variational formulation associated with y will provide Ω σ y does not depend on x, we have∂ ∂x 2λ m -(σ x + σ y ) + (σ x -σ y ) 2λ m -(σ x + σ y ) u ⊥ m = 2(λ m -σ y ) ∂ ∂x u ⊥ m 2λ m -σ x + σ y = (-λ m + σ y ) m -(σ x + σ y ) -(σ x -σ y ) 2λ -(σ x + σ y ) ε b (-λ m + σ x )(-λ m + σ y )u mTherefore, the unknown u m satisfies the following variational formulationΩ ε b (-λ + σ x ) (-λ + σ y ) u m ϕ + µ m • ∇ϕ dΩ = 0which is the same variational formulation satisfied by ũm (first component of the eigenvector of Kx m = λ m Mx m ). u m is proportional to ũm if λ m is a simple eigenvalue. In order to have u ⊥ m = ũm in the physical domain, we will divide by -λ m . We therefore have the following relationu ⊥ m = 1 -σ x + σ y 2λ m ũm We infer that u * ,⊥ m = σ x -σ y 2λ m ũmFor the last component v ⊥ m , it is computed from u ⊥ m and u * ,⊥ Computation of biorthogonal vector for 3-D PMLWe have the eigenvalue problem KU = λ m MU

  3,1 )E ⊥ -(-λ m + T 1,2,3 )E * ,⊥ = 0 µ b (-λ m + T 2,3,1 )H ⊥ -(-λ m + T 1,2,3 )H * ,⊥ = 0 (-λ m + T 3,1,2 )E * ,⊥ + ∇ × H ⊥ = 0 (-λ m + T 3,1,2 )H * ,⊥ -∇ × E ⊥ = 0 (36)

  ε b (-λ m + T 2,3,1 ) -λ m + T 1,2,3 E ⊥ = -1 -λ m + T 3,1,2 ∇ × --λ m + T 1,2,3 µ b (-λ m + T 2,3,1 )(-λ m + T 3,1,2 ) ∇ × E ⊥Since E * verifies the same eigenvalue problem, we can choose the constant such thatE ⊥ = E *Next, we will show thatH ⊥ = -H * . Using H * ,⊥ = ∇ × E * -λ m + T 3,1,2 and H * ,⊥ = µ b (-λ m + T 2,3,1 ) -λ m + T 1,2,3 H ⊥ ,we can show thatH ⊥ = -λ m + T 1,2,3 µ b (-λ m + T 2,3,1 )(-λ m + T 3,1,2 ) ∇ × E * = -H *and from there, it can easily be shown thatE * ,⊥ = -λ m + T 2,3,1 -λ m + T 3Therefore we have obtained the left eigenvector given in formula (33). RR n°9348 RESEARCH CENTRE BORDEAUX -SUD-OUEST 200 avenue de la Vieille Tour 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr

Table 3 :

 3 Constants ε ∞ , σ k , Ω k for silicon.

Table 5 :

 5 , we have represented the Constants ε ∞ , η k , γk , σ k , Ω k for silver. measured ε(ω) and the computed model. Fourth-order edge elements are used, the numerical Figure 18: Real and imaginary part of ε versus ω. Comparison between the measured permittivity and the Drude-Lorentz model.eigenpulsations obtained for this case are displayed in figure19. The matrices M h , K h have 36 089 rows. Among the 9836 stored eigenvectors, we have 486 eigenvectors associated with a degenerate eigenvalue. For this case, we observe in figure19that degenerate eigenvalues appear only for poles and zeros of the permittivity ε(ω). Degenerate eigenvalues associated with PMLs appear most of the time because of the symmetry of the PML layout. For example, if we replace

	ε ∞	η 1	γ1	η 2	γ2
	1.2135 883.46 0 -883.46 0.0866887
	σ 1	Ω 1		σ 2	Ω 2
	3.5366 + 4.5884i 4.249 -2.1153i 0.74248e -0.19950i 4.0097 -0.30102i

RR n°9348


domain: sdv sdv.bv sde

text:
Management of the collection

Around 27 000 accessions of major agronomical species, wheat, barley, oat, rye, triticale and wild relatives, are conserved in the small grain cereals BRC. Activities of collection's management have been certified under the French

WHAT IS BEHIND THE MANAGEMENT OF SMALL GRAIN CEREALS BIOLOGICAL RESOURCES CENTRE?

A Biological Resources Centre (BRC) gathers collections which are, in the case of agronomical BRC, most of the time patrimonial genetic resources (landraces, breeding lines and registered cultivars). Some of its activities are well known as the preservation, propagation and diffusion of genetic resources.

Strongly linked to its first mission, BRC is also involved in genetic resources valorisation. These activities are carried out in two ways: by being a member of genetic resources networks and being involved directly in scientific programs at national and European level.

Diversity's valorisation

Genetic diversity held in BRC collection is characterized and analyzed in several project to better understand its structure.

Analysis of Agronomical diversity in French Networks

  .inra.fr www.clermont.inra.fr/umr1095/ http://urgi.versailles.inra.fr/siregal

  96-900 specific of Biological Resources Centre for these activities of receipt, preparation, propagation and provision of small grain cereal biological resources.

  and IRD AFNOR quality committee To adapt NF S 96-900 at ISO level Knowledge on structure at spatial and temporal scales of collections New diversity

  Audrey Didier, Cindy Haller, Sandra Robert, Marion Deloche, Florence Exbrayat, Lionel Bardy,

				François Balfourier
	Quality system
			Regulatory
	Equipments Monitoring system Metrology	Nagoya protocol European regulation 511/2014 French law in process
	Laboratory expertise
	Germination test Relative humidity test	Database Associated data	Diffusion
			Stock quantity	Around 4 200
	Satisfaction survey		Seed quality	accessions
	Process review			provided by year
	Quality policy			with MTA
		National INRA Plant database SIREGAL
			Communications to Publics and Scientifics purposes


1

A few remarks on models and international agreements by Jean-Marc Boussard I have been recently involved in the appraisal of several documents regarding international agreements between CEDEAO's member states as well as between the CEDEAO and the European union. Authors were preoccupied with the consequences of such agreements for the agricultural sectors in their countries, especially from the point of view of the "rice economy". In fact, all of them were starting with excellent analysis of paddy production and processing in each countries, so that I took much pleasure reading them. Indeed, I learned a lot of things with these papers. At the same time, all of them were ending with considerations regarding the overall consequences of the envisaged treatises for the whole economy. These considerations were derived from either partial or general equilibrium models which, apparently, had been generously let at the disposal of the authors by some international institution, perhaps the World Bank.

I am a model addict myself, so that I certainly would not reject the idea of making use of model s in designing policy instruments or negotiating international treatises. Quite the contrary, in the course of my career, I often tried to convince decision makers not to neglect the numerous lessons one can derive from a careful use of such tools. Yet, I have been extremely shocked by the way specific models had been used in these specific cases.

By definition, a model is the translation of a theory into a mathematical equation (or a set of equations). Of course, it cannot be better than the theory, so that, it does not bring much by itself. Yet, in all sciences, building models is useful by providing numerical predictions which can be compared with observed reality in order to check the validity of the underlying theory. This is also true in economics, where theories are often qualitative and loosely specified . Indeed, it has been contended that models provide economics with the same degree of rigor as, for instance, astronomy.

In practice, however, models builder are aware of the fact that many incontrollable events are often occurring in economics, thus making this pretence unrealistic. The fact that everything depends upon everything, with hundreds of variables acting together in opposite directions implies that making "true" predictions is virtually impossible. This is a major reason not to have too much confidence into the validity of models results.

Yet, models are useful not so much by their results -always discussible -but by the fact that, in the process of building a model, the author must sharpen theory, and take a large number of decisions, each of which leading to a different version of the model -and, of course, to different results. Most of the time, these "variants", so plausible as they can be, lead to obviously fancy results. And getting such fancy results help in identifying the important details in the theory, be it the magnitude of such or such parameter, or the neglect of an aspect which was deemed unimportant at first sight, and which reveal itself crucial after all.

Thus, using a model for whatever purpose requires a long process of trial and error with different versions, ending with qualitative (rather than quantitative) conclusions on what is important and what is not… Nothing of this sort seems to have been present here. Reading theses documents, one get the impression that the authors were provided with a set of equations they did really not understand, a set of data they did not collect themselves and the exact meaning of which they did not know, and that they only push a button on their computers to get results they took at their facial value, as if revealed by some divinity… Of course, they were confident that the authorities which had supplied these instruments were perfectly honest… in such case, was the model really necessary ? This is just like having recourse to a wizard : if the wizard is honest, let him (or her) take the lead… otherwise, what could happen? Now, it turns out that in the present case, I am not sure the wizard is honest. All these discussions on liberalisation are so much loaded with ideology -not to speak of very real interests -that I would be very cautious before accepting these results. A major reason for this reluctance is that the trick has been widely used during the discussions known as the "Uruguay Round ", with huge models of the world economy promising world benefits as large as more than 600 billion $ from the agricultural liberalisation as decided by the Marrakech treatise. In practice, first, 600 billion $ is not a large sum at the world level, and, second, it is sure we never got it after more than 20 years. .. The models used in this occasion were in fact unanimous in their predictions only because they were all the same, with practically the same data, and the same equations, relying upon the same oversimplified theory of equilibrium, which says that, because of competition, marginal costs equate prices everywhere. Although reflecting a part of reality, this theory neglects important aspects, especially the fact that producers not only react to the mean level but also to the variability of prices. In addition, production levels are decided on the basis of guessed rather than actual price. For these (and a few other ) reasons, in a completely liberal economy, the variability of prices is likely to increase, and the production by cautious farmers to decrease by comparison with the situation where prices are administered, or at least, regulated. Then, agricultural prices, instead of decreasing, for the benefit of consumer, tend to increase. But this increase does not benefit farmers, because they reduced production, loosing on quantities what they gained in prices1 … This is actually what happened since the Uruguay Round on the world markets, culminating with the famous staple food crisis in 2008. Thus, the predictions of all these « Uruguay round models" were perfectly wrong, and by a large magnitude. Yet, they had been pushed on the front stage by the promoters of the liberalisation, who saw them as a way of putting a bright "scientific" paint over a fairly clumsy vessel… I have therefore a strong suspicion that the same scenario is going to be played again with these recent "APE" projects. It does not imply that the UE/CEDEAO treatise should not be signed : I have no opinion on this point, because I do not know what is precisely going on. I would only advise my African friends not to sign anything which would prevent them to make use of emergency measures in case of prices really too high or really too low -because my own models abundantly show the high probability of such events in the next few years, although I cannot say when exactly. For the rest, I am -just as some of my colleagues in INRA or CIRAD -ready to help volunteers to set up their own models… After such an exercise, conducted as defined above, I am sure they would never say the same kind of nonsense as before . .

(see for instance BOUSSARD, J.M., F. GERARD et M.G. PIKETTY (2005) : Libéraliser l'agriculture ? Théorie, modèles et réalités . Editions Quae, Paris,

2005.This book exists in English -Should agriculture be liberalized ? Science publisher, Enfield , New Hampshire-and in Spanish -? Se deberia liberalisar la agricultura mundial ?teorias,modelos y realidad, orientation graphica editoria, Buenos aires.-)