# Long time behaviour of the solution of Maxwell's equations in dissipative generalized Lorentz materials (I)
## A frequency dependent Lyapunov function approach

Maxence Cassier, Patrick Joly, Luis Alejandro Rosas Martínez

► To cite this version:

Maxence Cassier, Patrick Joly, Luis Alejandro Rosas Martínez. Long time behaviour of the solution of Maxwell's equations in dissipative generalized Lorentz materials (I) A frequency dependent Lyapunov function approach. Zeitschrift für Angewandte Mathematik und Physik, 2023, 74 (3), pp.115. 10.1007/s00033-023-01989-9. hal-03850472

HAL Id: hal-03850472

https://inria.hal.science/hal-03850472

Submitted on 13 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Long time behaviour of the solution of Maxwell's equations in dissipative generalized Lorentz materials (I) A frequency dependent Lyapunov function approach

Maxence Cassier<sup>a</sup>, Patrick Joly<sup>b</sup> and Luis Alejandro Rosas Martínez<sup>b</sup>

<sup>a</sup> Aix Marseille Univ, CNRS, Centrale Marseille, Institut Fresnel, Marseille, France

<sup>b</sup> ENSTA / POEMS<sup>1</sup>, 32 Boulevard Victor, 75015 Paris, France

(maxence.cassier@fresnel.fr, patrick.joly@inria.fr, alejandro.rosas@ensta-paris.fr)

October 19, 2022

**Abstract**

It is well-known that electromagnetic dispersive structures such as metamaterials can be modelled by generalized Drude-Lorentz models. The present paper is the first of two articles dedicated to dissipative generalized Drude-Lorentz open structures. We wish to quantify the loss in such media in terms of the long time decay rate of the electromagnetic energy for the corresponding Cauchy problem. By using an approach based on frequency dependent Lyapounov estimates, we show that this decay is polynomial in time. These results extend to an unbounded structure the ones obtained for bounded media in [18] via a quite different method based on the notion of cumulated past history and semi-group theory. A great advantage of the approach developed here is to be less abstract and directly connected to the physics of the system via energy balances.

**Keywords:** Maxwell's equations, passive electromagnetic media, dissipative generalized Lorentz models, long time electromagnetic energy decay rate, frequency-dependent Lyapunov estimates.

# 1 Introduction and motivation

The study of the long time behaviour of solutions of dispersive and dissipative models for linear wave propagation has already been extensively studied in the literature, primarily for applications in visco-elasticity and more recently in electromagnetism. The subject has recently known a regain of interest related to metamaterials. We can refer for instance, in electromagnetism, to the article [5] in which we presented a systematic construction of mathematical models compatible with physically motivated principle such as causality and passivity (see also [14, 6, 23]). The common point to all these models lies in that the constitutive laws include memory effects corresponding to time convolution nonlocal operators that induce dispersion (the velocity of waves is frequency dependent) and dissipation (the energy decay of the solution) that are in often intimately related.

For such models one of the most natural question is the study of the long time behaviour of the corre-
sponding Cauchy problem: prove that the energy of the solution tends to 0 when $t$ tends to $+\infty$ and study
the rate of decay. This is of course closely related to the theory of control and stabilization of dynamical

<footnote n="1"><sup>1</sup>POEMS (Propagation d'Ondes: Etude Mathématique et Simulation) is a mixed research team (UMR 7231) between CNRS (Centre National de la Recherche Scientifique), ENSTA Paris (Ecole Nationale Supérieure de Techniques Avancées) and INRIA (Institut National de Recherche en Informatique et en Automatique).</footnote>
systems where one commonly distinguishes the notion of exponential stability (which corresponds to an exponential decay of the energy) and polynomial stability (the energy decays as the inverse of a positive power of t).

# 1.1 Maxwell's equations in dispersive media

## 1.1.1 General features

Maxwell's equations relate the electric and magnetic inductions $D(\mathbf{x}, t)$ and $B(\mathbf{x}, t)$ ($\mathbf{x} \in \mathbb{R}^3$ and $t > 0$ are respectively the space and variables) to the electric and magnetic fields $E(\mathbf{x}, t)$ and $H(\mathbf{x}, t)$:

$$ \left\{ \begin{array}{l} \partial_t \mathbf{D} - \nabla \times \mathbf{H} = 0, \\ \partial_t \mathbf{B} + \nabla \times \mathbf{E} = 0. \end{array} \right. \qquad (1.1) $$

On the other hand, one defines the electric polarization and magnetization by

$$ \left\{ \begin{array}{ll} \mathbf{D} = \varepsilon_0 \mathbf{E} + \mathbf{P}_{\text{tot}}, & \mathbf{P}_{\text{tot}} : \text{electric polarization,} \\ \mathbf{B} = \mu_0 \mathbf{E} + \mathbf{M}_{\text{tot}}, & \mathbf{M}_{\text{tot}} : \text{magnetization.} \end{array} \right. \qquad (1.2) $$

where $\varepsilon_0 > 0$ and $\mu_0 > 0$ are the vacuum permittivity and permeability. The above equations are completed by the following non local constitutive laws (we consider the case of a homogeneous medium)

$$ \left\{ \begin{array}{l} \mathbf{P}_{\text{tot}}(\cdot, t) = \varepsilon_0 \int_0^t \chi_e(t-s) \mathbf{E}(\cdot, s) \, ds, \\ \mathbf{M}_{\text{tot}}(\cdot, t) = \mu_0 \int_0^t \chi_m(t-s) \mathbf{H}(\cdot, s) \, ds, \end{array} \right. \qquad (1.3) $$

where $\chi_e$ and $\chi_m$ are the electrical and magnetic susceptibilities of the material (convolutions products being understood in the distributional sense, see for e.g. [7, 24], for $(\chi_e, \chi_m)$ not in $L^1$).

In the Fourier-Laplace domain

$$ \mathbf{E}(\cdot, t) \quad \rightarrow \quad \hat{\mathbf{E}}(\cdot, \omega) = \int_{0}^{+\infty} \mathbf{E}(\cdot, t) e^{i\omega t} \, dt, \quad \text{Im } \omega > 0, $$

(1.2) and (1.3) reduce to

$$ \left\{ \begin{array}{l} \hat{\mathbf{D}}(\cdot, \omega) = \varepsilon(\omega) \hat{\mathbf{E}}(\cdot, \omega), \\ \hat{\mathbf{B}}(\cdot, \omega) = \mu(\omega) \hat{\mathbf{H}}(\cdot, \omega), \end{array} \right. \qquad (1.4) $$

where the complex permittivity $\varepsilon(\omega)$ and the complex permeability $\mu(\omega)$ are given in terms of the Fourier-Laplace transform of the susceptibility functions:

$$ \varepsilon(\omega) = \varepsilon_0 (1 + \hat{\chi}_e(\omega)) \quad \text{and} \quad \mu(\omega) = \mu_0 (1 + \hat{\chi}_m(\omega)). \qquad (1.5) $$

where $\varepsilon(\omega) \rightarrow \varepsilon_0$ and $\mu(\omega) \rightarrow \mu_0$ when $\omega \rightarrow \infty$ in $\mathbb{C}^+ := \{\omega \in \mathbb{C} / \text{Im } \omega > 0\}$. In other words, the material behaves as the vacuum at high frequencies. In the frequency domain, passivity, causality and the high frequency behaviour are traduced by the fact that (see [1, 5, 6, 23, 24] for more details)

$$ \omega \mapsto \omega \varepsilon(\omega) \text{ and } \omega \mapsto \omega \mu(\omega) \text{ are Herglotz functions,} \qquad (1.6) $$
that is to say analytic functions from $\mathbb{C}^+$ into its closure $\overline{\mathbb{C}}^+$. Furthermore as the susceptibilities $\chi_e$ and $\chi_m$ are real-valued functions in the time domain, the permittivity and permeability satisfy $\overline{\varepsilon}(\omega) = \varepsilon(-\overline{\omega})$ and $\overline{\mu}(\omega) = \mu(-\overline{\omega})$, $\forall \omega \in \mathbb{C}^+$.

**Remark 1.1.** [About the notion of passivity] The condition (1.6) is the condition which is most often used to define passive materials: we called it mathematical passivity in [5]. In the same article, we define the related notion physical passivity which is associated to the Cauchy problem associated to (1.1., 1.2, 1.3), seen as an evolution problem with respect to the electromagnetic field (E, H). In other words, we look at the free evolution of the system i.e in the absence of external sources. More precisely a material is physical passive if and only if the electromagnetic energy

$$ \mathcal{E}(t) \equiv \mathcal{E}(\mathbf{E}, \mathbf{H}, t) := \frac{1}{2} \left( \varepsilon_0 \int_{\mathbb{R}^3} |\mathbf{E}(\mathbf{x}, \mathbf{t})|^2 d\mathbf{x} + \mu_0 \int_{\mathbb{R}^3} |\mathbf{H}(\mathbf{x}, \mathbf{t})|^2 d\mathbf{x} \right). \quad (1.7) $$

can never exceeds its value at $t = 0$, that is to say

$$ \forall t \ge 0, \quad \mathcal{E}(t) \le \mathcal{E}(0). \quad (1.8) $$

It is emphasized in [5] that the above property does not mean that the electromagnetic energy is a decreasing function of time.

### 1.1.2 The generalized Lorentz media

In this paper, we shall concentrate of the most well-known subclass of models: the (dissipative) generalized Lorentz media. Such model will be called local because of the relationship between D and E or B and H can be written with ordinary differential equations. More precisely, these correspond to

$$ \mathbf{P}_{\text{tot}} = \varepsilon_0 \sum_{j=1}^{N_e} \Omega_{e,j}^2 \mathbf{P}_j, \quad \mathbf{M}_{\text{tot}} = \mu_0 \sum_{\ell=1}^{N_m} \Omega_{m,\ell}^2 \mathbf{M}_{\ell}, \quad (1.9) $$

where each $\mathbf{P}_j$ (resp. each $\mathbf{M}_\ell$) is related to E (resp. H) by an ordinary differential equation

$$ \left\{ \begin{array}{ll} \partial_t^2 \mathbf{P}_j + \alpha_{e,j} \partial_t \mathbf{P}_j + \omega_{e,j}^2 \mathbf{P}_j = \mathbf{E}, & 1 \le j \le N_e, \\ \partial_t^2 \mathbf{M}_\ell + \alpha_{m,\ell} \partial_t \mathbf{M}_\ell + \omega_{m,\ell}^2 \mathbf{M}_\ell = \mathbf{H}, & 1 \le \ell \le N_m, \end{array} \right. \quad (1.10) $$

completed by 0 initial conditions

$$ \left\{ \begin{array}{ll} \mathbf{P}_j(\cdot, 0) = \partial_t \mathbf{P}_j(\cdot, 0) = 0, & 1 \le j \le N_e, \\ \mathbf{M}_\ell(\cdot, 0) = \partial_t \mathbf{M}_\ell(\cdot, 0) = 0, & 1 \le \ell \le N_m. \end{array} \right. \quad (1.11) $$

In the above equations, the (real) coefficients $\Omega_{e,j}$, $\Omega_{m,\ell}$, $\omega_{e,j}$, $\omega_{m,\ell}$ are supposed to satisfy

$$ \Omega_{e,j} > 0, \quad \omega_{e,j} \ge 0, \quad 1 \le j \le N_e, \quad \Omega_{m,\ell} > 0, \quad \omega_{m,\ell} \ge 0, \quad 1 \le \ell \le N_m, \quad (1.12) $$

while for stability/dissipation issues the coefficients $(\alpha_{m,\ell}, \alpha_{m,\ell})$ must be positive

$$ \alpha_{e,j} \ge 0, \quad 1 \le j \le N_e, \quad \alpha_{m,\ell} \ge 0 \quad 1 \le \ell \le N_m. \quad (1.13) $$

Moreover, the reader will easily check that one can assume without any loss of generality that the couples $(\alpha_{e,j}, \omega_{e,j})$ (resp. $(\alpha_{m,\ell}, \omega_{m,\ell})$) are all distinct the ones from the others.
Note that (1.9, 1.10) corresponds to

$$ (a) \quad \varepsilon(\omega) = \varepsilon_0 \left( 1 - \sum_{j=1}^{N_e} \frac{\Omega_{e,j}^2}{\omega^2 + i \alpha_{e,j} \omega - \omega_{e,j}^2} \right), \quad (b) \quad \mu(\omega) = \mu_0 \left( 1 - \sum_{\ell=1}^{N_m} \frac{\Omega_{m,\ell}^2}{\omega^2 + i \alpha_{m,\ell} \omega - \omega_{m,\ell}^2} \right). \quad (1.14) $$

Straightforward calculations show that (1.9, 1.10) are equivalent to (1.3) with

$$ \chi_e = \sum_{j=1}^{N_e} \Omega_{e,j}^2 \chi_{e,j}, \quad \chi_m = \sum_{\ell=1}^{N_m} \Omega_{m,\ell}^2 \chi_{m,\ell} \quad (1.15) $$

where the expression of each $\chi_{\nu,j}$ for $\nu = e, m$ and $j \in \{1, \dots, N_\nu\}$ is given by

$$ \begin{align} (i) \quad & \chi_{\nu,j}(t) = 2 \delta_{\nu,j}^{-1} \sinh \left( \delta_{\nu,j} t/2 \right) e^{-(\alpha_{\nu,j} t/2)}, && \text{if } \alpha_{\nu,j} > 2\omega_{\nu,j}, \nonumber \\ (ii) \quad & \chi_{\nu,j}(t) = 2 \delta_{\nu,j}^{-1} \sin \left( \delta_{\nu,j} t/2 \right) e^{-(\alpha_{\nu,j} t/2)}, && \text{if } \alpha_{\nu,j} < 2\omega_{\nu,j}, \tag{1.16} \\ (iii) \quad & \chi_{\nu,j}(t) = t e^{-(\alpha_{\nu,j} t/2)}, && \text{if } \alpha_{\nu,j} = 2\omega_{\nu,j}, \nonumber \end{align} $$

where we have set $\delta_{\nu,j} = \sqrt{\alpha_{\nu,j}^2 - 4\omega_{\nu,j}^2}$ if $\alpha_{\nu,j} \ge 2\omega_{\nu,j}$ and $\delta_{\nu,j} = \sqrt{4\omega_{\nu,j}^2 - \alpha_{\nu,j}^2}$ if $\alpha_{\nu,j} < 2\omega_{\nu,j}$.

Note that each kernel $\chi_{\nu,j}$ is not monotonous with respect to time as soon as $\omega_{\nu,j} > 0$ or $\alpha_{\nu,j} > 0$ and tends to 0 when $t \to +\infty$ if (and only if) $\alpha_{\nu,j} > 0$ (see figure 1.1.2, first two pictures). As a consequence $\chi_\nu$ does not tend to 0 at infinity as soon as one of the $\alpha_{\nu,j}$ vanishes (see figure 1.1.2, third picture).

<figure><img src="image_2.png" /><figcaption>Figure 1: Kernels as functions of time. Left: $\chi_{\nu,j}$ for $\alpha_{\nu,j} > 2\omega_{\nu,j} > 0$. Center: $\chi_{\nu,j}$ for $0 < \alpha_{\nu,j} < 2\omega_{\nu,j}$. Right: $\chi_\nu$ for $N_\nu = 2$, $\alpha_{\nu,1} = 0$, $\alpha_{\nu,2} > 0$.</figcaption></figure>

## 1.2 A brief review of the literature

As said in introduction, there are already many existing results on the long time behaviour of the solution of dissipative dispersive systems. In this paragraph, we discuss in some detail some of the most significant contributions that are in close connection with the present work.

In the article [13], the authors considered a very abstract evolution model that in particular includes (1.3) with $\chi_m = 0$ and a function $\chi_e$ for which $\omega \hat{\chi}_e(\omega)$ is a Herglotz function. If one assumes that this function satisfies the additional assumption that

$$ \text{for a. e. } \omega \in \mathbb{R}, \quad \gamma_e(\omega) := \lim_{\zeta \in \mathbb{C}^+ \to \omega} \operatorname{Im} \zeta \hat{\chi}_e(\zeta) \text{ exists,} \quad (1.17) $$
(which is satisfied by most of dispersive materials in physics and in particular by generalized Lorentz materials) then the sufficient dissipation condition (6.4) of [13] is equivalent to

$$ \text{for a. e. } \omega \in \mathbb{R}, \quad \gamma_e(\omega) > 0, \quad \gamma_e^{-1} \in L_{\text{loc}}^1(\mathbb{R}). \qquad (1.18) $$

When applied to the generalized Lorentz model (see appendix A.1), namely when $\varepsilon(\omega)$ is given by (1.14) the above condition corresponds to

$$ \exists 1 \le j \le N_e \text{ such that } \omega_{e,j} = 0 \text{ and } \alpha_{e,j} > 0. \qquad (1.19) $$

Under this condition, it is proven that the electromagnetic energy (see definition 1.7) tends to 0 for any initial data $(\mathbf{E}_0, \mathbf{H}_0)$ in $L^2(\mathbb{R}^3)^3 \times L^2(\mathbb{R}^3)^3$:

$$ \forall (\mathbf{E}_0, \mathbf{H}_0) \in L^2(\mathbb{R}^3)^3 \times L^2(\mathbb{R}^3)^3, \quad \lim_{t \to +\infty} \mathcal{E}(t) = 0. \qquad (1.20) $$

This result is proven in [5] (section 4.4) in a much more pedestrian way on a toy problem corresponding the Drude model with $N_e = N_m = 1$, $\omega_{e,1} = \omega_{m,1} = 0$ and $\alpha_{e,1}, \alpha_{e,m} > 0$.

In the above references, the question of the rate of convergence to 0 of the electromagnetic energy is not discussed. This question is addressed in a series of work by S. Nicaise and her collaborator C. Pignotti [17], [18] (which generalizes [17]) (see also [19] for local dissipation models). These works consider the initial value problem in a bounded domain $\Omega \subset \mathbb{R}^3$ with perfectly conducting boundary conditions for which they prove polynomial stability in the sense mentioned below (see estimate (1.25)). The conditions for polynomial stability in [18] are two fold:

(i) The first condition is expressed in the time domain, more precisely in terms of regularity and decay properties of the kernels $\chi_e$ and $\chi_m$:

$$ \chi = \chi_e \text{ or } \chi_m \text{ satisfies } \chi \in C^2(\mathbb{R}^+), \quad \lim_{t \to +\infty} \chi'(t) = 0, \quad |\chi''(t)| \le C e^{-\delta t} \quad (\text{with } C, \delta > 0.) \qquad (1.21) $$

Note that (1.21) implies that it exists $C_1 > 0$ such that $|\chi'(t)| \le C_1 \delta^{-1} e^{-\delta t}$ and $|\chi(t)| \le C_1$ for $t \ge 0$ (see e.g. the appendix of [18] for the details). Hence, it follows (by integrations by parts of the Fourier-Laplace integral) that

$$ \omega \hat{\chi}(\omega) = i \hat{\chi}'(\omega) + i \chi(0) \quad \text{and} \quad \omega \hat{\chi}(\omega) = -(\hat{\chi}''(\omega) + \chi'(0)) \omega^{-1} + i \chi(0), \quad \forall \omega \in \mathbb{C}^+ \qquad (1.22) $$

where, as the Laplace-Fourier transform of a $L^1$ causal function, $\omega \mapsto \hat{\chi}'(\omega)$ and $\omega \mapsto \hat{\chi}''(\omega)$ are analytic on $\mathbb{C}^+$, continuous on $\overline{\mathbb{C}}^+$ and decay to 0 when $\omega \to \infty$ in $\overline{\mathbb{C}}^+$. Thus, using (1.22), one observes that (1.21) implies that the functions $\omega \mapsto \omega \hat{\chi}_e(\omega)$ and $\omega \mapsto \omega \hat{\chi}_m(\omega)$ are analytic on $\mathbb{C}^+$ and can be extended as continuous and bounded functions in the closed upper half-plane $\overline{\mathbb{C}}^+$. Furthermore, $\hat{\chi} = \hat{\chi}_e, \hat{\chi}_m$, one has $\omega \hat{\chi}(\omega) = i\chi(0) - \chi'(0)\omega^{-1} + o(\omega^{-1})$ when $\omega \to \infty$ in $\overline{\mathbb{C}}^+$.

(ii) The second condition is expressed in the frequency domain for real frequencies. It also has two parts. The first one is a strict positivity condition

$$ \forall \omega \in \mathbb{R}^*, \quad \gamma_e(\omega) := \operatorname{Im} \omega \hat{\chi}_e(\omega) > 0, \quad \gamma_m(\omega) := \operatorname{Im} \omega \hat{\chi}_m(\omega) > 0. \qquad (1.23) $$

which is completed the additional assumption

$$ \exists \omega_0, q, C > 0 \text{ such that } |\omega| \ge \omega_0 \quad \Rightarrow \quad \gamma_e(\omega) \ge C |\omega|^{-q}, \gamma_m(\omega) \ge C \omega^{-q}, \qquad (1.24) $$
that means the the (strictly positive) functions $\gamma_e(\omega)$ and $\gamma_e(\omega)$ do not decay too fast at infinity.

Under assumptions (1.21), (1.23) and (1.24), the authors of [18] prove, for $H^1$ initial data, decay estimates of the form

$$ \mathcal{E}(t) \le C_q t^{-\frac{2}{q}} (\|\mathbf{E}_0\|_{H^1(\Omega)}^2 + \|\mathbf{H}_0\|_{H^1(\Omega)}^2). \quad (1.25) $$

When specified to the case of the Lorentz kernels, it is easy to see (with (1.15) and (1.16)) that, in addition to the non-negativity of the coefficients $(\alpha_{e,j}, \alpha_{m,l})$, the conditions (1.21), (1.23) and (1.24) (which is then satisfied for $q = 2$) correspond to the following strong dissipation condition.

**Definition 1.2** (Strong Dissipation for Lorentz models).

$$ \forall 1 \le j \le N_e, \quad \alpha_{e,j} > 0, \quad \forall 1 \le l \le N_m, \quad \alpha_{m,l} > 0. \quad (1.26) $$

It is worthwhile mentioning briefly the techniques of proof for (1.20) in [13] and (1.25) in [18].

The technique used in [13] is difficult to describe in a few lines but we can give some of the main ideas. The authors use an augmented formulation of the evolution problem where, typically, the convolutions (1.3) are hidden behind the introduction of new unknowns. In the very abstract framework of [13], these new unknowns can be seen as elements of an adequately constructed Hilbert space. In the case of the Maxwell's equations, these new unknowns are functions of $(\mathbf{x}, t)$ but also of additional variable $\xi$ that varies in $\mathbb{R}$ (or a subset of $\mathbb{R}$), see [13] and also [5, 14]. The fundamental property of the obtained "augmented" system is that it is conservative: in other words, it is an evolution problem associated with a self-adjoint operator $\mathcal{A}_c$ in the augmented Hilbert space:

$$ \frac{d\mathcal{U}_c}{dt} + i \mathcal{A}_c \mathcal{U}_c = 0. $$

This allows to use tools from spectral theory of self-adjoint operators and to obtain an adequate (Fourier-like) integral representation of the solution of the original problem. The convergence result (1.20) then appears as a consequence of the spectral theorem and the use of Riemann-Lebesgue theorem on compact sets, that is justified by the assumption (1.18).

Although quite different, the approach of [18] also starts from another augmented formulation of the original system in which the convolutions (1.3) are again replaced by additional purely differential equations. The construction of this augmented model relies of the very nice trick of Dafermos [9] for treating viscoelasticity. This implies to introduce an additional time variable $s$ (that plays a similar role than $\xi$ in [5]) and additional unknowns: the so called cumulated past histories of the fields $\mathbf{E}$ and $\mathbf{H}$. The convolution operators are then replaced by non homogeneous transport equations in the $(s, t)$ plane. Contrary to [13], the augmented system is not conservative and is written as an autonomous evolution problem involving an operator $\mathcal{A}_a$ which is not self-adjoint but maximal dissipative:

$$ \frac{d\mathcal{U}_a}{dt} + i \mathcal{A}_a \mathcal{U}_a = 0 $$

This problem arises from the application of the theory of semi-groups. In particular, estimates (1.25) are obtained by the application of theorem 2.4 of [4]. Applying this theorem essentially requires to establish localization results for the spectrum of $\mathcal{A}_a$ (inside $\mathbb{C}^-$ = $\{\omega \in \mathbb{C} | \mathrm{Im}(\omega) < 0\}$) and appropriate estimates for its resolvent, using the conditions (1.22,1.23,1.24).
## 1.3 About Lyapunov techniques

To conclude this short bibliographical review, it is worthwhile mentioning that other stability results for dispersive/dissipative results have also been obtained via the Lyapunov technique: roughly speaking, the idea is to derive some differential inequality (in time) for a certain functional of the solution, namely the Lyapunov function $\mathcal{L}$, which dominates the energy (or some function of the energy). The stability estimates are then obtained from the time integration of the differential inequality.

In the context of dissipative systems with memory, this type of technique was introduced in particular to show exponential stability in the theory of linear viscoelasticity [20] (see also [21, 16] in the context of electromagnetism) and used in [16] where the authors establish polynomial stability estimates associated to various damping phenomena in electromagnetism to model, for instance, a rigid electric conductor or the ionized atmosphere (see also remark 1.3). However their technique cannot be applied to our case because it requires sign properties of the derivatives of convolution kernels, which clearly prevents from time oscillations as in (1.15). To give an idea about why this kind of assumption appears, let us come back the system (1.1, 1.2, 1.3) in the whole space $\mathbb{R}^3$, provided that the kernels $\chi_e$ and $\chi_m$ are of class $C^3$ on $\mathbb{R}^+$. Adapting the techniques of [20] to dispersive Maxwell's equations, one can show (formally) the following identity (see appendix A.2)

$$ \frac{d}{dt} \mathcal{L}(\mathbf{E}, \mathbf{H})(t) + \mathcal{D}(\mathbf{E}, \mathbf{H})(t) = 0 \qquad (1.27) $$

where the Lyapunov function $\mathcal{L}(\mathbf{E}, \mathbf{H})$ is a kind of augmented energy, namely

$$ \mathcal{L}(\mathbf{E}, \mathbf{H})(t) = \mathcal{E}(\mathbf{E}, \mathbf{H})(t) + \mathcal{E}_{ad}(\mathbf{E}, \mathbf{H})(t) \qquad (1.28) $$

with the additional energy

$$ \left| \begin{array}{l} \mathcal{E}_{ad}(\mathbf{E}, \mathbf{H})(t) = \frac{\varepsilon_0}{2} \chi_e'(t) \int_{\mathbb{R}^3} |\mathbf{E}_p(\mathbf{x}, t)|^2 \mathrm{d}\mathbf{x} + \frac{\mu_0}{2} \chi_m'(t) \int_{\mathbb{R}^3} |\mathbf{H}_p(\mathbf{x}, t)|^2 \mathrm{d}\mathbf{x} \\ \qquad - \frac{\varepsilon_0}{2} \int_0^t \chi_e''(t-s) \left( \int_{\mathbb{R}^3} |\mathbf{E}_p(\mathbf{x}, t) - \mathbf{E}_p(\mathbf{x}, s)|^2 \mathrm{d}\mathbf{x} \right) \mathrm{d}s \\ \qquad - \frac{\mu_0}{2} \int_0^t \chi_m''(t-s) \left( \int_{\mathbb{R}^3} |\mathbf{H}_p(\mathbf{x}, t) - \mathbf{H}_p(\mathbf{x}, s)|^2 \mathrm{d}\mathbf{x} \right) \mathrm{d}s. \end{array} \right. \qquad (1.29) $$

and the functional $\mathcal{D}(\mathbf{E}, \mathbf{H})$ is given by

$$ \left| \begin{array}{l} \mathcal{D}(\mathbf{E}, \mathbf{H})(t) = \varepsilon_0 \chi_e(0) \int_{\mathbb{R}^3} |\mathbf{E}(\mathbf{x}, t)|^2 \mathrm{d}\mathbf{x} + \mu_0 \chi_m(0) \int_{\mathbb{R}^3} |\mathbf{H}(\mathbf{x}, t)|^2 \mathrm{d}\mathbf{x} \\ \qquad - \frac{\varepsilon_0}{2} \chi_e''(t) \int_{\mathbb{R}^3} |\mathbf{E}_p(\mathbf{x}, t)|^2 \mathrm{d}\mathbf{x} - \frac{\mu_0}{2} \chi_m''(t) \int_{\mathbb{R}^3} |\mathbf{H}_p(\mathbf{x}, t)|^2 \mathrm{d}\mathbf{x} \\ \qquad + \frac{\varepsilon_0}{2} \int_0^t \chi_e'''(t-s) \left( \int_{\mathbb{R}^3} |\mathbf{E}_p(\mathbf{x}, t) - \mathbf{E}_p(\mathbf{x}, s)|^2 \mathrm{d}\mathbf{x} \right) \mathrm{d}s \\ \qquad + \frac{\mu_0}{2} \int_0^t \chi_m'''(t-s) \left( \int_{\mathbb{R}^3} |\mathbf{H}_p(\mathbf{x}, t) - \mathbf{H}_p(\mathbf{x}, s)|^2 \mathrm{d}\mathbf{x} \right) \mathrm{d}s. \end{array} \right. \qquad (1.30) $$

Sufficient conditions to ensure stability estimates and time decay results simply amount to check that

$$ \mathcal{L}(\mathbf{E}, \mathbf{H}) \geq \mathcal{E}(\mathbf{E}, \mathbf{H}) \quad (\iff \mathcal{E}_{ad}(\mathbf{E}, \mathbf{H}) \geq 0) \quad \text{and} \quad \mathcal{D}(\mathbf{E}, \mathbf{H}) \geq 0. \qquad (1.31) $$
It is clear on (1.29) and (1.30) that (1.31) lead to the (sufficient) conditions

$$ \chi_{\nu}(0) \ge 0, \quad \chi'_{\nu}(t) \ge 0, \quad \chi''_{\nu}(t) \le 0 \quad \text{and} \quad \chi^{(3)}_{\nu}(t) \ge 0 \quad \forall t \in \mathbb{R}^{+}, \quad \nu = e, m. \qquad (1.32) $$

Such conditions are fulfilled for instance if one has

$$ \chi_{\nu}(0) > 0, \quad -\chi''_{\nu}(t) \ge \beta_{\nu} \chi'_{\nu}(t) \quad \text{and} \quad \chi'''_{\nu}(t) \ge -\beta_{\nu} \chi''_{\nu}(t) \qquad (1.33) $$

which yields the existence of a constant $\delta > 0$ such that $\mathcal{D}(\mathbf{E}, \mathbf{H}) \ge \delta \mathcal{L}(\mathbf{E}, \mathbf{H})$ so that immediately implies the exponential stability of Maxwell's equations in the sense that

$$ \mathcal{E}(t) \le C e^{-\delta t}. \qquad (1.34) $$

An elementary example of susceptibility kernels satisfying (1.33) is given by $\chi_{\nu}(t) = 2 - e^{-t}$.

Unfortunately, the conditions (1.31) are useless for analysing the stability of Maxwell's generalized Lorentz materials given by (1.10) and (1.12), the conditions (1.32) are only satisfied by non-dissipative Drude materials for which $\chi_e(t) = \Omega_e^2 t$ and $\chi_m(t) = \Omega_m^2 t$.

**Remark 1.3.** It is worthwhile to come back here to what we said in the remark 1.1 and more precisely on the possible equivalence (generally conjectured) between the two notions of mathematical passivity and physical passivity. What we show above is that the conditions (1.32) are sufficient conditions for physical passivity. As a consequence, finding functions $\chi_e$ and $\chi_m$ satisfying (1.32) but such that the Herglotz property (1.6) would not hold, would provide a counter example to the equivalence.

The conditions (1.32) are clearly reminiscent of the notion of Bernstein functions [2, 22], i.e. positive continuous function $f: [0, \infty[ \to (0, \infty)$, $C^{\infty}$ on $(0, +\infty)$ and whose derivative $f'$ is a completely monotonous functions, which means that the sign of the successive derivatives of $f$ alternate with the order of derivation:

$$ \forall j \in \mathbb{N} \setminus \{0\} \quad (-1)^{j+1} f^{(j)}(t) \ge 0. \qquad (1.35) $$

Completely monotone functions are also characterized as the Laplace transforms of positive Borel measures on $[0, +\infty)$ (see, e.g. [2, 22]). It is known (see [15], Theorem 3.2 and corollary 3.14) that Bernstein functions satisfy the Herglotz property

$$ \omega \hat{f}(\omega) \text{ is an Herglotz function, where } \hat{f}(\omega) \text{ is the Laplace-Fourier transform of } f. \qquad (1.36) $$

However, when the alternating sign property (1.35) is only true for $j \le 3$, which corresponds to (1.32), (1.36) could a priori fail.

## 1.4 Objectives and outline of our work

We revisit in this series of two papers is to revisit the stability theory of (1.1, 1.2, 1.9, 1.10), that is to say Maxwell's equations in generalized dissipative Lorentz media. For the simplicity of exposition, we shall consider the problem posed in the whole space $\mathbb{R}^3$ that authorizes the use of the Fourier transform in space. We have a double objective

* propose new constructive proofs of stability estimates based on elementary tools that avoids any use of "black box" results of abstract mathematical theory,

* extend the existing results with less restrictive assumptions than those appearing in [18] or [19], namely to the case of the weak dissipativity assumption (1.37).
**Definition 1.4** (Weak Dissipation for Lorentz models).

$$ \sum_{j=1}^{N_e} \alpha_{e,j} + \sum_{\ell=1}^{N_m} \alpha_{m,\ell} > 0. \qquad (1.37) $$

In this first paper, we shall restrict ourselves to the strong dissipation assumption (1.26) and wish to recover in a quite explicit form the results from [18] with a technique inspired by the Lyapunov approach. Compared to more standard Lyapunov methods, we introduce frequency dependent Lyapunov functions (understand spatial frequency or wave numbers) that in particular allows to distinguish the respective roles of low and high frequencies (the role of low frequencies, that does not appear in a bounded domain as in [18] is due to the fact that we work in an unbounded domain).

In the second paper, we shall use, in addition to the Fourier transform, a spectral representation of the solution that will permit us to derive sharp asymptotic long time estimates. This approach is less tricky than the frequency dependent Lyapunov approach but is technically more involved because non self-adjoint operators have to be handled. It also has the interest to only assume the weak dissipativity condition (1.37) and to provide optimal results. We also think that, in both papers, the arguments we shall use are quite close to physical notions (plane waves, dispersion analysis, energy balance) which should make these papers more accessible to physicists.

The outline of the present paper is as follows. For pedagogical purpose and to emphasize the main ideas that guided our computations, we first consider in section 2 the case of the (single) dissipative Drude model that corresponds to the particular case of (1.10) when $N_e = N_m = 1$ and $\omega_{e,1} = \omega_{e,m} = 0$ (this corresponds to the toy problem considered in the section 4.4 of [5]). In section 3, we shall extend the technical developments of the previous section to the generalized Lorentz models (1.10), emphasizing the changes to be done in order to treat this more general model. Our stability results are compared with the ones of [18]. In section 4, we present how to apply our method to bounded domains (section 4.1) and extend our results to generalized Drude-Lorentz models (section 4.2). Finally, the appendix section A gives the proofs of technical results used through the paper.

## 2 The case of the Drude model

In this section we are interested in studying the behaviour for long times of the solutions of the electric and magnetic fields, respectively, **E** and **H**, of the Drude model whose dissipative formulation is obtained by introducing the time-derivatives of the polarization term **P** and the magnetization term **M** (where for the Drude model $\mathbf{P}_{\text{tot}} = \varepsilon_0 \Omega_e^2 \mathbf{P}$ and $\mathbf{M}_{\text{tot}} = \mu_0 \Omega_m^2 \mathbf{M}$). The unknowns of the problem are

$$ \left\{ \begin{array}{l} \mathbf{E}(\mathbf{x}, t) : \mathbb{R}^3 \times \mathbb{R}^+ \rightarrow \mathbb{R}^3, \quad \mathbf{H}(\mathbf{x}, t) : \mathbb{R}^3 \times \mathbb{R}^+ \rightarrow \mathbb{R}^3, \\ \partial_t \mathbf{P}(\mathbf{x}, t) : \mathbb{R}^3 \times \mathbb{R}^+ \rightarrow \mathbb{R}^3, \quad \partial_t \mathbf{M}(\mathbf{x}, t) : \mathbb{R}^3 \times \mathbb{R}^+ \rightarrow \mathbb{R}^3, \end{array} \right. $$

and satisfy the governing equations

$$ \left\{ \begin{array}{l} \varepsilon_0 \partial_t \mathbf{E} - \nabla \times \mathbf{H} + \varepsilon_0 \Omega_e^2 \partial_t \mathbf{P} = 0, \qquad (\mathbf{x}, t) \in \mathbb{R}^3 \times \mathbb{R}^+ \qquad (2.1a) \\ \mu_0 \partial_t \mathbf{H} + \nabla \times \mathbf{E} + \mu_0 \Omega_m^2 \partial_t \mathbf{M} = 0, \qquad (\mathbf{x}, t) \in \mathbb{R}^3 \times \mathbb{R}^+ \qquad (2.1b) \\ \partial_t^2 \mathbf{P} + \alpha_e \partial_t \mathbf{P} = \mathbf{E}, \qquad (\mathbf{x}, t) \in \mathbb{R}^3 \times \mathbb{R}^+ \qquad (2.1c) \\ \partial_t^2 \mathbf{M} + \alpha_m \partial_t \mathbf{M} = \mathbf{H}, \qquad (\mathbf{x}, t) \in \mathbb{R}^3 \times \mathbb{R}^+ \qquad (2.1d) \end{array} \right. $$
completed by initial conditions

$$ \left\{ \begin{array}{l} \mathbf{E}(\cdot, 0) = \mathbf{E}_0, \quad \mathbf{H}(\cdot, 0) = \mathbf{H}_0, \\ \partial_t \mathbf{P}(\cdot, 0) = \partial_t \mathbf{M}(\cdot, 0) = 0. \end{array} \right. \qquad (2.2) $$

In (2.1), the coefficients ($\varepsilon_0, \mu_0, \Omega_e, \Omega_m$) are strictly positive. The coefficients ($\alpha_e, \alpha_m$) are strictly positive damping coefficients.

Setting $\mathcal{H} = L^2(\mathbb{R}^3)^3 \times L^2(\mathbb{R}^3)^3 \times L^2(\mathbb{R}^3)^3 \times L^2(\mathbb{R}^3)^3$, the proposition (A.3) in the appendix (see also remark A.6) insures that for $(\mathbf{E}_0, \mathbf{H}_0) \in L^2(\mathbb{R}^3)^3 \times L^2(\mathbb{R}^3)^3$, the system admits a unique mild solution $\mathbf{U} = (\mathbf{E}, \mathbf{H}, \partial_t \mathbf{P}, \partial_t \mathbf{M})$ in $C^0(\mathbb{R}^+, \mathcal{H})$ which is a strong solution in $C^1(\mathbb{R}^+, \mathcal{H})$ as soon as $(\mathbf{E}_0, \mathbf{H}_0)$ belongs to $H^1(\mathbb{R}^3)^3 \times H^1(\mathbb{R}^3)^3$. The goal of what follows is to analyze their influence on the long time behaviour (decay) of the solution. More precisely, our goal is in particular to obtain decay rates for the standard electromagnetic energy defined as following

$$ \mathcal{E}(t) = \mathcal{E}(\mathbf{E}, \mathbf{H}, t) := \frac{1}{2} \left( \varepsilon_0 \|\mathbf{E}(\cdot, t)\|_{L^2(\mathbb{R}^3)}^2 + \mu_0 \|\mathbf{H}(\cdot, t)\|_{L^2(\mathbb{R}^3)}^2 \right). \qquad (2.3) $$

This will be done through the following (augmented) energy

$$ \mathcal{L}(t) = \mathcal{L}(\mathbf{E}, \mathbf{H}, \mathbf{P}, \mathbf{M}, t) := \mathcal{E}(\mathbf{E}, \mathbf{H}, t) + \frac{1}{2} \left( \varepsilon_0 \Omega_e^2 \|\mathbf{P}(\cdot, t)\|_{L^2(\mathbb{R}^3)}^2 + \mu_0 \Omega_m^2 \|\mathbf{M}(\cdot, t)\|_{L^2(\mathbb{R}^3)}^2 \right), \qquad (2.4) $$

that is a decreasing function of time, according to the energy identity

$$ \frac{d}{dt} \mathcal{L}(t) + \alpha_e \varepsilon_0 \Omega_e^2 \int_{\mathbb{R}^3} |\mathbf{P}(\mathbf{x}, t)|^2 d\mathbf{x} + \mu_0 \alpha_m \Omega_m^2 \int_{\mathbb{R}^3} |\mathbf{M}(\mathbf{x}, t)|^2 d\mathbf{x} = 0, \qquad (2.5) $$

that is easily demonstrated by standard arguments (see [5], section 4.4.1 and remark 2.1).

For this purpose we shall use a Lyapunov function approach, based on the use the 3D spatial Fourier transform $\mathcal{F}$ defined by:

$$ \mathcal{G}(\mathbf{k}) = \mathcal{F}(\mathbf{G})(\mathbf{k}) = \frac{1}{(2\pi)^{\frac{3}{2}}} \int_{\mathbb{R}^3} \mathbf{G}(\mathbf{x}) e^{-i\mathbf{k}\cdot\mathbf{x}} d\mathbf{x} \quad \forall \mathbf{G} \in L^1(\mathbb{R}^3)^3 \cap L^2(\mathbb{R}^3)^3, $$

where we denote $\mathbf{k} \in \mathbb{R}^3$ the dual variable of $\mathbf{x}$ (or wave vector). We recall that $\mathcal{F}$ extends by density as a unitary transformation from $L^2(\mathbb{R}^3)^3$ to $L^2(\mathbb{R}^3)^3$. We set $\mathcal{E}(\mathbf{k}, t) := \mathcal{F}\{\mathbf{E}(\cdot, t)\}(\mathbf{k})$ and analogously $\mathbf{H}, \mathbf{P}$ and $\mathbf{M}$ the Fourier transforms of $\mathbf{H}, \mathbf{P}$ and $\mathbf{M}$, respectively. Accordingly we denote $\mathcal{E}_0$ the Fourier transform of $\mathbf{E}_0$ and $\mathcal{H}_0$ the Fourier transform of $\mathbf{H}_0$. Then $(\mathcal{E}, \mathbf{H}, \mathbf{P}, \mathbf{M})$ satisfy

$$ \varepsilon_0 \partial_t \mathcal{E} - i (\mathbf{k} \times \mathcal{H}) + \varepsilon_0 \Omega_e^2 \partial_t \mathcal{P} = 0, \qquad (2.6a) $$

$$ \mu_0 \partial_t \mathcal{H} + i (\mathbf{k} \times \mathcal{E}) + \mu_0 \Omega_m^2 \partial_t \mathcal{M} = 0, \qquad (2.6b) $$

$$ \partial_t^2 \mathcal{P} + \alpha_e \partial_t \mathcal{P} = \mathcal{E}, \qquad (2.6c) $$

$$ \partial_t^2 \mathcal{M} + \alpha_m \partial_t \mathcal{M} = \mathcal{H}. \qquad (2.6d) $$

According to (2.2), the above system in complemented by

$$ \mathcal{E}(\mathbf{k}, 0) = \mathcal{E}_0(\mathbf{k}), \quad \mathcal{H}(\mathbf{k}, 0) = \mathcal{H}_0(\mathbf{k}), \quad \partial_t \mathcal{P}(\mathbf{k}, 0) = \partial_t \mathcal{M}(\mathbf{k}, 0) = 0. \qquad (2.7) $$
Let us note that by multiplying (2.6a) by $\overline{\mathbb{E}}$, using (2.6c) and taking the real part, we obtain that

$$ \frac{1}{2} \frac{d}{dt} (\varepsilon_0 |\mathbb{E}(\mathbf{k}, t)|^2 + \varepsilon_0 \Omega_e^2 |\partial_t \mathbb{P}(\mathbf{k}, t)|^2) + \alpha_e \varepsilon_0 \Omega_e^2 |\partial_t \mathbb{P}(\mathbf{k}, t)|^2 - \mathrm{Re}(i(\mathbf{k} \times \mathbb{H}(\mathbf{k}, t)) \cdot \overline{\mathbb{E}(\mathbf{k}, t)}) = 0. \quad (2.8) $$

Analogously we have from (2.6b) and (2.6d) that

$$ \frac{1}{2} \frac{d}{dt} (\mu_0 |\mathbb{H}(\mathbf{k}, t)|^2 + \mu_0 \Omega_m^2 |\partial_t \mathbb{M}(\mathbf{k}, t)|^2) + \alpha_m \mu_0 \Omega_m^2 |\partial_t \mathbb{M}(\mathbf{k}, t)|^2 + \mathrm{Re}(i(\mathbf{k} \times \mathbb{E}(\mathbf{k}, t)) \cdot \overline{\mathbb{H}(\mathbf{k}, t)}) = 0. \quad (2.9) $$

Finally by adding up (2.8) and (2.9) and using $\mathrm{Re}(i(\mathbf{k} \times \mathbb{E}(\mathbf{k}, t)) \cdot \overline{\mathbb{H}(\mathbf{k}, t)}) - \mathrm{Re}(i(\mathbf{k} \times \mathbb{H}(\mathbf{k}, t)) \cdot \overline{\mathbb{E}(\mathbf{k}, t)}) = 0$, we show that

$$ \frac{d}{dt} \mathcal{L}_{\mathbf{k}} + \mathcal{D}_{\alpha, \mathbf{k}} = 0, \quad (2.10) $$

where we have set

$$ \left\{ \begin{array}{ll} \mathcal{L}_{\mathbf{k}}(t) = \mathcal{E}_{\mathbf{k}}(t) + \mathcal{E}_{\Omega, \mathbf{k}}(t), & \text{the Lyapunov density,} \\ \\ \mathcal{E}_{\mathbf{k}}(t) := \frac{1}{2}(\varepsilon_0 |\mathbb{E}(\mathbf{k}, t)|^2 + \mu_0 |\mathbb{H}(\mathbf{k}, t)|^2), & \text{the energy density,} \\ \\ \mathcal{E}_{\Omega, \mathbf{k}}(t) := \frac{1}{2}(\varepsilon_0 \Omega_e^2 |\partial_t \mathbb{P}(\mathbf{k}, t)|^2 + \mu_0 \Omega_m^2 |\partial_t \mathbb{M}(\mathbf{k}, t)|^2), & \text{the additional energy density,} \end{array} \right. \quad (2.11) $$

and the decay density (the index $\alpha$ is here to emphasize the fact that this is the term in (2.10) which involves the damping coefficients $\alpha_e$ and $\alpha_m$)

$$ \mathcal{D}_{\alpha, \mathbf{k}}(t) := \alpha_e \varepsilon_0 \Omega_e^2 |\partial_t \mathbb{P}(\mathbf{k}, t)|^2 + \alpha_m \mu_0 \Omega_m^2 |\partial_t \mathbb{M}(\mathbf{k}, t)|^2. \quad (2.12) $$

We employ the term "density" to refer to the fact that one works at fixed $\mathbf{k}$: the electromagnetic energy $\mathcal{E}$, for instance, is obtained, via Plancherel's theorem, by integration over $\mathbf{k}$ of the energy density $\mathcal{E}_{\mathbf{k}}$

$$ \mathcal{E}(t) = \int \mathcal{E}_{\mathbf{k}}(t) d\mathbf{k}. \quad (2.13) $$

**Remark 2.1.** The reader will note that (2.5) can be recovered by integrating the above identity over $\mathbf{k} \in \mathbb{R}^3$ and applying Plancherel's theorem.

We can not exploit only (2.10) for studying the long time behaviour of $\mathcal{L}_{\mathbf{k}}(t)$ because we can not estimate $\mathcal{L}_{\mathbf{k}}(t)$ with the help of $\mathcal{D}_{\alpha, \mathbf{k}}(t)$, which does not involve $\mathbb{E}$ and $\mathbb{H}$. On the other hand, we see on (2.6) that, in order to control $\mathbb{E}$ and $\mathbb{H}$, we need the second order derivatives $\partial_t^2 \mathbb{P}$, $\partial_t^2 \mathbb{M}$. This suggests to look at the system obtained after time differentiation of (2.6). Then in the same way that we obtained (2.10), we have (with obvious notation)

$$ \frac{d}{dt} \mathcal{L}_{\mathbf{k}}^1 + \mathcal{D}_{\alpha, \mathbf{k}}^1 = 0, \quad (2.14) $$

where we have introduced the first order densities (where "first order" refers to that fact that first order derivatives of the electromagnetic field are involved)

$$ \left\{ \begin{array}{ll} \mathcal{L}_{\mathbf{k}}^{1}(t) = \mathcal{E}_{\mathbf{k}}^{1}(t) + \mathcal{E}_{\Omega, \mathbf{k}}^{1}(t), & \text{first order Lyapunov density,} \\ \\ \mathcal{E}_{\mathbf{k}}^{1}(t) := \frac{1}{2}(\varepsilon_0 |\partial_t \mathbb{E}(\mathbf{k}, t)|^2 + \mu_0 |\partial_t \mathbb{H}(\mathbf{k}, t)|^2) & \text{first order energy density,} \\ \\ \mathcal{E}_{\Omega, \mathbf{k}}^{1}(t) := \frac{1}{2}(\varepsilon_0 \Omega_e^2 |\partial_t^2 \mathbb{P}(\mathbf{k}, t)|^2 + \mu_0 \Omega_m^2 |\partial_t^2 \mathbb{M}(\mathbf{k}, t)|^2), & \text{first order add. energy density,} \end{array} \right. \quad (2.15) $$
and the first order decay density

$$ \mathcal{D}_{\alpha, k}^{1}(t) := \alpha_{e} \varepsilon_{0} \Omega_{e}^{2} |\partial_{t}^{2} \mathbb{P}(\mathbf{k}, t)|^{2} + \alpha_{m} \mu_{0} \Omega_{m}^{2} |\partial_{t}^{2} \mathbb{M}(\mathbf{k}, t)|^{2}, \qquad (2.16) $$

which precisely involves $\partial_t^2 \mathbb{P}$, $\partial_t^2 \mathbb{M}$. The next step is to combine (2.10) and (2.14). In what follows, we shall use some standard notation, to begin with

$$ \langle \mathbf{k} \rangle = (1 + |\mathbf{k}|^2)^{\frac{1}{2}} \qquad (2.17) $$

and to compare two positive functions $f(\mathbf{k}, t)$ and $g(\mathbf{k}, t)$, the notation

$$ f(\mathbf{k}, t) \le g(\mathbf{k}, t) \iff \exists C > 0 \text{ such that } \forall \mathbf{k} \in \mathbb{R}^3, \forall t \ge 0, \quad f(\mathbf{k}, t) \le C g(\mathbf{k}, t). \qquad (2.18) $$

Performing the linear combination (2.10) + $\langle \mathbf{k} \rangle^{-2}$ (2.14) we obtain (cf. remark 2.2)

$$ \frac{d}{dt} \mathcal{L}_{\mathbf{k}}^{(1)} + \mathcal{D}_{\alpha, \mathbf{k}}^{(1)} = 0, \qquad (2.19) $$

where we have introduced the first order cumulated densities for the Lyapunov function, the energy and the additional energy, (note the difference of notation between $\mathcal{L}_{\mathbf{k}}^{(1)}$ and $\mathcal{L}_{\mathbf{k}}^{1}$, etc ...)

$$ \mathcal{L}_{\mathbf{k}}^{(1)} := \mathcal{L}_{\mathbf{k}} + \langle \mathbf{k} \rangle^{-2} \mathcal{L}_{\mathbf{k}}^{1}, \quad \mathcal{E}_{\mathbf{k}}^{(1)} := \mathcal{E}_{\mathbf{k}} + \langle \mathbf{k} \rangle^{-2} \mathcal{E}_{\mathbf{k}}^{1}, \quad \mathcal{E}_{\Omega, \mathbf{k}}^{(1)} := \mathcal{E}_{\Omega, \mathbf{k}}^{1} + \langle \mathbf{k} \rangle^{-2} \mathcal{E}_{\Omega, \mathbf{k}}^{1}, \qquad (2.20) $$

and the first order cumulated decay density

$$ \mathcal{D}_{\alpha, \mathbf{k}}^{(1)}(t) := \mathcal{D}_{\alpha, \mathbf{k}}(t) + \langle \mathbf{k} \rangle^{-2} \mathcal{D}_{\alpha, \mathbf{k}}^{1}(t). \qquad (2.21) $$

**Remark 2.2.** The weight $\langle \mathbf{k} \rangle^{-2}$ is here to compensate the time derivation in the expressions of $\mathcal{L}_{\mathbf{k}}^{1}$ and $\mathcal{D}_{\alpha, \mathbf{k}}^{1}$ so that the quantities summed in the definitions (2.20) and (2.21) have the same “homogeneity”.

The key point is that (2.19) can be exploited thanks to the following lemma.

**Lemma 2.3.** Assume that $\alpha_e > 0$ and $\alpha_m > 0$. Then the following estimate holds

$$ \mathcal{L}_{\mathbf{k}}^{(1)}(t) \le \langle \mathbf{k} \rangle^2 \mathcal{D}_{\alpha, \mathbf{k}}^{(1)}(t). \qquad (2.22) $$

*Proof.* In what follows, for conciseness, we omit to mention the (implicit) dependence of various quantities with respect to $t$ and/or $\mathbf{k}$. We recall that $\mathcal{L}_{\mathbf{k}}^{(1)} = \mathcal{E}_{\mathbf{k}}^{(1)} + \mathcal{E}_{\Omega, \mathbf{k}}^{(1)}$.

Since both $\alpha_e$ and $\alpha_m$ are strictly positive, we immediately observe that we can control the cumulated energy additional $\mathcal{E}_{\Omega, \mathbf{k}}^{(1)}$ with $\mathcal{D}_{\alpha, \mathbf{k}}^{(1)}$:

$$ \mathcal{E}_{\Omega, \mathbf{k}} \le \mathcal{D}_{\alpha, \mathbf{k}} \quad \text{and} \quad \mathcal{E}_{\Omega, \mathbf{k}}^{1} \le \mathcal{D}_{\alpha, \mathbf{k}}^{1} \quad \text{which implies, by (2.20, 2.21),} \quad \mathcal{E}_{\Omega, \mathbf{k}}^{(1)} \le \mathcal{D}_{\alpha, \mathbf{k}}^{(1)}. \qquad (2.23) $$

For estimating the energy density $\mathcal{E}_{\mathbf{k}}$, it is natural to use the constitutive equations (2.6c) and (2.6d) to deduce, again because $\alpha_e$ and $\alpha_m$ are strictly positive,

$$ |\mathbb{E}|^2 \le |\partial_t \mathbb{P}|^2 + |\partial_t^2 \mathbb{P}|^2 \le \mathcal{D}_{\alpha, \mathbf{k}} + \mathcal{D}_{\alpha, \mathbf{k}}^{1}, \quad |\mathbb{H}|^2 \le |\partial_t \mathbb{M}|^2 + |\partial_t^2 \mathbb{M}|^2 \le \mathcal{D}_{\alpha, \mathbf{k}} + \mathcal{D}_{\alpha, \mathbf{k}}^{1}. \qquad (2.24) $$

Thus, by definition of $\mathcal{E}_{\mathbf{k}}$,

$$ \mathcal{E}_{\mathbf{k}} \le \mathcal{D}_{\alpha, \mathbf{k}} + \mathcal{D}_{\alpha, \mathbf{k}}^{1}. \qquad (2.25) $$

For estimating the first order energy density $\mathcal{E}_{\mathbf{k}}^{1}$, we need to bound $\partial_t \mathbb{E}$ and $\partial_t \mathbb{H}$ for which it is natural to use Maxwell equations (2.6a) and (2.6b). This is where we need to introduce $\mathbf{k}$-dependent coefficients
in our estimates. Indeed, from (2.6a) and (2.6b),

$$ egin{aligned} & |\partial_t \mathbb{E}|^2 \le |\mathbf{k} \times \mathbb{H}|^2 + |\partial_t \mathbb{P}|^2 \le |\mathbf{k}|^2 |\mathbb{H}|^2 + |\partial_t \mathbb{P}|^2, \\ & |\partial_t \mathbb{H}|^2 \le |\mathbf{k} \times \mathbb{E}|^2 + |\partial_t \mathbb{M}|^2 \le |\mathbf{k}|^2 |\mathbb{E}|^2 + |\partial_t \mathbb{M}|^2. \end{aligned} \quad (2.26) $$

Adding the two inequalities yields, by definition of $\mathcal{E}_{\mathbf{k}}^{1}$ (2.15) and $\mathcal{D}_{\alpha, \mathbf{k}}$ (2.16) (we use again $\alpha_e, \alpha_m > 0$)

$$ \mathcal{E}_{\mathbf{k}}^{1} \le |\mathbf{k}|^2 \mathcal{E}_{\mathbf{k}} + \mathcal{D}_{\alpha, \mathbf{k}} \le \langle \mathbf{k} \rangle^2 \mathcal{D}_{\alpha, \mathbf{k}} + |\mathbf{k}|^2 \mathcal{D}_{\alpha, \mathbf{k}}^{1} \quad (2.27) $$

where, for the second inequality, we have used (2.25). Finally, performing (2.25) + $\langle \mathbf{k} \rangle^{-2}$ (2.27), gives

$$ \mathcal{E}_{\mathbf{k}}^{(1)} \le \mathcal{D}_{\alpha, \mathbf{k}} + \mathcal{D}_{\alpha, \mathbf{k}}^{1} \le \langle \mathbf{k} \rangle^2 \mathcal{D}_{\alpha, \mathbf{k}}^{(1)} \quad (\text{since } 1 \le \langle \mathbf{k} \rangle^2). \quad (2.28) $$

Finally, (2.22) results from $\mathcal{L}_{\mathbf{k}}^{(1)} = \mathcal{E}_{\mathbf{k}}^{(1)} + \mathcal{E}_{\Omega, \mathbf{k}}^{(1)}$, (2.23) and (2.28). $\square$

Lemma (2.22) means the existence of a constant $\sigma = \sigma(\varepsilon_0, \mu_0, \alpha_e, \alpha_m, \Omega_e, \Omega_m) > 0$ such that

$$ \mathcal{D}_{\alpha, \mathbf{k}}^{(1)} \ge \sigma \langle \mathbf{k} \rangle^{-2} \mathcal{L}_{\mathbf{k}}^{(1)}. $$

Combined with (2.14), this gives the differential inequality

$$ \frac{d}{dt} \mathcal{L}_{\mathbf{k}}^{(1)} + \sigma \langle \mathbf{k} \rangle^{-2} \mathcal{L}_{\mathbf{k}}^{(1)} \le 0 $$

which gives by integration in time:

$$ \mathcal{L}_{\mathbf{k}}^{(1)}(t) \le \mathcal{L}_{\mathbf{k}}^{(1)}(0) e^{-\frac{\sigma t}{\langle \mathbf{k} \rangle^2}}. \quad (2.29) $$

Next, we show how to control $\mathcal{L}_{\mathbf{k}}^{(1)}(0)$ in terms on the initial data as

$$ \mathcal{L}_{\mathbf{k}}^{(1)}(0) \le |\mathbb{E}(\mathbf{k}, 0)|^2 + |\mathbb{H}(\mathbf{k}, 0)|^2. \quad (2.30) $$

Indeed, by definition and initial conditions (2.7),

$$ \mathcal{L}_{\mathbf{k}}(0) \le |\mathbb{E}(\mathbf{k}, 0)|^2 + |\mathbb{H}(\mathbf{k}, 0)|^2. \quad (2.31) $$

On the other hand, from (2.15), one has

$$ \mathcal{L}_{\mathbf{k}}^{1}(0) \le |\partial_t \mathbb{E}(\mathbf{k}, 0)|^2 + |\partial_t \mathbb{H}(\mathbf{k}, 0)|^2 + |\partial_t^2 \mathbb{P}(\mathbf{k}, 0)|^2 + |\partial_t^2 \mathbb{M}(\mathbf{k}, 0)|^2. $$

By (2.6a, 2.6b) at $t = 0$, $|\partial_t \mathbb{E}(\mathbf{k}, 0)|^2 + |\partial_t \mathbb{H}(\mathbf{k}, 0)|^2 \le |\mathbf{k}|^2 (|\mathbb{E}(\mathbf{k}, 0)|^2 + |\mathbb{H}(\mathbf{k}, 0)|^2)$ since $\partial_t \mathbb{P}(\mathbf{k}, 0)$ and $\partial_t \mathbb{M}(\mathbf{k}, 0)$ vanish, while (2.6c,2.6d) give $|\partial_t^2 \mathbb{P}(\mathbf{k}, 0)|^2 + |\partial_t^2 \mathbb{M}(\mathbf{k}, 0)|^2 \le |\partial_t \mathbb{E}(\mathbf{k}, 0)|^2 + |\partial_t \mathbb{H}(\mathbf{k}, 0)|^2$. Thus

$$ \mathcal{L}_{\mathbf{k}}^{1}(0) \le |\mathbf{k}|^2 (|\mathbb{E}(\mathbf{k}, 0)|^2 + |\mathbb{H}(\mathbf{k}, 0)|^2). \quad (2.32) $$

Finally (2.30) results from (2.31) + $\langle \mathbf{k} \rangle^{-2}$ (A.13). At last, substituting (2.30) into (2.29), as $\mathcal{L}_{\mathbf{k}} \le \mathcal{L}_{\mathbf{k}}^{(1)}$

$$ \mathcal{L}_{\mathbf{k}}(t) \le (|\mathbb{E}(\mathbf{k}, 0)|^2 + |\mathbb{H}(\mathbf{k}, 0)|^2) e^{-\frac{\sigma t}{\langle \mathbf{k} \rangle^2}}. \quad (2.33) $$

This inequality says that the Lyapunov density decays exponentially in time for each $|\mathbf{k}|$. It also says the decay rate depends on $|\mathbf{k}|$ and tends to 0 when $|\mathbf{k}|$ tends to $+\infty$. This is the reason why, when coming back to the augmented energy $\mathcal{L}(t)$ (cf. (2.4)) we shall obtain only polynomial decay. This is summarized in the main theorem of this section:
**Theorem 2.4.** For any $(\mathbf{E}_0, \mathbf{H}_0) \in L^2(\mathbb{R}^3)^3 \times L^2(\mathbb{R}^3)^3$, the total energy tends to 0 when $t$ tends to $+\infty$:

$$ \lim_{t \to +\infty} \mathcal{L}(t) = 0. \qquad (2.34) $$

If $(\mathbf{E}_0, \mathbf{H}_0) \in H^m(\mathbb{R}^3)^3 \times H^m(\mathbb{R}^3)^3$ for some integer $m > 0$, one has a polynomial decay rate

$$ \mathcal{L}(t) \lesssim \left( \|\mathbf{E}_0\|_{H^m(\mathbb{R}^3)}^2 + \|\mathbf{H}_0\|_{H^m(\mathbb{R}^3)}^2 \right) t^{-m}. \qquad (2.35) $$

*Proof.* From the respective definitions of $\mathcal{L}$ and $\mathcal{L}_k$, by Plancherel's identity and (2.33), we have

$$ \mathcal{L}(t) = \int_{\mathbb{R}^3} \mathcal{L}_k(t) \, d\mathbf{k} \lesssim \int_{\mathbb{R}^3} (|\mathbf{E}_0(\mathbf{k})|^2 + |\mathbf{H}_0(\mathbf{k})|^2) e^{-\frac{\sigma t}{(\mathbf{k})^2}} \, d\mathbf{k}. \qquad (2.36) $$

From Lebesgue's dominated convergence theorem, we first conclude that (2.34) holds for any initial data $(\mathbf{E}_0, \mathbf{H}_0) \in L^2(\mathbb{R}^3)^3 \times L^2(\mathbb{R}^3)^3$.

Next, in order to exploit the Sobolev regularity of the initial for obtaining (2.35), we rewrite (2.36) as follows (we simply make appear artificially the factor $\langle \mathbf{k} \rangle^{2m}/t^m$)

$$ \begin{aligned} \mathcal{L}(t) & \lesssim t^{-m} \int_{\mathbb{R}^3} \langle \mathbf{k} \rangle^{2m} (|\mathbf{E}_0(\mathbf{k})|^2 + |\mathbf{H}_0(\mathbf{k})|^2) (t/\langle \mathbf{k} \rangle^2)^m e^{-\frac{\sigma t}{(\mathbf{k})^2}} \, d\mathbf{k} \\ & \lesssim t^{-m} \int_{\mathbb{R}^3} \langle \mathbf{k} \rangle^{2m} (|\mathbf{E}_0(\mathbf{k})|^2 + |\mathbf{H}_0(\mathbf{k})|^2) F(t/\langle \mathbf{k} \rangle^2) \, d\mathbf{k}. \end{aligned} $$

where we have set $F_m(r) := r^m e^{-\sigma r}$, $r \ge 0$ which is clearly bounded on $\mathbb{R}^+$.

Setting $C_m := \sup_{r \ge 0} F_m(r) = (m/(\sigma e))^m$, by Fourier characterization of Sobolev norms, we get

$$ \mathcal{L}(t) \lesssim C_m t^{-m} \left( \|\mathbf{E}_0\|_{H^{m+1}(\mathbb{R}^3)}^2 + \|\mathbf{H}_0\|_{H^{m+1}(\mathbb{R}^3)}^2 \right). \quad \square $$

# 3 The case of the generalized Lorentz model

In this section, our goal is to extend the results in the latest section to the case of the (generalized) Lorentz model. The evolution (Cauchy) problem reads as follows

$$ \text{Find} \quad \left\{ \begin{array}{l} \mathbf{E}(\mathbf{x}, t) : \mathbb{R}^3 \times \mathbb{R}^+ \to \mathbb{R}^3 \\ \mathbf{P}_j(\mathbf{x}, t) : \mathbb{R}^3 \times \mathbb{R}^+ \to \mathbb{R}^3, \quad 1 \le j \le N_e, \end{array} \right. \quad \left\{ \begin{array}{l} \mathbf{H}(\mathbf{x}, t) : \mathbb{R}^3 \times \mathbb{R}^+ \to \mathbb{R}^3 \\ \mathbf{M}_\ell(\mathbf{x}, t) : \mathbb{R}^3 \times \mathbb{R}^+ \to \mathbb{R}^3, \quad 1 \le \ell \le N_m, \end{array} \right. $$

such that (for all $1 \le j \le N_e, 1 \le \ell \le N_m$)

$$ \varepsilon_0 \partial_t \mathbf{E} - \nabla \times \mathbf{H} + \varepsilon_0 \sum_{j=1}^{N_e} \Omega_{e,j}^2 \partial_t \mathbf{P}_j = 0, \qquad (\mathbf{x}, t) \in \mathbb{R}^3 \times \mathbb{R}^+, \qquad (3.1a) $$

$$ \mu_0 \partial_t \mathbf{H} + \nabla \times \mathbf{E} + \mu_0 \sum_{\ell=1}^{N_m} \Omega_{m,\ell}^2 \partial_t \mathbf{M}_\ell = 0, \qquad (\mathbf{x}, t) \in \mathbb{R}^3 \times \mathbb{R}^+, \qquad (3.1b) $$

$$ \partial_t^2 \mathbf{P}_j + \alpha_{e,j} \partial_t \mathbf{P}_j + \omega_{e,j}^2 \mathbf{P}_j = \mathbf{E}, \qquad (\mathbf{x}, t) \in \mathbb{R}^3 \times \mathbb{R}^+, \qquad (3.1c) $$

$$ \partial_t^2 \mathbf{M}_\ell + \alpha_{m,\ell} \partial_t \mathbf{M}_\ell + \omega_{m,\ell}^2 \mathbf{M}_\ell = \mathbf{H}, \qquad (\mathbf{x}, t) \in \mathbb{R}^3 \times \mathbb{R}^+, \qquad (3.1d) $$
completed by the following divergence free initial conditions

$$ \left\{ \begin{array}{l} \mathbf{E}(\cdot, 0) = \mathbf{E}_0, \quad \mathbf{H}(\cdot, 0) = \mathbf{H}_0 \quad \text{with} \quad \nabla \cdot \mathbf{E}_0 = \nabla \cdot \mathbf{H}_0 = 0, \\ \mathbf{P}(\cdot, 0) = \mathbf{M}(\cdot, 0) = \partial_t \mathbf{P}(\cdot, 0) = \partial_t \mathbf{M}(\cdot, 0) = 0. \end{array} \right. \qquad (3.2) $$

where for simplifying some formulas, we treat the $\mathbf{P}_j$ and the $\mathbf{M}_\ell$ collectively setting

$$ \mathbf{P} = (\mathbf{P}_j)_{j=1}^{N_e} \quad \text{and} \quad \mathbf{M} = (\mathbf{M}_\ell)_{\ell=1}^{N_m}. \qquad (3.3) $$

In the above equations, the coefficients $(\Omega_{e,j}, \Omega_{m,\ell}, \Omega_{e,j}, \Omega_{m,\ell})$ are supposed to satisfy (see remark 3.1)

$$ \begin{array}{l} 0 < \Omega_{e,1} \le \cdots \le \Omega_{e,N_e}, \quad 0 < \Omega_{m,1} \le \cdots \le \Omega_{m,N_m}, \\ \omega_{e,j} > 0, \quad 1 \le j \le N_e, \quad \omega_{m,\ell} > 0, \quad 1 \le \ell \le N_m. \end{array} \qquad (3.4) $$

and one can assume without any loss of generality that the couples $(\alpha_{e,j}, \omega_{e,j})$ (resp. $(\alpha_{m,\ell}, \omega_{m,\ell})$) are all distinct the ones from the others.

**Remark 3.1.** Note that one recovers the Drude model of section with $N_e = N_m = 1$ if $\omega_{e,1} = \omega_{m,1} = 0$.

Setting $\mathcal{H} = L^2(\mathbb{R}^3)^3 \times L^2(\mathbb{R}^3)^3 \times L^2(\mathbb{R}^3)^{3N_e} \times L^2(\mathbb{R}^3)^{3N_e} \times L^2(\mathbb{R}^3)^{3N_m} \times L^2(\mathbb{R}^3)^{3N_m}$, the proposition A.3 insures that for $(\mathbf{E}_0, \mathbf{H}_0) \in L^2(\mathbb{R}^3)^3$, the system (3.1) admits a unique mild solution $\mathbf{U} = (\mathbf{E}, \mathbf{H}, \partial_t \mathbf{P}, \partial_t \mathbf{M})$ in $C^0(\mathbb{R}^+, \mathcal{H})$, which is a strong solution in $C^1(\mathbb{R}^+, \mathcal{H})$ as soon as $(\mathbf{E}_0, \mathbf{H}_0) \in H^1(\mathbb{R}^3)^3 \times H^1(\mathbb{R}^3)^3$.

The equations (3.1) are completed by the initial conditions (3.2). If the initial electric and magnetic fields are divergence free, all vector fields appearing in (3.1) are divergence free at any time, see Proposition A.5 of appendix A.4.

$$ \nabla \cdot \mathbf{E}(\cdot, t) = \nabla \cdot \mathbf{H}(\cdot, t) = \nabla \cdot \mathbf{P}_j(\cdot, t) = \nabla \cdot \mathbf{M}_\ell(\cdot, t) = 0, \quad \forall t > 0, \quad \forall j, m. \qquad (3.5) $$

The equivalent for Lorentz of the identity (2.5) (for Drude) is (see [5])

$$ \frac{d}{dt} \mathcal{L}(t) + \varepsilon_0 \sum_{j=1}^{N_e} \alpha_{j,e} \Omega_{j,e}^2 \int_{\mathbb{R}^3} |\mathbf{P}_j(\mathbf{x}, t)|^2 \mathrm{d}\mathbf{x} + \mu_0 \sum_{\ell=1}^{N_m} \alpha_{m,\ell} \Omega_{m,\ell}^2 \int_{\mathbb{R}^3} |\mathbf{M}_\ell(\mathbf{x}, t)|^2 \mathrm{d}\mathbf{x} = 0. \qquad (3.6) $$

We assume the strong dissipation assumption (1.26), namely all the damping coefficients are positive: $\alpha_{e,j}, \alpha_{m,\ell} > 0$ for all $1 \le j \le N_e, 1 \le \ell \le N_m$.

Next, we observe that the space Fourier transforms $(\mathbb{E}, \mathbb{H}, \mathbf{P}_j, \mathbf{M}_\ell)$ of $(\mathbf{E}, \mathbf{H}, \mathbf{P}_j, \mathbf{M}_\ell)$ satisfy

$$ \varepsilon_0 \partial_t \mathbb{E} - i \mathbf{k} \times \mathbb{H} + \varepsilon_0 \sum_{j=1}^{N_e} \Omega_{e,j}^2 \partial_t \mathbf{P}_j = 0, \qquad (3.7a) $$

$$ \mu_0 \partial_t \mathbb{H} + i \mathbf{k} \times \mathbb{E} + \mu_0 \sum_{\ell=1}^{N_m} \Omega_{m,\ell}^2 \partial_t \mathbf{M}_\ell = 0, \qquad (3.7b) $$

$$ \partial_t^2 \mathbf{P}_j + \alpha_{e,j} \partial_t \mathbf{P}_j + \omega_{e,j}^2 \mathbf{P}_j = \mathbb{E}, \qquad (3.7c) $$

$$ \partial_t^2 \mathbf{M}_\ell + \alpha_{m,\ell} \partial_t \mathbf{M}_\ell + \omega_{m,\ell}^2 \mathbf{M}_\ell = \mathbb{H}. \qquad (3.7d) $$
According to (3.3), we shall set

$$ \mathbb{P} := (\mathbb{P}_j)_{j=1}^{N_e}, \quad \mathbb{M} := (\mathbb{M}_\ell)_{\ell=1}^{N_m}, \qquad (3.8) $$

and will also use the condensed notation

$$ |\mathbb{P}|^2 := \sum_{j=1}^{N_e} |\mathbb{P}_j|^2, \quad \text{and} \quad |\mathbb{M}|^2 := \sum_{\ell=1}^{N_e} |\mathbb{P}_\ell|^2. \qquad (3.9) $$

To study the long time behaviour of the solution of (3.7), it is natural to try to use the same approach than for the Drude model in section 2. As a matter of fact proceeding as for obtaining (2.10), we can derive from (3.7), and the identity $\mathrm{Re}(i(\mathbf{k} \times \mathbb{E}(\mathbf{k}, t)) \cdot \overline{\mathbb{H}(\mathbf{k}, t)}) - \mathrm{Re}(i(\mathbf{k} \times \mathbb{H}(\mathbf{k}, t)) \cdot \overline{\mathbb{E}(\mathbf{k}, t)}) = 0$, the relation

$$ \frac{d}{dt} \mathcal{L}_{\mathbf{k}} + \mathcal{D}_{\alpha, \mathbf{k}} = 0, \qquad (3.10) $$

where we have introduced the energy densities

$$ \left\{ \begin{aligned} \mathcal{L}_{\mathbf{k}}(t) &= \mathcal{E}_{\mathbf{k}}(t) + \mathcal{E}_{\Omega, \mathbf{k}}(t), \\ \mathcal{E}_{\mathbf{k}}(t) &:= \frac{1}{2} (\varepsilon_0 |\mathbb{E}(\mathbf{k}, t)|^2 + \mu_0 |\mathbb{H}(\mathbf{k}, t)|^2) \\ \mathcal{E}_{\Omega, \mathbf{k}}(t) &:= \frac{1}{2} \left( \varepsilon_0 \sum_{j=1}^{N_e} \Omega_{e,j}^2 |\partial_t \mathbb{P}_j(\mathbf{k}, t)|^2 + \mu_0 \sum_{\ell=1}^{N_m} \Omega_{m,\ell}^2 |\partial_t \mathbb{M}_\ell(\mathbf{k}, t)|^2 \right) \\ &\quad + \frac{1}{2} \left( \varepsilon_0 \sum_{j=1}^{N_e} \omega_{e,j}^2 \Omega_{e,j}^2 |\mathbb{P}_j(\mathbf{k}, t)|^2 + \mu_0 \sum_{\ell=1}^{N_m} \omega_{m,\ell}^2 \Omega_{m,\ell}^2 |\mathbb{M}_\ell(\mathbf{k}, t)|^2 \right), \end{aligned} \right. \qquad (3.11) $$

and the decay density

$$ \mathcal{D}_{\alpha, \mathbf{k}}(t) := \varepsilon_0 \sum_{j=1}^{N_e} \alpha_{e,j} \Omega_{e,j}^2 |\partial_t \mathbb{P}_j(\mathbf{k}, t)|^2 + \mu_0 \sum_{\ell=1}^{N_m} \alpha_{e,\ell} \Omega_{m,\ell}^2 |\partial_t \mathbb{M}_\ell(\mathbf{k}, t)|^2. \qquad (3.12) $$

The main novelty, with respect to the Drude case, is the apparition of the second line in the definition of $\mathcal{E}_{\Omega, \mathbf{k}}$, that involves the fields $\mathbb{P}$ and $\mathbb{M}$ and not only their time derivative.

Reasoning with the time derivatives of the fields as for the Drude case, we also have

$$ \frac{d}{dt} \mathcal{L}_{\mathbf{k}}^1 + \mathcal{D}_{\alpha, \mathbf{k}}^1 = 0, \qquad (3.13) $$

having defined the first order energy densities

$$ \left\{ \begin{aligned} \mathcal{L}_{\mathbf{k}}^{1}(t) &= \mathcal{E}_{\mathbf{k}}^{1}(t) + \mathcal{E}_{\Omega, \mathbf{k}}^{1}(t), \\ \mathcal{E}_{\mathbf{k}}^{1}(t) &:= \frac{1}{2} (\varepsilon_0 |\partial_t \mathbb{E}(\mathbf{k}, t)|^2 + \mu_0 |\partial_t \mathbb{H}(\mathbf{k}, t)|^2) \\ \mathcal{E}_{\Omega, \mathbf{k}}^{1}(t) &:= \frac{1}{2} \left( \varepsilon_0 \sum_{j=1}^{N_e} \Omega_{e,j}^2 |\partial_t \mathbb{P}_j(\mathbf{k}, t)|^2 + \mu_0 \sum_{\ell=1}^{N_m} \Omega_{m,\ell}^2 |\partial_t \mathbb{M}_\ell(\mathbf{k}, t)|^2 \right) \\ &\quad + \frac{1}{2} \left( \varepsilon_0 \sum_{j=1}^{N_e} \omega_{e,j}^2 \Omega_{e,j}^2 |\partial_t \mathbb{P}_j(\mathbf{k}, t)|^2 + \mu_0 \sum_{\ell=1}^{N_m} \omega_{m,\ell}^2 \Omega_{m,\ell}^2 |\partial_t \mathbb{M}_\ell(\mathbf{k}, t)|^2 \right), \end{aligned} \right. \qquad (3.14) $$
and the first order decay density

$$ \mathcal{D}_{\alpha, \mathbf{k}}^{1}(t) := \varepsilon_0 \sum_{j=1}^{N_e} \alpha_{e,j} \Omega_{e,j}^2 |\partial_t^2 \mathbb{P}_j(\mathbf{k}, t)|^2 + \mu_0 \sum_{\ell=1}^{N_m} \alpha_{e,\ell} \Omega_{m,\ell}^2 |\partial_t^2 \mathbb{M}_\ell(\mathbf{k}, t)|^2. \quad (3.15) $$

However, this time, (3.10) and (3.13) will not be sufficient to proceed as in the Drude case because we need to control the term in the second line of the definition (3.11) of $\mathcal{E}_{\Omega, \mathbf{k}}$, in other words $\mathbb{P}$ and $\mathbb{M}$. Because these fields only appear the constitutive laws (3.7c) and (3.7d), we need to adopt a different strategy with respect to section 2:

(i) This time, the constitutive laws (3.7c) and (3.7d) are used to control $\mathbb{P}$ and $\mathbb{M}$ (and no longer $\mathbb{E}$ and $\mathbb{H}$) in function of $\mathbb{E}$, $\partial_t \mathbb{P}$, $\partial_t^2 \mathbb{P}$, $\mathbb{H}$, $\partial_t \mathbb{M}$ and $\partial_t^2 \mathbb{M}$.

(ii) we then need to control $\mathbb{E}$ and $\mathbb{M}$ in another manner: this will be done by using the Maxwell's equations (3.7a,3.7b), via $\mathbf{k} \times \mathbb{E}$ (resp. $\mathbf{k} \times \mathbb{M}$) (this control will thus degenerate when $|\mathbf{k}|$ tends to 0) in function of $\partial_t \mathbb{H}$ and $\partial_t \mathbb{M}$ (resp. $\partial_t \mathbb{E}$ and $\partial_t \mathbb{P}$). This will use $|\mathbf{k} \times \mathbb{E}| = |\mathbf{k}| |\mathbb{E}|$ (resp. $|\mathbf{k} \times \mathbb{H}| = |\mathbf{k}| |\mathbb{H}|$), which is the counterpart in Fourier space of the free divergence property (3.5).

(iii) Finally, to control $\partial_t \mathbb{E}$ and $\partial_t \mathbb{H}$, the idea is to use again (3.7c) and (3.7d), but this time after time differentiation. Doing so, we control $\partial_t \mathbb{E}$ and $\partial_t \mathbb{H}$ with $\partial_t^2 \mathbb{P}$ and $\partial_t^2 \mathbb{M}$, which do appear in the definitions of $\mathcal{D}_{\alpha, \mathbf{k}}$ and $\mathcal{D}_{\alpha, \mathbf{k}}^1$, but also the third order derivatives $\partial_t^3 \mathbb{P}$ and $\partial_t^3 \mathbb{M}$.

That is why, in order to make appear a damping function containing the third order derivatives, we differentiate the equations of the problem once more in time, which leads to the identity

$$ \frac{d}{dt} \mathcal{L}_{\mathbf{k}}^2 + \mathcal{D}_{\alpha, \mathbf{k}}^2 = 0, \quad (3.16) $$

having defined the second order energy densities

$$ \left\{ \begin{aligned} \mathcal{L}_{\mathbf{k}}^2(t) &= \mathcal{E}_{\mathbf{k}}^2(t) + \mathcal{E}_{\Omega, \mathbf{k}}^2(t), \\ \mathcal{E}_{\mathbf{k}}^2(t) &:= \frac{1}{2} (\varepsilon_0 |\partial_t^2 \mathbb{E}(\mathbf{k}, t)|^2 + \mu_0 |\partial_t^2 \mathbb{H}(\mathbf{k}, t)|^2) \\ \mathcal{E}_{\Omega, \mathbf{k}}^2(t) &:= \frac{1}{2} \left( \varepsilon_0 \sum_{j=1}^{N_e} \Omega_{e,j}^2 |\partial_t^3 \mathbb{P}_j(\mathbf{k}, t)|^2 + \mu_0 \sum_{\ell=1}^{N_m} \Omega_{m,\ell}^2 |\partial_t^3 \mathbb{M}_\ell(\mathbf{k}, t)|^2 \right) \\ &\quad + \frac{1}{2} \left( \varepsilon_0 \sum_{j=1}^{N_e} \omega_{e,j}^2 \Omega_{e,j}^2 |\partial_t^2 \mathbb{P}_j(\mathbf{k}, t)|^2 + \mu_0 \sum_{\ell=1}^{N_m} \omega_{m,\ell}^2 \Omega_{m,\ell}^2 |\partial_t^2 \mathbb{M}_\ell(\mathbf{k}, t)|^2 \right), \end{aligned} \right. \quad (3.17) $$

and the second order decay density

$$ \mathcal{D}_{\alpha, \mathbf{k}}^2(t) := \varepsilon_0 \sum_{j=1}^{N_e} \alpha_{e,j} \Omega_{e,j}^2 |\partial_t^3 \mathbb{P}_j(\mathbf{k}, t)|^2 + \mu_0 \sum_{\ell=1}^{N_m} \alpha_{e,\ell} \Omega_{m,\ell}^2 |\partial_t^3 \mathbb{M}_\ell(\mathbf{k}, t)|^2. \quad (3.18) $$

Finally, proceeding as for obtaining (2.19), we deduce from equations (3.10) to (3.18) that

$$ \frac{d}{dt} \mathcal{L}_{\mathbf{k}}^{(2)} + \mathcal{D}_{\alpha, \mathbf{k}}^{(2)} = 0, \quad (3.19) $$

where we have introduced the second order cumulated energy densities

$$ \mathcal{L}_{\mathbf{k}}^{(2)} := \sum_{j=0}^{2} \langle \mathbf{k} \rangle^{-2j} \mathcal{L}_{\mathbf{k}}^{j} \equiv \mathcal{E}_{\mathbf{k}}^{(2)} + \mathcal{E}_{\Omega, \mathbf{k}}^{(2)}, \quad \mathcal{E}_{\mathbf{k}}^{(2)} := \sum_{j=0}^{2} \langle \mathbf{k} \rangle^{-2j} \mathcal{E}_{\mathbf{k}}^{j}, \quad \mathcal{E}_{\Omega, \mathbf{k}}^{(2)} := \sum_{j=0}^{2} \langle \mathbf{k} \rangle^{-2j} \mathcal{E}_{\Omega, \mathbf{k}}^{j}, \quad (3.20) $$
and the second order cumulated decay density

$$ \mathcal{D}_{\alpha, \mathbf{k}}^{(2)} := \sum_{j=0}^{2} \langle \mathbf{k} \rangle^{-2j} \mathcal{D}_{\alpha, \mathbf{k}}^{j}. \qquad (3.21) $$

In (3.20, 3.21), by convention, $\mathcal{L}_{\mathbf{k}}^0 = \mathcal{L}_{\mathbf{k}}$, etc ... The key point is that, according to the process (i)(ii)(iii) described above, we can bound $\mathcal{L}_{\mathbf{k}}^{(2)}$ in terms of $\mathcal{D}_{\alpha, \mathbf{k}}^{(2)}$:

**Lemma 3.2.** Assume that the strong dissipation assumption holds. Then, one has the following estimate

$$ \mathcal{L}_{\mathbf{k}}^{(2)}(t) \lesssim (\langle \mathbf{k} \rangle^2 + |\mathbf{k}|^{-2}) \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}(t). \qquad (3.22) $$

*Proof.* Before entering the technical details, let us first give the main ideas and steps of the proof. The goal is to control $\mathcal{H}^{(2)}$ with the help of $\mathcal{D}_{\alpha}^{(2)}$. Towards this goal, we observe that

$$ \mathcal{L}_{\mathbf{k}}^{(2)} \quad \text{is a norm in} \quad \mathcal{U} := (\mathbb{E}, \partial_t \mathbb{E}, \partial_t^2 \mathbb{E}, \mathbb{H}, \partial_t \mathbb{H}, \partial_t^2 \mathbb{H}, \mathbb{P}, \partial_t \mathbb{P}, \partial_t^2 \mathbb{P}, \partial_t^3 \mathbb{P}, \mathbb{M}, \partial_t \mathbb{M}, \partial_t^2 \mathbb{M}, \partial_t^3 \mathbb{M}), $$

$$ \mathcal{D}_{\alpha, \mathbf{k}}^{(2)} \quad \text{is a norm in} \quad \mathcal{V} := (\partial_t \mathbb{P}, \partial_t^2 \mathbb{P}, \partial_t^3 \mathbb{P}, \partial_t \mathbb{M}, \partial_t^2 \mathbb{M}, \partial_t^3 \mathbb{M}), $$

The idea that we shall develop is that, roughly speaking, $\mathcal{D}_{\alpha}^{(2)}$ is also a norm with respect to in $\mathcal{U}$ along the linear manifold (in the $\mathcal{U}$-space) generated by the equations (3.7a, 3.7b, 3.7c, 3.7d) and their time-derivatives. More precisely, we shall be able to control the terms which are missing in $\mathcal{V}$ namely $\mathbb{E}, \partial_t \mathbb{E}, \partial_t^2 \mathbb{E}, \mathbb{H}, \partial_t \mathbb{H}, \partial_t^2 \mathbb{H}, \mathbb{P}, \mathbb{M}$ with the terms appearing in $\mathcal{D}_{\alpha}^{(2)}$ namely $\partial_t \mathbb{P}, \partial_t^2 \mathbb{P}, \partial_t^3 \mathbb{P}, \partial_t \mathbb{M}, \partial_t^2 \mathbb{M}, \partial_t^3 \mathbb{M}$. This will be done in the following order

(a) $\left\{ \begin{array}{l} \text{Control } \partial_t \mathbb{E} \text{ with } \partial_t \mathbb{P}, \partial_t^2 \mathbb{P}, \partial_t^3 \mathbb{P} \text{ using } \partial_t (3.7c), \\ \text{Control } \partial_t \mathbb{H} \text{ with } \partial_t \mathbb{M}, \partial_t^2 \mathbb{M}, \partial_t^3 \mathbb{M} \text{ using } \partial_t (3.7d), \end{array} \right.$

(b) $\left\{ \begin{array}{l} \text{Control } \mathbb{E} \text{ with } \partial_t \mathbb{H}, \partial_t \mathbb{M} \text{ using (3.7b)}, \\ \text{Control } \mathbb{H} \text{ with } \partial_t \mathbb{E}, \partial_t \mathbb{P} \text{ using (3.7a)}, \end{array} \right.$

(c) $\left\{ \begin{array}{l} \text{Control } \partial_t^2 \mathbb{E} \text{ with } \partial_t \mathbb{H}, \partial_t \mathbb{P} \text{ using } \partial_t (3.7a), \\ \text{Control } \partial_t^2 \mathbb{H} \text{ with } \partial_t \mathbb{E}, \partial_t \mathbb{M} \text{ using } \partial_t (3.7b), \end{array} \right.$

(d) $\left\{ \begin{array}{l} \text{Control } \mathbb{P} \text{ with } \partial_t \mathbb{P}, \partial_t^2 \mathbb{P} \text{ and } \mathbb{E} \text{ using (3.7c)}, \\ \text{Control } \mathbb{M} \text{ with } \partial_t \mathbb{M}, \partial_t^2 \mathbb{M} \text{ and } \mathbb{H} \text{ using (3.7d)}. \end{array} \right.$

Of course the constants in the estimates issued from using Maxwell's equations (3.7a, 3.7b), that is to say the ones of steps (b) and (c), will be $\mathbf{k}$-dependent. It is worth while mentioning that these equations are used in a different manner in (b) and (c). In step (b), they are used to estimate lower order time derivatives of the field $(\mathbb{E}, \mathbb{H})$ with higher order time derivatives: as a consequence, the constants in the estimates will blow up when $\mathbf{k}$ tend to 0. At the contrary, in step (c), they are used to estimate higher order time derivatives of the field $(\mathbb{E}, \mathbb{H})$ with lower order time derivatives: as a consequence, the constants in the estimates will blow up when $|\mathbf{k}|$ tend to $+\infty$.

Let us now enter in the details of the proof.

**Step 1 : Control of the energy density $\mathcal{E}_{\mathbf{k}}^{(2)}$**. We control below each of the terms $\langle \mathbf{k} \rangle^{-2j} \mathcal{E}_{\mathbf{k}}^j$ appearing in the sum (3.20) defining $\mathcal{E}_{\mathbf{k}}^{(2)}$.
(a) Case $j = 1$: control of $\langle \mathbf{k} \rangle^{-2} \mathcal{E}_{\mathbf{k}}^1$.

Using the strong dissipation assumption (1.26), we obtain from differentiating (3.7c) with respect to time and a summation over $j$ according to (3.8, 3.9)

$$|\partial_t E|^2 \le |\partial_t P|^2 + |\partial_t^2 P|^2 + |\partial_t^3 P|^2. \qquad (3.23)$$

On the other hand, by definition of each $\mathcal{D}_{\alpha, \mathbf{k}}^j$, and making appear at the right hand sides the terms in the sum that defines $\mathcal{D}_{\alpha, \mathbf{k}}^{(2)}$, we have

$$\begin{align*}
|\partial_t^3 \mathbb{P}|^2 \le \mathcal{D}_{\alpha, \mathbf{k}}^2 &\implies \langle \mathbf{k} \rangle^{-2} |\partial_t^3 \mathbb{P}|^2 \le \langle \mathbf{k} \rangle^{-2} \mathcal{D}_{\alpha, \mathbf{k}}^2 = \langle \mathbf{k} \rangle^2 (\langle \mathbf{k} \rangle^{-4} \mathcal{D}_{\alpha, \mathbf{k}}^2) \\
|\partial_t^2 \mathbb{P}|^2 \le \mathcal{D}_{\alpha, \mathbf{k}}^1 &\implies \langle \mathbf{k} \rangle^{-2} |\partial_t^2 \mathbb{P}|^2 \le \langle \mathbf{k} \rangle^{-2} \mathcal{D}_{\alpha, \mathbf{k}}^1 \le \langle \mathbf{k} \rangle^2 (\langle \mathbf{k} \rangle^{-2} \mathcal{D}_{\alpha, \mathbf{k}}^1) \quad (\langle \mathbf{k} \rangle \ge 1) \\
|\partial_t \mathbb{P}|^2 \le \mathcal{D}_{\alpha, \mathbf{k}} &\implies \langle \mathbf{k} \rangle^{-2} |\partial_t^2 \mathbb{P}|^2 \le \langle \mathbf{k} \rangle^{-2} \mathcal{D}_{\alpha, \mathbf{k}} \le \langle \mathbf{k} \rangle^2 (\mathcal{D}_{\alpha, \mathbf{k}}) \quad (\langle \mathbf{k} \rangle \ge 1)
\end{align*}$$

After summation of the above three inequalities and by definition of $\mathcal{D}_{\alpha, \mathbf{k}}^{(2)}$, we deduce from (3.23) that

$$\langle \mathbf{k} \rangle^{-2} |\partial_t E|^2 \le \langle \mathbf{k} \rangle^2 \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}. \qquad (3.24)$$

Analogously from differentiating (3.7d) with respect to time, we get $\langle \mathbf{k} \rangle^{-2} |\partial_t \mathbb{H}|^2 \le \langle \mathbf{k} \rangle^2 \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}$ which, combined with (3.24) and the definition (3.14) of $\mathcal{E}_{\mathbf{k}}^1$, leads to

$$\langle \mathbf{k} \rangle^{-2} \mathcal{E}_{\mathbf{k}}^1 \le \langle \mathbf{k} \rangle^2 \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}. \qquad (3.25)$$

(b) Case $j = 0$: control of $\mathcal{E}_{\mathbf{k}}$.

Using the divergence free property (3.5) (i.e. $\mathbf{k} \cdot \mathbb{H} = 0$), (3.7a) and (3.24), we get

$$|\mathbf{k} \times \mathbb{H}|^2 = |\mathbf{k}|^2 |\mathbb{H}|^2 \le |\partial_t E|^2 + |\partial_t P|^2 \le \langle \mathbf{k} \rangle^4 \mathcal{D}_{\alpha, \mathbf{k}}^{(2)} + \mathcal{D}_{\alpha, \mathbf{k}} \le \langle \mathbf{k} \rangle^4 \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}$$

since $\mathcal{D}_{\alpha, \mathbf{k}} \le \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}$ and $1 \le \langle \mathbf{k} \rangle^4$. Thus $|\mathbb{H}|^2 \le |\mathbf{k}|^{-2} \langle \mathbf{k} \rangle^4 \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}$, or equivalently

$$|\mathbb{H}|^2 \le (\langle \mathbf{k} \rangle^2 + |\mathbf{k}|^{-2}) \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}. \qquad (3.26)$$

Similarly, from (3.7b) we obtain $|E|^2 \le (\langle \mathbf{k} \rangle^2 + |\mathbf{k}|^{-2}) \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}$ which, combined with (3.26) and the definition (3.11) of $\mathcal{E}_{\mathbf{k}}$, leads to

$$\mathcal{E}_{\mathbf{k}} \le (\langle \mathbf{k} \rangle^2 + |\mathbf{k}|^{-2}) \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}. \qquad (3.27)$$

(c) Case $j = 2$: control of $\langle \mathbf{k} \rangle^{-4} \mathcal{E}_{\mathbf{k}}^2$.

Differentiating (3.7a) in time and then using (3.25), we get

$$|\partial_t^2 E|^2 \le |\mathbf{k}|^2 |\partial_t \mathbb{H}|^2 + |\partial_t^2 P|^2 \le \langle \mathbf{k} \rangle^6 \mathcal{D}_{\alpha, \mathbf{k}}^{(2)} + \mathcal{D}_{\alpha, \mathbf{k}}^1$$

Thus, as $\mathcal{D}_{\alpha, \mathbf{k}}^1 \le \langle \mathbf{k} \rangle^2 \mathcal{D}_{\alpha, \mathbf{k}}^{(2)} \le \langle \mathbf{k} \rangle^6 \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}$, it follows that:

$$\langle \mathbf{k} \rangle^{-4} |\partial_t^2 E|^2 \le \langle \mathbf{k} \rangle^2 \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}. \qquad (3.28)$$
Analogously, from (3.7b), we get $\langle \mathbf{k} \rangle^{-4} |\partial_t^2 \mathbb{H}|^2 \lesssim \langle \mathbf{k} \rangle^2 \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}$ which, combined with (3.28) gives

$$ \langle \mathbf{k} \rangle^{-4} \mathcal{E}_{\mathbf{k}}^{2} \lesssim \langle \mathbf{k} \rangle^{2} \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}. \qquad (3.29) $$

Using (3.27, 3.25, 3.29) in the definition (3.20) of the cumulated energy density $\mathcal{E}_{\mathbf{k}}^{(2)}$, we get

$$ \mathcal{E}_{\mathbf{k}}^{(2)} \lesssim (\langle \mathbf{k} \rangle^2 + |\mathbf{k}|^{-2}) \mathcal{D}_{\alpha, \mathbf{k}}^{(2)} \qquad (3.30) $$

Step 2: Control of the energy density $\mathcal{E}_{\Omega, \mathbf{k}}^{(2)}$. We immediately observe from (3.11), (3.12) (def. of $\mathcal{L}_{\mathbf{k}}$ and $\mathcal{D}_{\mathbf{k}}$), (3.14), (3.15) (def. of $\mathcal{L}_{\mathbf{k}}^1$ and $\mathcal{D}_{\mathbf{k}}^1$), (3.17) and (3.18) (def. of $\mathcal{L}_{\mathbf{k}}^2$ and $\mathcal{D}_{\mathbf{k}}^2$) that

$$ \mathcal{E}_{\Omega, \mathbf{k}} \lesssim |\mathbb{P}|^2 + |\mathbb{M}|^2 + \mathcal{D}_{\alpha, \mathbf{k}}, \quad \mathcal{E}_{\Omega, \mathbf{k}}^1 \lesssim \mathcal{D}_{\alpha, \mathbf{k}} + \mathcal{D}_{\alpha, \mathbf{k}}^1, \quad \mathcal{E}_{\Omega, \mathbf{k}}^2 \lesssim \mathcal{D}_{\alpha, \mathbf{k}}^1 + \mathcal{D}_{\alpha, \mathbf{k}}^2 $$

thus, with the adequate linear combination, $\mathcal{E}_{\Omega, \mathbf{k}}^{(2)} \lesssim \mathcal{D}_{\alpha, \mathbf{k}}^{(2)} + |\mathbb{P}|^2 + |\mathbb{M}|^2 + \langle \mathbf{k} \rangle^{-2} \mathcal{D}_{\alpha, \mathbf{k}} + \langle \mathbf{k} \rangle^{-4} \mathcal{D}_{\alpha, \mathbf{k}}^1$, i.e.

$$ \mathcal{E}_{\Omega, \mathbf{k}}^{(2)} \lesssim \mathcal{D}_{\alpha, \mathbf{k}}^{(2)} + |\mathbb{P}|^2 + |\mathbb{M}|^2. \qquad (3.31) $$

It remains to control $\mathbb{P}$ and $\mathbb{M}$ which corresponds to the point (d) above. Using the constitutive equation (3.7c) and the definitions (3.12, 3.15) of $(\mathcal{D}_{\alpha, \mathbf{k}}, \mathcal{D}_{\alpha, \mathbf{k}}^1)$, together with (3.27), yields to

$$ |\mathbb{P}|^2 \lesssim |\partial_t \mathbb{P}|^2 + |\partial_t^2 \mathbb{P}|^2 + |\mathbb{E}|^2 \lesssim \mathcal{D}_{\alpha, \mathbf{k}} + \mathcal{D}_{\alpha, \mathbf{k}}^1 + (\langle \mathbf{k} \rangle^2 + |\mathbf{k}|^{-2}) \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}, $$

thus, as $\mathcal{D}_{\alpha, \mathbf{k}} \le \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}$ and $\mathcal{D}_{\alpha, \mathbf{k}}^1 \le \langle \mathbf{k} \rangle^2 \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}$, (see (3.21)),

$$ |\mathbb{P}|^2 \lesssim (\langle \mathbf{k} \rangle^2 + |\mathbf{k}|^{-2}) \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}. \qquad (3.32) $$

Likewise from (3.7d), we have

$$ |\mathbb{M}|^2 \lesssim (\langle \mathbf{k} \rangle^2 + |\mathbf{k}|^{-2}) \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}. \qquad (3.33) $$

Using (3.32) and (3.33) in (3.31), we finally obtained

$$ \mathcal{E}_{\Omega, \mathbf{k}}^{(2)} \lesssim (\langle \mathbf{k} \rangle^2 + |\mathbf{k}|^{-2}) \mathcal{D}_{\alpha, \mathbf{k}}^{(2)}. \qquad (3.34) $$

The announced estimate (3.22) follows from (3.30), (3.31) and the definition of $\mathcal{L}_{\mathbf{k}}^{(2)}$. $\square$

Proceeding as in section 2 for obtaining (2.29), we deduce that, for some constant $\sigma > 0$,

$$ \forall \mathbf{k} \neq 0, \quad \mathcal{L}_{\mathbf{k}}^{(2)}(t) \le \mathcal{L}_{\mathbf{k}}^{(2)}(0) e^{-\frac{\sigma t}{\langle \mathbf{k} \rangle^2 + |\mathbf{k}|^{-2}}}. \qquad (3.35) $$

We shall use an estimate for the initial value $\mathcal{L}_{\mathbf{k}}^{(2)}(0)$ which is the equivalent of the estimate for the Drude case. As getting this estimate is slightly more lengthy and tedious than for the Drude case, we give it in a lemma whose proof is delayed in the appendix, section A.3.

**Lemma 3.3.**

$$ \mathcal{L}_{\mathbf{k}}^{(2)}(0) \le |\mathbb{E}_0(\mathbf{k})|^2 + |\mathbb{H}_0(\mathbf{k})|^2. \qquad (3.36) $$

Therefore, using (3.36) in (3.35), and as $\mathcal{L}_{\mathbf{k}} \le \mathcal{L}_{\mathbf{k}}^{(2)}$, we have

$$ \mathcal{L}_{\mathbf{k}}(t) \le (|\mathbb{E}(\mathbf{k}, 0)|^2 + |\mathbb{H}(\mathbf{k}, 0)|^2) e^{-\frac{\sigma t}{\langle \mathbf{k} \rangle^2 + |\mathbf{k}|^{-2}}}. \qquad (3.37) $$

The exponential decay rate in (3.35) degenerates when $|\mathbf{k}| \to +\infty$, as in the Drude case, but also when $|\mathbf{k}| \to 0$ reason why the low (spatial) frequencies will need a special treatment leading to new assumptions
involving the moments of the initial data. For this reason, we introduce the spaces, for $p \in \mathbb{N}$,

$$ \left\{ \begin{array}{l} L_p^1(\mathbb{R}^d) := \{u \in L^1(\mathbb{R}^d) / (1 + |x|)^p u \in L^1(\mathbb{R}^d)\} \\ L_{p,0}^1(\mathbb{R}^d) := \{u \in L_p^1(\mathbb{R}^d) / \forall \alpha \text{ such that } |\alpha| \le p - 1, \int_{\mathbb{R}^d} x^\alpha u(x)dx = 0\} \end{array} \right. \quad (3.38) $$

where $\alpha = (\alpha_1, \alpha_2, \alpha_3)$ denotes a multi-index with “length” $|\alpha| = \alpha_1 + \alpha_2 + \alpha_3$ and

$$ \partial^{\alpha} := \partial_{1}^{\alpha_1} \partial_{2}^{\alpha_2} \partial_{1}^{\alpha_3} \quad \text{and} \quad x^{\alpha} := x_{1}^{\alpha_1} x_{2}^{\alpha_2} x_{3}^{\alpha_3} \quad \text{for any } x = (x_1, x_2, x_3) \in \mathbb{R}^3. \quad (3.39) $$

We point out that, from the definition (3.38), one has $L_0^1(\mathbb{R}^d) = L_{0,0}^1(\mathbb{R}^d) = L^1(\mathbb{R}^d)$.

We are now in position to state the main theorem which is expressed in terms of the augmented energy $\mathcal{L}(t)$ that naturally replaces the one defined by (2.4) for the Drude case, namely

$$ \left| \begin{array}{l} \mathcal{L}(t) := \mathcal{E}(\mathbf{E}, \mathbf{H}, t) + \frac{1}{2} \left( \varepsilon_0 \sum_{j=1}^{N_r} \Omega_{e,j}^2 \| \partial_t \mathbf{P}_j(\cdot, t) \|_{L^2(\mathbb{R}^3)}^2 + \mu_0 \sum_{\ell=1}^{N_m} \Omega_{m,\ell}^2 \| \partial_t \mathbf{M}_\ell(\cdot, t) \|_{L^2(\mathbb{R}^3)}^2 \right), \\ \qquad + \frac{1}{2} \left( \varepsilon_0 \sum_{j=1}^{N_e} \omega_{e,j}^2 \Omega_{e,j}^2 \| \mathbf{P}_j(\cdot, t) \|_{L^2(\mathbb{R}^3)}^2 + \mu_0 \sum_{\ell=1}^{N_m} \omega_{m,\ell}^2 \Omega_{m,\ell}^2 \| \mathbf{M}_\ell(\cdot, t) \|_{L^2(\mathbb{R}^3)}^2 \right). \end{array} \right. \quad (3.40) $$

**Theorem 3.4.** For any $(\mathbf{E}_0, \mathbf{H}_0) \in L^2(\mathbb{R}^3)^3 \times L^2(\mathbb{R}^3)^3$ satisfying the free divergence condition, the total energy tends to 0 when $t$ tends to $+\infty$:

$$ \lim_{t \to +\infty} \mathcal{L}(t) = 0. \quad (3.41) $$

Moreover if for some integers $m \ge 0$ and $p \ge 0$,

$$ (\mathbf{E}_0, \mathbf{H}_0) \in H^m(\mathbb{R}^3)^3 \times H^m(\mathbb{R}^3)^3, \quad (\mathbf{E}_0, \mathbf{H}_0) \in L_{p,0}^1(\mathbb{R}^3)^3 \times L_{p,0}^1(\mathbb{R}^3)^3, \quad (3.42) $$

one has a polynomial decay rate

$$ \mathcal{L}(t) \le \frac{C_{\mathrm{HF}}^{m}(\mathbf{E}_0, \mathbf{H}_0)}{t^m} + \frac{C_{\mathrm{LF}}^{p}(\mathbf{E}_0, \mathbf{H}_0)}{t^{p+\frac{d}{2}}} \quad (3.43) $$

where the above constants satisfy

$$ \left\{ \begin{array}{l} C_{\mathrm{HF}}^{m}(\mathbf{E}_0, \mathbf{H}_0) \lesssim \|\mathbf{E}_0\|_{H^m(\mathbb{R}^3)}^2 + \|\mathbf{H}_0\|_{H^m(\mathbb{R}^3)}^2, \\ C_{\mathrm{LF}}^{p}(\mathbf{E}_0, \mathbf{H}_0) \lesssim \sup_{|\alpha|=p} \|x^{\alpha}\mathbf{E}_0\|_{L^1} + \sup_{|\alpha|=p} \|x^{\alpha}\mathbf{H}_0\|_{L^1}. \end{array} \right. $$

*Proof.* From the respective definitions of $\mathcal{L}$ and $\mathcal{L}_k$, see (3.11) and (3.40), we can use Plancherel's identity in (3.37) to obtain

$$ \mathcal{L}(t) = \int_{\mathbb{R}^3} \mathcal{L}_k(t) \, d\mathbf{k} \lesssim \int_{\mathbb{R}^3} \left( |\mathbb{E}_0(\mathbf{k})|^2 + |\mathbb{H}_0(\mathbf{k})|^2 \right) e^{-\frac{\sigma t}{(\mathbf{k})^2 + |\mathbf{k}|^{-2}}} \, d\mathbf{k}. \quad (3.44) $$

Next, as announced before, we treat low and high (space) frequencies separately. We begin with high frequencies who can be treated as for the Drude model.
(i) If $|\mathbf{k}| \ge 1$, then $\langle \mathbf{k} \rangle^2 + |\mathbf{k}|^{-2} \le \langle \mathbf{k} \rangle^2 + 1 \le 2 \langle \mathbf{k} \rangle^2$. Thus

$$ \int_{|\mathbf{k}|\ge 1} \mathcal{L}_{\mathbf{k}}(t) \mathrm{d}\mathbf{k} \le \int_{|\mathbf{k}|\ge 1} \left( |\mathbb{E}_0(\mathbf{k})|^2 + |\mathbb{H}_0(\mathbf{k})|^2 \right) e^{-\frac{\sigma t}{2\langle \mathbf{k} \rangle^2}} \mathrm{d}\mathbf{k}. $$

Bounding, in the right hand side, the integral over $|\mathbf{k}| \ge 1$ by the one over $\mathbb{R}^3$, we can proceed as in the proof of Theorem 2.4 (with $\sigma/2$ instead of $\sigma$, cf. (2.36)) to obtain

$$ \int_{|\mathbf{k}|\ge 1} \mathcal{L}_{\mathbf{k}}(t) \mathrm{d}\mathbf{k} \le \left( \frac{2m}{\sigma e t} \right)^m \left( \|\mathbf{E}_0\|_{H^m(\mathbb{R}^3)}^2 + \|\mathbf{H}_0\|_{H^m(\mathbb{R}^3)}^2 \right). \quad (3.45) $$

(ii) If $|\mathbf{k}| \le 1$, then $\langle \mathbf{k} \rangle^2 + |\mathbf{k}|^{-2} \le 3 |\mathbf{k}|^{-2}$, (both 1 and $|\mathbf{k}|^2$ are smaller than $|\mathbf{k}|^{-2}$) thus

$$ \int_{|\mathbf{k}|\le 1} \mathcal{L}_{\mathbf{k}}(t) \mathrm{d}\mathbf{k} \le \int_{\mathbb{R}^3} \left( |\mathbb{E}_0(\mathbf{k})|^2 + |\mathbb{H}_0(\mathbf{k})|^2 \right) e^{-\frac{|\mathbf{k}|^2 \sigma t}{3}} \mathrm{d}\mathbf{k}. \quad (3.46) $$

The behaviour of the right hand side is obviously dominated by what happens when $|\mathbf{k}|$ tends to 0.

The condition $(\mathbf{E}_0, \mathbf{H}_0) \in L_{p,0}^1(\mathbb{R}^3)^3 \times L_{p,0}^1(\mathbb{R}^3)^3$ implies in particular that

$$ \forall \alpha / |\alpha| \le p-1, \quad \partial^\alpha \mathbb{E}_0(0) = \partial^\alpha \mathbb{H}_0(0). $$

Furthermore, as $(\mathbf{E}_0, \mathbf{H}_0) \in L_{p,0}^1(\mathbb{R}^3)^3$, their Fourier transform $\mathbb{E}_0$ and $\mathbb{H}_0$ are bounded functions of class $C^p$ whose partial derivatives are bounded up to the order $p$. Consequently, using a Taylor expansion at 0 truncated at order $p$, we have

$$ |\mathbb{E}_0(\mathbf{k})| \le |\mathbf{k}|^p \sup_{|\alpha|=p} \|\partial_\alpha \mathbb{E}_0\|_{L^\infty}, \quad |\mathbb{H}_0(\mathbf{k})| \le |\mathbf{k}|^p \sup_{|\alpha|=p} \|\partial_\alpha \mathbb{H}_0\|_{L^\infty} $$

which implies, using well known properties of the Fourier transform,

$$ |\mathbb{E}_0(\mathbf{k})| \le |\mathbf{k}|^p \sup_{|\alpha|=p} \| |x|^\alpha \mathbf{E}_0 \|_{L^1}, \quad |\mathbb{H}_0(\mathbf{k})| \le |\mathbf{k}|^p \sup_{|\alpha|=p} \| |x|^\alpha \mathbf{H}_0 \|_{L^1} $$

Substituting the above in (3.46) yields

$$ \int_{|\mathbf{k}|\le 1} \mathcal{L}_{\mathbf{k}}(t) \mathrm{d}\mathbf{k} \le \left( \int_{\mathbb{R}^d} |\mathbf{k}|^{2p} e^{-\frac{2|\mathbf{k}|^2 \sigma t}{3}} \mathrm{d}\mathbf{k} \right) \left( \sup_{|\alpha|=p} \|x^\alpha \mathbf{E}_0\|_{L^1} + \sup_{|\alpha|=p} \|x^\alpha \mathbf{H}_0\|_{L^1} \right). $$

With the change of variable $\xi = \sqrt{\sigma t} \mathbf{k}$, we compute that, for some constant $C(p, d) > 0$,

$$ \int_{\mathbb{R}^d} |\mathbf{k}|^{2p} e^{-\frac{2|\mathbf{k}|^2 \sigma t}{3}} \mathrm{d}\mathbf{k} = \frac{1}{(\sigma t)^{p+\frac{d}{2}}} \int_{\mathbb{R}^d} |\xi|^{2p} e^{-\frac{2|\xi|^2}{3}} \mathrm{d}\xi \equiv \frac{C(p, d)}{(\sigma t)^{p+\frac{d}{2}}}. $$

Finally, with another constant $C$ that only depends on $d, p$ and the parameters of the model, we get

$$ \int_{|\mathbf{k}|\le 1} \mathcal{L}_{\mathbf{k}}(t) \mathrm{d}\mathbf{k} \lesssim \frac{C}{t^{m+\frac{d}{2}}} \left( \sup_{|\alpha|=m} \|x^\alpha \mathbf{E}_0\|_{L^1} + \sup_{|\alpha|=m} \|x^\alpha \mathbf{H}_0\|_{L^1} \right). \quad (3.47) $$

At the end, the final estimate (3.43) is obtained by joining (3.45) and (3.47). $\square$

**Remark 3.5** (Comparison with the estimates of [18]). The reader will check that, if we drop the second term in the estimate (3.43), which is specific to the problem in the whole space, our results for $m = 1$ coincide qualitatively with the ones of [18], cf. (1.25) with $q = 2$).
# 4 Extensions

## 4.1 The problem in a bounded domain

One can consider the evolution problem associated to equations (3.1) (or (2.1)) but posed in a bounded Lipschitz domain $\Omega$ of $\mathbb{R}^3$ and completed, for instance, with perfectly conducting conditions (in some sense the “Dirichlet” problem for Maxwell’s equations)

$$ \mathbf{E} \times \mathbf{n} = 0 \quad \text{on } \partial\Omega, \quad (\text{with } \mathbf{n} \text{ the unit normal vector to } \partial\Omega), \qquad (4.1) $$

where $\partial\Omega$ denotes the boundary of $\Omega$ and $\mathbf{n}$ the unit outward normal vector of $\partial\Omega$. In such a case, the analysis of sections 2 and 3 can be extended. The main difference is the use of the Fourier transform in space has to be replaced by an adequate modal expansion. More precisely, we introduce the cavity eigenvalue problem:

$$ \left\{ \begin{array}{l} \text{Find } k \in \mathbb{R} \text{ and } (\mathbf{u}, \mathbf{v}) \neq 0 \in L^2(\Omega)^3 \times L^2(\Omega)^3 \text{ such that} \\ \text{i} \left( \begin{array}{cc} 0 & \nabla \times \\ -\nabla \times & 0 \end{array} \right) \left( \begin{array}{c} \mathbf{u} \\ \mathbf{v} \end{array} \right) = k \left( \begin{array}{c} \mathbf{u} \\ \mathbf{v} \end{array} \right) \quad \text{in } \Omega, \\ \nabla \cdot \mathbf{u} = \nabla \cdot \mathbf{v} = 0 \quad \text{in } \Omega, \quad \mathbf{u} \times \mathbf{n} = \mathbf{v} \cdot \mathbf{n} = 0 \quad \text{on } \partial\Omega, \end{array} \right. \qquad (4.2) $$

which corresponds to find the eigenvalues of the self-adjoint Maxwell operator $\mathcal{A}$ in the closed subspace of $L^2(\Omega)^3$, $\mathcal{H} := \{(\mathbf{u}, \mathbf{v}) \in L^2(\Omega)^3 / \nabla \cdot \mathbf{u} = \nabla \cdot \mathbf{v} = 0 \text{ and } \mathbf{v} \cdot \mathbf{n} = 0 \text{ on } \partial\Omega\}$, namely

$$ \mathcal{A} \begin{pmatrix} \mathbf{u} \\ \mathbf{v} \end{pmatrix} = i \begin{pmatrix} \nabla \times \mathbf{v} \\ -\nabla \times \mathbf{u} \end{pmatrix}, \quad \forall (\mathbf{u}, \mathbf{v}) \in D(\mathcal{A}) := \mathcal{H} \cap (H_0(\text{rot}; \Omega) \times H(\text{rot}; \Omega)), \qquad (4.3) $$

where $H(\text{rot}; \Omega) := \{\mathbf{v} \in L^2(\Omega)^3 | \nabla \times \mathbf{v} \in L^2(\Omega)^3\}$ and $H_0(\text{rot}; \Omega) := \{\mathbf{u} \in H(\text{rot}; \Omega) | \mathbf{u} \times \mathbf{n} = 0 \text{ on } \partial\Omega\}$. One show, see e.g. [11] chapter IX, that the operator $\mathcal{A}$ has a compact resolvent. From the theory of selfadjoint operators with compact resolvent, and using the symmetries of Maxwell’s equations, one knows that there is a countable infinity of cavity modes

$$ (\pm k_n, \mathbf{u}_n^{\pm}, \mathbf{v}_n^{\pm}) \in \mathbb{R} \times L^2(\Omega)^3 \times L^2(\Omega)^3, \quad n \in \mathbb{N}, \quad \text{with } k_n \ge 0, \quad k_n \to +\infty \qquad (4.4) $$

with $\mathbf{u}_n^+ = \mathbf{u}_n^-$ and $\mathbf{v}_n^+ = -\mathbf{v}_n^-$ in such a way that $\{(\mathbf{u}_n^+, \mathbf{v}_n^+), n \in \mathbb{N}\}$ form an orthonormal basis of the Hilbert space $\mathcal{H}$. Then, one can decompose the electromagnetic field as

$$ \mathbf{E}(\cdot, t) = \sum_{\pm} \sum_{n=0}^{+\infty} \mathcal{E}_{n}^{\pm}(t) \mathbf{u}_{n}^{\pm}, \quad \mathbf{H}(\cdot, t) = \sum_{\pm} \sum_{n=0}^{+\infty} \mathbb{H}_{n}^{\pm}(t) \mathbf{v}_{n}^{\pm}, $$

and the auxiliary field $\mathbf{P}_j$ and $\mathbf{M}_\ell$ (understood as in (3.1)) accordingly

$$ \mathbf{P}_j(\cdot, t) = \sum_{\pm} \sum_{n=0}^{+\infty} \mathbb{P}_{j,n}^{\pm}(t) \mathbf{u}_n^{\pm}, \quad \mathbf{M}_{\ell}(\cdot, t) = \sum_{\pm} \sum_{n=0}^{+\infty} \mathbb{M}_{\ell,n}^{\pm}(t) \mathbf{v}_n^{\pm}. $$

The rest of the analysis follows exactly the same lines as for $\Omega = \mathbb{R}^3$, modulo the following substitutions

$$ \mathbf{k} \in \mathbb{R}^3 \rightarrow \{(\pm, n), n \in \mathbb{N}\}, \quad \int_{\mathbb{R}^3} d\mathbf{k} \rightarrow \sum_{\pm} \sum_{n=0}^{+\infty}, \quad \dots \text{ etc.} $$
The results are then similar to the one of Theorem 3.4 without the second term in the right hand side of the estimate (3.43) and the hypothesis of the second line of (3.42), provided some modifications on the assumptions for the initial data which must now satisfy

$$ (\mathbf{E}_0, \mathbf{H}_0) \in \mathcal{H}_0 := (\mathrm{Ker} \, \mathcal{A})^{\perp}, \quad (\mathrm{orthogonality} \, \mathrm{in} \, \mathcal{H}) \quad (\mathrm{see} \, \mathrm{also} \, \mathrm{remark} \, 4.1) \qquad (4.5) $$

and the Sobolev regularity (3.42)(first line) must be replaced by

$$ (\mathbf{E}_0, \mathbf{H}_0) \in D(\mathcal{A}^m) \cap \mathcal{H}_0, \quad (\mathrm{see} \, \mathrm{also} \, \mathrm{remark} \, 4.2). \qquad (4.6) $$

**Remark 4.1.** The (finitely dimensional) space Ker $\mathcal{A}$ is the space of electromagnetic static fields

$$ \mathrm{Ker} \, \mathcal{A} = \{ (\mathbf{u}, \mathbf{v}) \in L^2(\Omega)^3 / \nabla \cdot \mathbf{u} = \nabla \cdot \mathbf{v} = 0, \nabla \times \mathbf{u} = \nabla \times \mathbf{v} = 0 \text{ on } \Omega \text{ and } \mathbf{u} \times \mathbf{n} = 0, \mathbf{v} \cdot \mathbf{n} = 0 \text{ on } \partial\Omega \}. $$

When $\Omega$ in simply connected, it is known that Ker $\mathcal{A} = \{0\}$ and that its dimension increases with the complexity of the topology of $\Omega$. For the proof of these assertions, we refer to [11], chapter IX. In some sense, the condition (4.5) can be seen as a substitute to the condition (3.42)(second line).

**Remark 4.2.** The condition $(\mathbf{E}_0, \mathbf{H}_0) \in D(\mathcal{A}^m)$ implies $(\mathbf{E}_0, \mathbf{H}_0) \in H_{\mathrm{loc}}^m(\Omega)^3 \times H_{\mathrm{loc}}^m(\Omega)^3$. Furthermore, for $C^\infty$ domains $\Omega$, the condition $(\mathbf{E}_0, \mathbf{H}_0) \in D(\mathcal{A}^m)$ (see e.g. [11], chapter IX) is equivalent to $(\mathbf{E}_0, \mathbf{H}_0) \in H^m(\Omega)^3 \times H^m(\Omega)^3$ for any integer $m > 0$.

## 4.2 The case of mixed Drude-Lorentz models

In the sums (1.14) defining $\varepsilon(\omega)$ and $\mu(\omega)$ the resonance frequency $\omega_{e,j}$ or $\omega_{m,l}$ are either strictly positive or zero. For instance, for $\varepsilon(\omega)$, we shall say that

$$ \frac{\Omega_{e,j}^2}{\omega^2 + i \alpha_{e,j} \omega - \omega_{e,j}^2}, \omega_j > 0 \text{ is a Lorentz term,} \qquad \frac{\Omega_{e,j}^2}{\omega^2 + i \alpha_{e,j} \omega} \text{ is a Drude term.} \qquad (4.7) $$

In section 2 (standard Drude model), we consider the case where $\varepsilon(\omega)$ and $\mu(\omega)$ contained a single Drude term while in section 3, we consider the case where $\varepsilon(\omega)$ and $\mu(\omega)$ only contained Lorentz terms, because of assumption (3.4, second line).

It is natural to look at the cases where $\varepsilon(\omega)$ (resp. $\mu(\omega)$) contains Lorentz terms but also Drude terms whose number is $N_{d,e}$ (respectively $N_{d,m}$). It appears that our Lyapunov approach can easily handle these cases (modulo minor additional manipulations) with the following results:

* If $N_{d,e} > 0$ and $N_{d,m} > 0$, the result is the same as for the Drude model (see theorem 2.4).

* If $N_{d,e} > 0$ and $N_{d,m} = 0$ (or the contrary), the result is the one for Lorentz (see theorem 3.4).

# A Appendix

## A.1 On the dissipation condition of [13] for Lorentz models

Let us recall that, when the limit (1.17) exist almost everywhere on the real axis (which is the case for Lorentz models), the sufficient dissipation condition (6.4) given in [13] reads

$$ \text{for a. e. } \omega \in \mathbb{R}, \quad \mathrm{Im} \, \omega \hat{\chi}_e(\omega) \geq \gamma(\omega)^{-1} > 0, \quad \gamma \in L_{\mathrm{loc}}^1(\mathbb{R}). \qquad (A.1) $$
By virtue of the expression (1.14)(a) of the complex permittivity $\varepsilon(\omega)$ and the formula (1.5), we have, as soon as $\omega \in \mathbb{R}$ and $\omega^2 + i\alpha_{e,j}\omega - \omega_{e,j}^2 \neq 0$ for all $j \in \{1, \dots, N_e\}$:

$$ \mathrm{Im} \, \omega \hat{\chi}_e(\omega) = \sum_{j=1}^{N_e} \frac{\alpha_{e,j} \Omega_{e,j}^2 \omega^2}{|\omega^2 + i \alpha_{e,j} \omega - \omega_{e,j}^2|^2}. $$

Let $J_0 := \{j \in \{1, \dots, N_e\} / \omega_{e,j} = 0\}$ and $J_+ := \{j \in \{1, \dots, N_e\} / \alpha_{e,j} > 0\}$.

(i) $J_+ = \emptyset$, i. e. for all $j \in \{1, \dots, N_e\}$, $\alpha_{e,j} = 0$ then $\mathrm{Im} \, \omega \hat{\chi}_e(\omega) = 0$ and (1.18)(i) can not hold.

(ii) $J_+ \neq \emptyset$. We distinguish two subcases

(a) $J_+ \cap J_0 \neq \emptyset$. This means that one $\omega_{e,j}$ vanishes, for instance $\omega_{e,1} = 0$, and $\alpha_{e,1} > 0$. Then

$$ \mathrm{Im} \, \omega \hat{\chi}_e(\omega) = \frac{\alpha_{e,1} \Omega_{e,1}^2 \omega^2}{|\omega^2 + i \alpha_{e,j} \omega|^2} + \sum_{j=2}^{N_e} \frac{\alpha_{e,j} \Omega_{e,j}^2 \omega^2}{|\omega^2 + i \alpha_{e,j} \omega - \omega_{e,j}^2|^2} \geq \frac{\alpha_{e,1} \Omega_{e,1}^2}{|\omega + i \alpha_{e,1}|} $$

in which case (1.18)(ii) holds true with $\gamma : \omega \mapsto (\alpha_{e,1} \Omega_{e,1}^2)^{-1} |\omega + i\alpha_{e,1}| \in L_{\mathrm{loc}}^1(\mathbb{R})$.

(b) $J_+ \cap J_0 = \emptyset$. In this case $\mathrm{Im} \, \omega \hat{\chi}_e(\omega) \sim_{\omega \to 0} \left( \sum_{j \in J_+} \frac{\alpha_{e,j} \Omega_{e,j}^2}{\omega_{e,j}^4} \right) \omega^2$, and (1.18)(i) can not hold.

## A.2 On the energy identity (1.27)

Let $\mathbf{E}, \mathbf{H} \in C^0(\mathbb{R}^+, H(\mathrm{rot}; \mathbb{R}^3)) \cap C^1(\mathbb{R}^+, L^2(\mathbb{R}^3)^3)$ and $\mathbf{E}, \mathbf{H}, \mathbf{D}, \mathbf{B} \in C^1(\mathbb{R}^+, L^2(\mathbb{R}^3)^3)$ be solutions of the equations (1.1), (1.2) and (1.3) where $\chi_\nu \in C^3(\mathbb{R}^+)$ for $\nu = e, m$. We explain in this appendix how to obtain the energy identity (1.27) with (1.29) and (1.30).

From equations (1.1), (1.2) and (1.3) it is straightforward to deduce the identity (with $\mathcal{E}(t)$ the electro-magnetic energy, see (1.7))

$$ \frac{d\mathcal{E}(t)}{dt} + \mathcal{I}(t) = 0, \quad \mathcal{I}(t) := \int_{\mathbb{R}^3} \partial_t \mathbf{P}_{\mathrm{tot}}(\mathbf{x}, t) \mathbf{E}(\mathbf{x}, t) d\mathbf{x} + \int_{\mathbb{R}^3} \partial_t \mathbf{M}_{\mathrm{tot}}(\mathbf{x}, t) \mathbf{H}(\mathbf{x}, t) d\mathbf{x}. \quad (A.2) $$

Differentiating the constitutive laws (1.3) in time, we have

$$ \partial_t \mathbf{P}_{\mathrm{tot}}(\mathbf{x}, t) = \varepsilon_0 \int_0^t \chi_e'(t-s) \mathbf{E}(\mathbf{x}, s) \, ds + \varepsilon_0 \chi_e(0) \mathbf{E}(\mathbf{x}, t) $$

$$ \partial_t \mathbf{M}_{\mathrm{tot}}(\mathbf{x}, t) = \mu_0 \int_0^t \chi_m'(t) \mathbf{H}(\mathbf{x}, t-s) \, ds + \mu_0 \chi_m(0) \mathbf{H}(\mathbf{x}, t) $$

from which we deduce that

$$ \mathcal{I}(t) = \mathcal{I}_e(t) + \mathcal{I}_m(t) \quad \text{with} $$

$$ \mathcal{I}_e(t) := \varepsilon_0 \chi_e(0) \int_{\mathbb{R}^3} |\mathbf{E}(\mathbf{x}, t)|^2 d\mathbf{x} + \varepsilon_0 \int_{\mathbb{R}^3} \left( \mathbf{E}(\mathbf{x}, t) \cdot \int_0^t \chi_e'(t-s) \mathbf{E}(\mathbf{x}, s) \, ds \right) d\mathbf{x}, \quad (A.3) $$

$$ \mathcal{I}_m(t) := \mu_0 \chi_m(0) \int_{\mathbb{R}^3} |\mathbf{H}(\mathbf{x}, t)|^2 d\mathbf{x} + \mu_0 \int_{\mathbb{R}^3} \left( \mathbf{H}(\mathbf{x}, t) \cdot \int_0^t \chi_m'(t-s) \mathbf{H}(\mathbf{x}, s) \, ds \right) d\mathbf{x}. $$

It remains to transform $\mathcal{I}_e(t)$ and $\mathcal{I}_m(t)$.
The basic technical ingredient concerns convolution type quadratic forms:

$$ Qu(t) := \int_{0}^{t} k'(t-s) u(s)u'(t) ds, \qquad (A.4) $$

where $k(t)$ is a given convolution kernel. When $k(t) = C\delta(t)$, formally $Qu(t) = C \frac{d}{dt}|u(t)|^2$. The next lemma (lemma 3.2 in [16]) generalizes this observation to a smooth kernel. For completeness, we provide here a constructive proof (which is not given in [16]).

**Lemma A.1.** Given $k \in C^2(\mathbb{R}^+)$ and $u \in C^1(\mathbb{R}^+)$ with $u(0) = 0$, $Qu(t)$ given by (A.4) satisfies

$$ \left| \begin{aligned} Qu(t) = & \frac{1}{2} \frac{d}{dt} \left[ (k(t) - k(0)) |u(t)|^2 - \int_0^t k'(t-s) (u(s) - u(t))^2 ds \right] \\ & - \frac{1}{2} k'(t) |u(t)|^2 + \frac{1}{2} \int_0^t k''(t-s) (u(s) - u(t))^2 ds. \end{aligned} \right. \qquad (A.5) $$

*Proof.* The guiding idea is to make appear time derivatives of square quantities in the expression of $Qu(t)$.

The main trick is to write $u(s) = u(t) + (u(s) - u(t))$ in (A.4), in such a way that

$$ Qu(t) = \frac{1}{2} (k(t) - k(0)) \frac{d}{dt} |u(t)|^2 + \int_{0}^{t} k'(t-s) (u(s) - u(t)) u'(t) ds. \qquad (A.6) $$

On the one hand, one has

$$ \frac{1}{2} (k(t) - k(0)) \frac{d}{dt} |u(t)|^2 = \frac{1}{2} \frac{d}{dt} \left( (k(t) - k(0)) |u(t)|^2 \right) - \frac{1}{2} k'(t) |u(t)|^2. \qquad (A.7) $$

On the other hand, observing that $(u(t) - u(s)) u'(t) = \frac{1}{2} \frac{d}{dt} [(u(s) - u(t))^2]$ we have

$$ \int_{0}^{t} k'(t-s) (u(s) - u(t)) u'(t) ds = -\frac{1}{2} \int_{0}^{t} k'(t-s) \frac{d}{dt} [(u(s) - u(t))^2] ds, $$

i.e., since $k'(t-s) \frac{d}{dt} [(u(s) - u(t))^2] = \frac{d}{dt} [k'(t-s) (u(s) - u(t))^2] - k''(t-s) (u(s) - u(t))^2$,

$$ \left| \begin{aligned} \int_{0}^{t} k'(t-s) (u(s) - u(t)) u'(t) ds = & -\frac{1}{2} \int_{0}^{t} \frac{d}{dt} [k'(t-s) (u(s) - u(t))^2] ds \\ & + \frac{1}{2} \int_{0}^{t} k''(t-s) (u(s) - u(t))^2 ds. \end{aligned} \right. $$

Finally $\int_{0}^{t} \frac{d}{dt} [k'(t-s) (u(s) - u(t))^2] ds = \frac{d}{dt} \int_{0}^{t} k'(t-s) (u(s) - u(t))^2 ds$, thus

$$ \left| \begin{aligned} \int_{0}^{t} k'(t-s) (u(s) - u(t)) u'(t) ds = & -\frac{1}{2} \frac{d}{dt} \int_{0}^{t} k'(t-s) (u(s) - u(t))^2 ds \\ & + \frac{1}{2} \int_{0}^{t} k''(t-s) (u(s) - u(t))^2 ds. \end{aligned} \right. \qquad (A.8) $$

Finally, substituting (A.7) and (A.8) in (A.6) leads to (A.5). □

Now, we wish to transform the integrand (in space) in the second term of the expression (A.3) of $I_e(t)$ by making appear a quantity of the form (A.4). For this, it is useful to introduce primitives (in time) of
the fields. More precisely, for $\mathbf{F} = \mathbf{E}, \mathbf{H}, \mathbf{D}, \mathbf{B}, \mathbf{P}, \mathbf{M}$, we define for $t \ge 0$

$$ \mathbf{F}_p(\mathbf{x}, t) = \int_0^t f(\mathbf{x}, s) \, ds, \quad \forall t \ge 0 \text{ and a.e. } \mathbf{x} \in \mathbb{R}^3. $$

In order to transform $\mathcal{I}_e(t)$, we first perform an integration by parts in time to get, since $\mathbf{E}_p(\mathbf{x}, 0) = 0$,

$$ \int_0^t \chi'_e(t-s) \mathbf{E}(\mathbf{x}, s) \, ds \equiv \int_0^t \chi'_e(t-s) \partial_t \mathbf{E}_p(\mathbf{x}, s) \, ds = \chi'_e(0) \mathbf{E}_p(\mathbf{x}, t) + \int_0^t \chi''_e(t-s) \mathbf{E}_p(\mathbf{x}, s) \, ds. $$

As a consequence, since $\mathbf{E} = \partial_t \mathbf{E}_p$, we have

$$ \mathbf{E}(\mathbf{x}, t) \cdot \int_0^t \chi'_e(t-s) \mathbf{E}(\mathbf{x}, s) \, ds = \frac{\chi'_e(0)}{2} \frac{d}{dt} |\mathbf{E}_p(\mathbf{x}, t)|^2 + \int_0^t \chi''_e(t-s) \mathbf{E}_p(\mathbf{x}, s) \cdot \partial_t \mathbf{E}_p(\mathbf{x}, t) \, ds. $$

The second term in the right hand side of the above expression is a sum of terms of the form (A.4) with $k = \chi'_e$. Thus, integrating in space and then using the lemma A.1, and substituting the resulting equality in the expression (A.3) of $\mathcal{I}_e(t)$ gives (note that the terms involving $\chi'_e(0)$ cancel each other)

$$ \begin{aligned} \mathcal{I}_e(t) = & \frac{\varepsilon_0}{2} \frac{d}{dt} \left( \chi'_e(t) \int_{\mathbb{R}^3} |\mathbf{E}_p(\mathbf{x}, t)|^2 \, d\mathbf{x} \right) - \frac{\varepsilon_0}{2} \frac{d}{dt} \left( \int_0^t \chi''_e(t-s) \left( \int_{\mathbb{R}^3} |\mathbf{E}_p(\mathbf{x}, t) - \mathbf{E}_p(\mathbf{x}, s)|^2 \, d\mathbf{x} \right) ds \right) \\ & + \varepsilon_0 \chi_e(0) \int_{\mathbb{R}^3} |\mathbf{E}(\mathbf{x}, t)|^2 \, d\mathbf{x} - \frac{\varepsilon_0}{2} \chi''_e(t) \int_{\mathbb{R}^3} |\mathbf{E}_p(\mathbf{x}, t)|^2 \, d\mathbf{x} \\ & + \frac{\varepsilon_0}{2} \int_0^t \chi'''_e(t-s) \left( \int_{\mathbb{R}^3} |\mathbf{E}_p(\mathbf{x}, t) - \mathbf{E}_p(\mathbf{x}, s)|^2 \, d\mathbf{x} \right) ds. \end{aligned} \quad (A.9) $$

Analogously, we have

$$ \begin{aligned} \mathcal{I}_m(t) = & \frac{\mu_0}{2} \frac{d}{dt} \left( \chi'_m(t) \int_{\mathbb{R}^3} |\mathbf{H}_p(\mathbf{x}, t)|^2 \, d\mathbf{x} \right) - \frac{\mu_0}{2} \frac{d}{dt} \left( \int_0^t \chi''_m(t-s) \left( \int_{\mathbb{R}^3} |\mathbf{H}_p(\mathbf{x}, t) - \mathbf{H}_p(\mathbf{x}, s)|^2 \, d\mathbf{x} \right) ds \right) \\ & + \mu_0 \chi_m(0) \int_{\mathbb{R}^3} |\mathbf{H}(\mathbf{x}, t)|^2 \, d\mathbf{x} - \frac{\mu_0}{2} \chi''_m(t) \int_{\mathbb{R}^3} |\mathbf{H}_p(\mathbf{x}, t)|^2 \, d\mathbf{x} \\ & + \frac{\mu_0}{2} \int_0^t \chi'''_m(t-s) \left( \int_{\mathbb{R}^3} |\mathbf{H}_p(\mathbf{x}, t) - \mathbf{H}_p(\mathbf{x}, s)|^2 \, d\mathbf{x} \right) ds. \end{aligned} \quad (A.10) $$

Finally, (1.27) is simply obtained by gathering (A.2), (A.3), (A.9) and (A.10).

### A.3 Estimating $\mathcal{L}_k^{(2)}(0)$ in the Lorentz case

Below, we use the notation of Section 3 and our goal is to proof the estimate of Lemma 3.3, namely

$$ \mathcal{L}_k^{(2)}(0) \le |\mathbb{E}_0(\mathbf{k})|^2 + |\mathbb{H}_0(\mathbf{k})|^2. \quad (A.11) $$

First, by definition of $\mathcal{L}_k(0)$ and since $\mathbb{P}(\mathbf{k}, 0) = \mathbb{M}(\mathbf{k}, 0) = \partial_t \mathbb{P}(\mathbf{k}, 0) = \partial_t \mathbb{M}(\mathbf{k}, 0) = 0$,

$$ \mathcal{L}_k(0) = \frac{1}{2} (\varepsilon_0 |\mathbb{E}_0(\mathbf{k})|^2 + \mu_0 |\mathbb{H}_0(\mathbf{k})|^2). \quad (A.12) $$

Next we estimate $\langle \mathbf{k} \rangle^{-2} \mathcal{L}_k^1(0) = \langle \mathbf{k} \rangle^{-2} \mathcal{E}_k^1(0) + \langle \mathbf{k} \rangle^{-2} \mathcal{E}_{\Omega, \mathbf{k}}^1(0)$.
From equations (3.7a) and (3.7b) at $t=0$, since $\partial_t P(k, 0) = \partial_t M(k, 0) = 0$,

$$ \langle \mathbf{k} \rangle^{-2} \mathcal{E}_{\mathbf{k}}^{1}(0) \lesssim \langle \mathbf{k} \rangle^{-2} (|\mathbf{k} \times \mathbb{E}_{0}(\mathbf{k})|^{2} + |\mathbf{k} \times \mathbb{H}_{0}(\mathbf{k})|^{2}) \lesssim |\mathbb{E}_{0}(\mathbf{k})|^{2} + |\mathbb{H}_{0}(\mathbf{k})|^{2}. $$

Since $P(k, 0) = M(k, 0) = \partial_t P(k, 0) = \partial_t M(k, 0) = 0$, from (3.7c) and (3.7d) at $t=0$, we have

$$ \langle \mathbf{k} \rangle^{-2} \mathcal{E}_{\Omega, \mathbf{k}}^{1}(0) \lesssim \langle \mathbf{k} \rangle^{-2} (|\partial_{t}^{2} P(\mathbf{k}, 0)|^{2} + |\partial_{t}^{2} M(\mathbf{k}, 0)|^{2}) \lesssim |\mathbb{E}_{0}(\mathbf{k})|^{2} + |\mathbb{H}_{0}(\mathbf{k})|^{2} $$

and as a consequence

$$ \langle \mathbf{k} \rangle^{-2} \mathcal{L}_{\mathbf{k}}^{1}(0) \lesssim |\mathbb{E}_{0}(\mathbf{k})|^{2} + |\mathbb{H}_{0}(\mathbf{k})|^{2}. \qquad (A.13) $$

Finally we estimate $\langle \mathbf{k} \rangle^{-4} \mathcal{L}_{\mathbf{k}}^{2}(0) = \langle \mathbf{k} \rangle^{-4} \mathcal{E}_{\mathbf{k}}^{2}(0) + \langle \mathbf{k} \rangle^{-4} \mathcal{E}_{\Omega, \mathbf{k}}^{2}(0)$.

For bounding $\langle \mathbf{k} \rangle^{-4} \mathcal{E}_{\mathbf{k}}^{2}(0)$, we differentiate (3.7a) (resp. (3.7b)) and evaluate the resulting equations at $t=0$ to express $\partial_{t}^{2}\mathbb{E}(\mathbf{k}, 0)$ (resp. $\partial_{t}^{2}\mathbb{H}(\mathbf{k}, 0)$) in terms of $\mathbf{k} \times \partial_{t}\mathbb{H}(\mathbf{k}, 0)$ and $\partial_{t}^{2}\mathbb{P}(\mathbf{k}, 0)$ (resp. $\mathbf{k} \times \partial_{t}\mathbb{E}(\mathbf{k}, 0)$ and $\partial_{t}^{2}\mathbb{M}(\mathbf{k}, 0)$). From this, we deduce

$$ \langle \mathbf{k} \rangle^{-4} \mathcal{E}_{\mathbf{k}}^{2}(0) \lesssim \langle \mathbf{k} \rangle^{-4} (|\mathbf{k}|^{2} \mathcal{E}_{\mathbf{k}}^{1}(0) + |\partial_{t}^{2} P(\mathbf{k}, 0)|^{2} + |\partial_{t}^{2} M(\mathbf{k}, 0)|^{2}) \lesssim \langle \mathbf{k} \rangle^{-2} \mathcal{L}_{\mathbf{k}}^{1}(0) \lesssim |\mathbb{E}_{0}(\mathbf{k})|^{2} + |\mathbb{H}_{0}(\mathbf{k})|^{2} $$

where for the last inequality we have used (A.13). For the last term we first observe that

$$ \langle \mathbf{k} \rangle^{-4} \mathcal{E}_{\Omega, \mathbf{k}}^{2}(0) \lesssim \langle \mathbf{k} \rangle^{-4} (|\partial_{t}^{3} P(\mathbf{k}, 0)|^{2} + |\partial_{t}^{3} M(\mathbf{k}, 0)|^{2} + |\partial_{t}^{2} P(\mathbf{k}, 0)|^{2} + |\partial_{t}^{2} M(\mathbf{k}, 0)|^{2}) $$

and we obtain $\partial_{t}^{3}P(\mathbf{k}, 0)$ (resp. $\partial_{t}^{3}M(\mathbf{k}, 0)$) in function of $\partial_{t}\mathbb{E}(\mathbf{k}, 0)$ and $\partial_{t}^{2}\mathbb{P}(\mathbf{k}, 0)$ (resp. $\partial_{t}\mathbb{H}(\mathbf{k}, 0)$ and $\partial_{t}^{2}\mathbb{M}(\mathbf{k}, 0)$) from equations (3.7c) and (3.7d) after time differentiation. This leads to

$$ \langle \mathbf{k} \rangle^{-4} \mathcal{E}_{\Omega, \mathbf{k}}^{2}(0) \lesssim \langle \mathbf{k} \rangle^{-4} \langle \mathbf{k} \rangle^{2} \mathcal{L}_{\mathbf{k}}^{1}(0) \lesssim \langle \mathbf{k} \rangle^{-2} \mathcal{L}_{\mathbf{k}}^{1}(0) \lesssim |\mathbb{E}_{0}(\mathbf{k})|^{2} + |\mathbb{H}_{0}(\mathbf{k})|^{2} $$

where uses (A.13) for the last inequality. Adding the last two inequalities we obtain

$$ \langle \mathbf{k} \rangle^{-4} \mathcal{L}_{\mathbf{k}}^{2}(0) \lesssim |\mathbb{E}_{0}(\mathbf{k})|^{2} + |\mathbb{H}_{0}(\mathbf{k})|^{2}. \qquad (A.14) $$

Finally, (A.11) is deduced from (A.12), (A.13) and (A.14).

## A.4 Well-posdness and regularity of the solutions of the Cauchy problem in generalized Lorentz media

The (Cauchy) problem (3.1) can be rewritten as a generalized Schrödinger evolution problem:

$$ \frac{dU}{dt} + i\mathcal{A}U = 0 \quad \text{with} \quad U(0) = U_0, \qquad (A.15) $$

where the Hamiltonian $\mathcal{A}$ is an unbounded operator on the Hilbert-space:

$$ \mathcal{H} := L^2(\mathbb{R}^3)^3 \times L^2(\mathbb{R}^3)^3 \times L^2(\mathbb{R}^3)^{3N_e} \times L^2(\mathbb{R}^3)^{3N_e} \times L^2(\mathbb{R}^3)^{3N_m} \times L^2(\mathbb{R}^3)^{3N_m}, \qquad (A.16) $$
endowed by the following inner product : for any $\mathbf{U} = (\mathbf{E}, \mathbf{H}, \mathbf{P}, \dot{\mathbf{P}}, \mathbf{M}, \dot{\mathbf{M}})$ and $\mathbf{U}' = (\mathbf{E}', \mathbf{H}', \mathbf{P}', \mathbf{M}', \dot{\mathbf{P}}', \dot{\mathbf{M}}')$, where $(\mathbf{P}, \mathbf{M})$ is defined as in (3.3), (and the same for $(\mathbf{P}', \mathbf{M}')$, $(\dot{\mathbf{P}}, \dot{\mathbf{M}})$ and $(\dot{\mathbf{P}}', \dot{\mathbf{M}}')$)

$$ \left| \begin{aligned} (\mathbf{U}, \mathbf{U}')_{\mathcal{H}} = \frac{\varepsilon_0}{2} (\mathbf{E}, \mathbf{E}')_{L^2} &+ \frac{\mu_0}{2} (\mathbf{H}, \mathbf{H}')_{L^2} + \frac{\varepsilon_0}{2} \sum_{j=1}^{N_e} \omega_{e,j}^2 \Omega_{e,j}^2 (\mathbf{P}_j, \mathbf{P}'_j)_{L^2} + \frac{\varepsilon_0}{2} \sum_{j=1}^{N_e} \Omega_{e,j}^2 (\dot{\mathbf{P}}_j, \dot{\mathbf{P}}'_j)_{L^2} \\ &+ \frac{\mu_0}{2} \sum_{\ell=1}^{N_m} \Omega_{e,\ell}^2 (\mathbf{M}_\ell, \mathbf{M}'_\ell)_{L^2} + \frac{\mu_0}{2} \sum_{\ell=1}^{N_m} \Omega_{m,\ell}^2 (\dot{\mathbf{M}}_\ell, \dot{\mathbf{M}}'_\ell)_{L^2} , \end{aligned} \right. $$

More precisely, if we introduce $\mathcal{A} : D(\mathcal{A}) \subset \mathcal{H} \to \mathcal{H}$ defined by

$$ \forall \mathbf{U} = (\mathbf{E}, \mathbf{H}, \mathbf{P}, \dot{\mathbf{P}}, \mathbf{M}, \dot{\mathbf{M}}) \in D(\mathcal{A}), \quad \mathcal{A}\mathbf{U} := -i \begin{pmatrix} \varepsilon_0^{-1} \nabla \times \mathbf{H} + \sum \Omega_{e,j}^2 \mathbf{P}_j \\ -\mu_0^{-1} \nabla \times \mathbf{E} + \sum \Omega_{m,\ell}^2 \mathbf{M}_\ell \\ -\dot{\mathbf{P}} \\ \alpha_{e,j} \dot{\mathbf{P}}_j + \omega_{e,j}^2 \mathbf{P}_j - \mathbf{E} \\ -\dot{\mathbf{M}} \\ \alpha_{m,\ell} \dot{\mathbf{M}}_m + \omega_{m,\ell}^2 \mathbf{M}_\ell - \mathbf{H} \end{pmatrix} \quad (A.17) $$

the domain $D(\mathcal{A})$ (dense in $\mathcal{H}$) being given by

$$ D(\mathcal{A}) := H(\mathrm{rot}; \mathbb{R}^3) \times H(\mathrm{rot}; \mathbb{R}^3) \times L^2(\mathbb{R}^3)^{3N_e} \times L^2(\mathbb{R}^3)^{3N_e} \times L^2(\mathbb{R}^3)^{3N_m} \times L^2(\mathbb{R}^3)^{3N_m}, \quad (A.18) $$

we observe that one can rewrite (3.1) as (A.15) with the initial condition $\mathbf{U}_0 = (\mathbf{E}_0, \mathbf{H}_0, 0, 0, 0, 0) \in \mathcal{H}$.

The well-posedness of (A.15) is ensured by the following lemma:

**Lemma A.2.** If the assumption (1.37) holds, the operator $-i\mathcal{A}$ is a maximal dissipative operator.

*Proof.* To show that $-i\mathcal{A}$ is a maximal dissipative operator (see [12], theorem 8 p 340) is equivalent to show one one hand that $-i\mathcal{A}$ is dissipative, i.e. $\mathrm{Im}(\mathcal{A}\mathbf{U}, \mathbf{U}) \le 0$ for all $\mathbf{U} \in D(\mathcal{A})$ and that there exists $\omega \in \mathbb{C}^+$ such that $\mathcal{A} - \omega\mathbf{I}$ is surjective. Thus the proof will consists in two steps.

**Step 1:** $-i\mathcal{A}$ is dissipative. Performing $(\mathcal{A}\mathbf{U}, \mathbf{U})_{\mathcal{H}}$ thanks to (A.4) and (A.17), after integration by parts, i. e. s $(\nabla \times \mathbf{H}, \mathbf{E})_{L^2} = (\mathbf{H}, \nabla \times \mathbf{E})_{L^2} = (\nabla \times \mathbf{E}, \mathbf{H})_{L^2}$, one finds that

$$ \mathrm{Im}(\mathcal{A}\mathbf{U}, \mathbf{U}) = - \sum_{j=1}^{N_e} \alpha_{e,j} \Omega_{e,j}^2 \|\dot{\mathbf{P}}_j\|_{L^2}^2 - \sum_{\ell=1}^{N_m} \alpha_{m,\ell} \Omega_{m,\ell}^2 \|\dot{\mathbf{M}}_\ell\|_{L^2}^2 \le 0. $$

**Step 2:** for any $\omega \in \mathbb{C}^+$, $(\mathcal{A} - \omega\mathbf{I})D(\mathcal{A}) = \mathcal{H}$. We prove also the injectivity of the operator $\mathcal{A} - \omega\mathbf{I}$ by showing that for any $\mathbf{F} = (\mathbf{e}, \mathbf{h}, \mathbf{p}, \dot{\mathbf{p}}, \mathbf{m}, \dot{\mathbf{m}}) \in \mathcal{H}$ the system:

$$ (\mathcal{A} - \omega \mathbf{I}) \mathbf{U} = \mathbf{F} \quad (A.19) $$

admits a unique solution $\mathbf{U} = (\mathbf{E}, \mathbf{H}, \mathbf{P}, \dot{\mathbf{P}}, \mathbf{M}, \dot{\mathbf{M}}) \in D(\mathcal{A})$. To prove this, one first eliminates $\dot{\mathbf{P}}$ and $\dot{\mathbf{M}}$ in the system (A.19) by using that $\dot{\mathbf{P}} = -i\omega \mathbf{P} - i\mathbf{p}$ and $\dot{\mathbf{M}} = -i\omega \mathbf{M} - i\mathbf{m}$ and obtain the following
expression for **P** and **M** in term of **E** and **H**:

$$ \left\{ \begin{array}{ll} \mathbf{P}_j = -\frac{\mathbf{E}}{q_{e,j}(\omega)} + \mathbf{F}_{p,j}(\omega) & \text{where } \mathbf{F}_{p,j}(\omega) = \frac{(-i\alpha_{e,j} - \omega)\mathbf{P}_j - i\dot{\mathbf{P}}_j}{q_{e,j}(\omega)}, & \text{(A.20a)} \\ \mathbf{M}_\ell = -\frac{\mathbf{H}}{q_{m,l}(\omega)} + \mathbf{F}_{m,\ell}(\omega) & \text{where } \mathbf{F}_{m,l}(\omega) = \frac{(-i\alpha_{m,\ell} - \omega)\mathbf{m}_\ell - i\dot{\mathbf{m}}_\ell}{q_{m,l}(\omega)}, & \text{(A.20b)} \end{array} \right. $$

$$ \mathbf{M}_{\ell} = -\frac{\mathbf{H}}{q_{m,l}(\omega)} + \mathbf{F}_{m,\ell}(\omega) \quad \text{where} \quad \mathbf{F}_{m,l}(\omega) = \frac{(-i\alpha_{m,\ell} - \omega)\mathbf{m}_{\ell} - i\dot{\mathbf{m}}_{\ell}}{q_{m,\ell}(\omega)}, \quad \text{(A.20b)} $$

where and $q_{e,j}(\omega) = \omega^2 + i\alpha_{e,j}\omega - \omega_{e,j}^2 \neq 0$ and $q_{m,l}(\omega) = \omega^2 + i\alpha_{m,l}\omega - \omega_{m,l}^2 \neq 0$ for $\omega \in \mathbb{C}^+$. We point out that if $\omega_*$ is zero of $q_e$, $-\bar{\omega}_*$ is also zero of $q_e$, thus $\mathrm{Im}(\omega_*) = -\alpha_{e,j}/2 \le 0$ and therefore $\omega \in \mathbb{C}^+$ is not a zero of $q_{e,j}$. The same holds for the polynomials $q_{m,l}$. Thus, it follows the expressions of $\dot{\mathbf{P}}$ and $\dot{\mathbf{M}}$ in terms of **E** and **H**:

$$ \left\{ \begin{array}{ll} \dot{\mathbf{P}}_j = \frac{i\omega \mathbf{E}}{q_{e,j}(\omega)} + \dot{\mathbf{F}}_{p,j}(\omega) & \text{where } \dot{\mathbf{F}}_{p,j}(\omega) = \frac{i\omega_{e,j}^2 \mathbf{P}_j - \omega \dot{\mathbf{P}}_j}{q_{e,j}(\omega)}, & \text{(A.21a)} \\ \dot{\mathbf{M}}_\ell = \frac{i\omega \mathbf{H}}{q_{m,l}(\omega)} + \dot{\mathbf{F}}_{m,\ell}(\omega) & \text{where } \dot{\mathbf{F}}_{m,l}(\omega) = \frac{i\omega_{m,l}^2 \mathbf{m}_\ell - \omega \dot{\mathbf{m}}_\ell}{q_{m,l}(\omega)}. & \text{(A.21b)} \end{array} \right. $$

$$ \dot{\mathbf{M}}_{\ell} = \frac{i\omega \mathbf{H}}{q_{m,l}(\omega)} + \dot{\mathbf{F}}_{m,\ell}(\omega) \quad \text{where} \quad \dot{\mathbf{F}}_{m,l}(\omega) = \frac{i\omega_{m,l}^2 \mathbf{m}_{\ell} - \omega \dot{\mathbf{m}}_{\ell}}{q_{m,l}(\omega)}. \quad \text{(A.21b)} $$

Eliminating $\dot{\mathbf{P}}_j$ and $\dot{\mathbf{M}}_\ell$ in the two first equations of (A.19) yields to the following system in **E** and **H**:

$$ \left\{ \begin{array}{ll} -\nabla \times \mathbf{H} - i\omega \varepsilon(\omega) \mathbf{E} = \mathbf{F}_e(\omega) & \text{with } \mathbf{F}_e(\omega) = \varepsilon_0 [i \mathbf{e} - \sum_{j=1}^{N_e} \Omega_{e,j}^2 \dot{\mathbf{F}}_{p,j}(\omega)], & \text{(A.22a)} \\ \nabla \times \mathbf{E} - i\omega \mu(\omega) \mathbf{H} = \mathbf{F}_h(\omega) & \text{with } \mathbf{F}_h(\omega) = \mu_0 [i \mathbf{h} - \sum_{l=1}^{N_m} \Omega_{m,l}^2 \dot{\mathbf{F}}_{m,l}(\omega)]. & \text{(A.22b)} \end{array} \right. $$

$$ \nabla \times \mathbf{E} - i\omega\mu(\omega)\mathbf{H} = \mathbf{F}_h(\omega) \quad \text{with} \quad \mathbf{F}_h(\omega) = \mu_0 [i\mathbf{h} - \sum_{l=1}^{N_m} \Omega_{m,l}^2 \dot{\mathbf{F}}_{m,l}(\omega)]. \quad \text{(A.22b)} $$

where $\varepsilon(\omega)$ and $\mu(\omega)$ are defined by (1.14). Substituting

$$ \mathbf{H} = \frac{i}{\omega\mu(\omega)} (-\nabla \times \mathbf{E} + \mathbf{F}_h(\omega)) \quad \text{(A.23)} $$

in (A.22a) yields

$$ -\frac{1}{\omega\mu(\omega)}\nabla \times (\nabla \times \mathbf{E}) + \omega\varepsilon(\omega)\mathbf{E} = -\frac{1}{\omega\mu(\omega)}\nabla \times \mathbf{F}_h(\omega) + i\mathbf{F}_e(\omega). \quad \text{(A.24)} $$

One shows using standard arguments that $\mathbf{E} \in H(\mathrm{rot}; \mathbb{R}^3)$ is solution of (A.24) if and only if $\mathbf{E}$ is solution to the following variational problem in $H_{\mathrm{curl}}(\mathbb{R}^3)$:

$$ a(\mathbf{E}, \psi) = l(\psi), \quad \forall \phi \in H(\mathrm{rot}; \mathbb{R}^3), $$

where the sequilinear form $a$ and the antilinear form $l$ are defined for all $\phi, \psi \in H(\mathrm{rot}; \mathbb{R}^3)$:

$$ a(\phi, \psi) = \int_{\mathbb{R}^3} \left( -\frac{(\nabla \times \phi) \cdot (\overline{\nabla \times \psi})}{\omega\mu(\omega)} + \omega\varepsilon(\omega)\phi \cdot \overline{\psi} \right) d\mathbf{x}, \quad l(\psi) = \int_{\mathbb{R}^3} \left( -\frac{\mathbf{F}_h(\omega)}{\omega\mu(\omega)} \cdot (\overline{\nabla \times \psi}) + i\mathbf{F}_e(\omega) \cdot \overline{\psi} \right) d\mathbf{x}. $$

By using the Cauchy-Schwarz inequality, it is clear that $a$ (resp. $l$) is continuous on $H(\mathrm{rot}; \mathbb{R}^3)^2$ (resp. on $H(\mathrm{rot}; \mathbb{R}^3)$). Furthermore, $\omega \mapsto \omega\varepsilon(\omega)$ and $\omega \mapsto \omega\mu(\omega)$ are Herglotz functions. Thus, $\omega \mapsto -1/(\omega\mu(\omega))$ is aslo an Herglotz function. Hence, one shows that for all $\omega \in \mathbb{C}^+$:

$$ |a(\phi, \phi)| \ge |\mathrm{Im} a(\phi, \phi)| \ge \gamma(\omega) \|\phi\|_{H(\mathrm{rot}; \mathbb{R}^3)}^2 \quad \text{where } \alpha(\omega) = \min \left( \mathrm{Im} \left( - (\omega\mu(\omega))^{-1} \right), \mathrm{Im} \left( \omega\varepsilon(\omega) \right) \right) > 0. $$
Thus, by the Lax-Milgram theorem, (A.24) admits a unique solution in $H(\text{rot}; \mathbb{R}^3)$. Then, from (A.23), (A.20a), (A.21a), (A.20b), (A.21b), $\mathbf{H}, \mathbf{P}, \dot{\mathbf{P}}, \mathbf{M}, \dot{\mathbf{M}}$ are defined uniquely in term of $\mathbf{E}$ as elements of $L^2(\mathbb{R}^3)^3$. Moreover, with (A.22a), $\nabla \times \mathbf{H} \in L^2(\mathbb{R}^3)^3$, thus (A.19) admits a unique solution $\mathbf{U}$ in $D(\mathbf{A})$. $\square$

Under the weak dissipation assumption (1.37), $-i\mathbf{A}$ is a maximal dissipative operator. Hence, $\mathbf{A}$ is a closed operator (see e.g. [12], theorem 8 page 340) and its spectrum $\sigma(\mathbf{A})$ is contained in the closure of the lower half plane $\overline{\mathbb{C}^-}$. Furthermore, it follows from the Lumer-Phillips theorem (see e.g. [12], theorem 7 pages 336-337) that $-i\mathbf{A}$ is the generator of a contraction semi-group $\{S(t)\}_{t \ge 0}$ of class $C^0$. This implies the following standard results on the well-posdness of (A.15) and the regularity of its solutions (where the Hilbert space $D(\mathbf{A})$ is endowed with the graph norm defined by: $\|\mathbf{U}\|_{D(\mathbf{A})}^2 = \|\mathbf{U}\|_{\mathcal{H}}^2 + \|\mathbf{A}\mathbf{U}\|_{\mathcal{H}}^2, \forall \mathbf{U} \in D(\mathbf{A})$.

**Proposition A.3.** Let $m$ be an integer satisfying $m \ge 1$. If the initial condition $\mathbf{U}_0 \in D(\mathbf{A}^m)$ then the Cauchy problem (A.15) admits a unique strong solution $\mathbf{U} \in C^m(\mathbb{R}^+, \mathcal{H}) \cap C^{m-1}(\mathbb{R}^+, D(\mathbf{A}))$ given by $\mathbf{U}(t) = S(t)\mathbf{U}_0$ for $t \ge 0$. If $\mathbf{U}_0 \in \mathcal{H}$, (A.15) admits a unique mild solution $\mathbf{U} \in C^0(\mathbb{R}^+, \mathcal{H})$.

For the proof of this classical result, we refer to [12] (theorem 1 page 399) for the case $m = 1$ and to [3] (Theorem 7.5 page 191) for the case $m > 1$. For the definition of mild solutions when the initial data $\mathbf{U}_0$ is in $\mathcal{H}$, we refer to [12] pages 404-405.

**Remark A.4.** Under (1.37), $-i\mathbf{A}$ is dissipative. The identity $d\mathcal{L}(t)/dt = 2\operatorname{Im}(\mathbf{A}\mathbf{U}(t), \mathbf{U}(t))_{\mathcal{H}} \le 0$ with $\mathcal{L}(t) = (\mathbf{U}(t), \mathbf{U}(t))_{\mathcal{H}}$ (the abstract version of (2.5) or (3.6)) implies the decay in time of $\mathcal{L}(t)$.

We now define $\mathcal{H}_{\text{div},0}$ the closed subspace of $\mathcal{H}$ given by

$$ \mathcal{H}_{\text{div},0} := \{ \mathbf{U} = (\mathbf{E}, \mathbf{H}, \mathbf{P}, \dot{\mathbf{P}}, \mathbf{M}, \dot{\mathbf{M}}) \in \mathcal{H} \mid \nabla \cdot \mathbf{E} = \nabla \cdot \mathbf{H} = \nabla \cdot \mathbf{P} = \nabla \cdot \dot{\mathbf{P}} = \nabla \cdot \mathbf{M} = \nabla \cdot \dot{\mathbf{M}} = 0 \}. \quad (\text{A.25}) $$

Finally, the following proposition explains that if the initial has divergence free, the fields solutions to (A.15) remains divergence free for any $t \ge 0$.

**Proposition A.5.** Let $\mathbf{U}_0 \in D(\mathbf{A}) \cap \mathcal{H}_{\text{div},0}$ and $\mathbf{U} : \mathbb{R}^+ \to \mathcal{H}$ be the unique strong solution of (A.15) with $\mathbf{U}_0$ as initial condition then $\mathbf{U}(t) \in D(\mathbf{A}) \cap \mathcal{H}_{\text{div},0}$ for any $t \ge 0$. If $\mathbf{U}_0 \in \mathcal{H}_{\text{div},0}$ then the mild solution $\mathbf{U}$ of (A.15) satisfies also $\mathbf{U}(t) \in \mathcal{H}_{\text{div},0}$ for any $t \ge 0$.

*Proof.* We assume first that the initial condition is regular enough, namely $\mathbf{U}_0 \in D(\mathbf{A}^2)$. Hence by Proposition A.3, the strong solution $\mathbf{U}$ of (A.15) belongs to $C^2(\mathbb{R}^+, \mathcal{H}) \cap C^1(\mathbb{R}^+, D(\mathbf{A}))$ and thus the evolution equation (A.15) (or equivalently the system of equations (3.1)) holds also at $t = 0$ and one can differentiate it for $t \ge 0$. Taking the divergence in the distributional sense of the equation (3.1a) for $t \ge 0$ and taking the divergence of the equations $\partial_t((3.1c))$ for $t \ge 0$ and $j \in \{1, \dots, N_e\}$ leads to:

$$ \partial_t \nabla \cdot \mathbf{E} + \sum_{j=1}^{N_e} \Omega_{e,j}^2 \partial_t \nabla \cdot \mathbf{P}_j = 0 \quad \text{and} \quad \partial_t^3 \nabla \cdot \mathbf{P}_j + \alpha_{e,j} \partial_t^2 \nabla \cdot \mathbf{P}_j + \omega_{e,j}^2 \partial_t \nabla \cdot \mathbf{P}_j = \partial_t \nabla \cdot \mathbf{E}. \quad (\text{A.26}) $$

Thus, substituting $\partial_t \nabla \cdot \mathbf{E}$ in the second equation of (A.26) gives

$$ \partial_t^3 \nabla \cdot \mathbf{P}_j + \alpha_{e,j} \partial_t^2 \nabla \cdot \mathbf{P}_j + \omega_{e,j}^2 \partial_t \nabla \cdot \mathbf{P}_j + \sum_{j=1}^{N_e} \Omega_{e,j}^2 \partial_t \nabla \cdot \mathbf{P}_j = 0 \quad \text{for } j \in \{1, \dots, N_e\} \text{ and } t \ge 0, \quad (\text{A.27}) $$

with initial conditions: $\nabla \cdot \mathbf{P}_j(\cdot, 0) = \partial_t \nabla \cdot \mathbf{P}_j(\cdot, 0) = 0$ (since $\mathbf{U}_0 \in \mathcal{H}_{\text{div},0}$). Then taking the divergence of equation (3.1c) evaluated at $t = 0$ yields

$$ \partial_t^2 \nabla \cdot \mathbf{P}_j(\cdot, 0) = -\alpha_{e,j} \partial_t \nabla \cdot \mathbf{P}_j(\cdot, 0) - \omega_{e,j}^2 \partial_t \nabla \cdot \mathbf{P}_j(\cdot, 0) + \partial_t \nabla \cdot \mathbf{E}(\cdot, 0) = 0. $$
Thus, from (A.27), one deduces $\nabla \cdot \mathbf{P}_j(\cdot, t) = \nabla \cdot \partial_t \mathbf{P}_j(\cdot, t) = 0$ for any $t \ge 0$ and $j \in \{1, \dots, N_e\}$. It follows with (A.26) that $\partial_t \nabla \cdot \mathbf{E}(\cdot, t) = 0$ with $\nabla \cdot \mathbf{E}(\cdot, 0) = 0$. Therefore, one gets also $\nabla \cdot \mathbf{E}(\cdot, t) = 0$. Similarly using (3.1b) and (3.1d) and their derivatives, one shows that: $\nabla \cdot \mathbf{M}_\ell(\cdot, t) = \nabla \cdot \partial_t \mathbf{M}_\ell(\cdot, t) = \nabla \cdot \mathbf{H}(\cdot, t) = 0$ $t \ge 0$ and $\ell \in \{1, \dots, N_e\}$ and concludes that $\mathbf{U}(t) \in D(\mathbb{A}) \cap \mathcal{H}_{\text{div},0}$ for any $t \ge 0$.

If $\mathbf{U}_0$ is less regular, i.e. $\mathbf{U}_0$ in $\mathcal{H}_{\text{div},0}$, one obtains that $\mathbf{U}(t) = S(t)\mathbf{U}_0 \in \mathcal{H}_{\text{div},0}$ from the previous reasoning by using a density argument on the initial condition and the fact that the elements $S(t)$ of the $C^0$ semi-group corresponding to (A.15) are contractions. Moreover, if $\mathbf{U}_0 \in D(\mathbb{A})$, one has, by Proposition A.3, $\mathbf{U}(t) \in D(\mathbb{A})$ and thus $\mathbf{U}(t) \in D(\mathbb{A}) \cap \mathcal{H}_{\text{div},0}$. $\square$

**Remark A.6.** The (Cauchy) problem (3.1) assumes that all the resonance frequencies satisfy $\omega_{e,j}, \omega_{m,\ell} > 0$, in other words that $\varepsilon(\omega)$ and $\mu(\omega)$ only contained Lorentz terms (in the sense of Section 4.2). However, all the results in this appendix still hold with a similar proof if they contain Drude terms. In that case, if $N_{d,e}$ and $N_{d,m}$ are the number of electric and magnetic Drude terms, one has to redefine $\mathcal{H}$ in (A.16) as

$$ \mathcal{H} := L^2(\mathbb{R}^3)^3 \times L^2(\mathbb{R}^3)^3 \times L^2(\mathbb{R}^3)^{3(N_e-N_{d,e})} \times L^2(\mathbb{R}^3)^{3N_e} \times L^2(\mathbb{R}^3)^{3(N_m-N_{d,m})} \times L^2(\mathbb{R}^3)^{3N_m} $$

($\mathbf{P}_j$ (resp. $\mathbf{M}_\ell$) are no longer unknowns of the evolution problem for Drude terms).

## References

[1] A. Bernland, A. Luger and M. Gustafsson, *Sum rules and constraints on passive systems*, J. Phys. A: Math. Theor. 44 (14) (2011), 145205.

[2] C. Berg, Stieltjes-Pick-Bernstein-Schoenberg and their connection to complete monotonicity. Positive Definite Functions: From Schoenberg to Space-Time Challenges, pp. 15-45, 2008.

[3] H. Brezis, Functional Analysis, Sobolev Spaces and Partial Differential Equations, Springer, 2011.

[4] A. Borichev and Y. Tomilov, Optimal polynomial decay of functions and operator semigroups, Math. Ann. 347 (2) (2010), pp. 455-478.

[5] M. Cassier, P. Joly and M. Kachanovska, Mathematical models for dispersive electromagnetic waves: an overview, Comput. Math. with Appl. 74 (11) (2017), pp. 2639-2928.

[6] M. Cassier and G. W. Milton, Bounds on Herglotz functions and fundamental limits of broadband passive quasistatic cloaking, J. Math. Phys. 58 (7) (2017), 071504.

[7] M. Cessenat, Mathematical methods in electromagnetism: linear theory and applications. Vol. 41., World scientific, 1996.

[8] M. Conti, S. Gatti, V. Pata, Uniform decay properties of linear Volterra integro-differential equations, Math. Models Methods Appl. Sci. 18 (1) (2008), pp. 21-45

[9] C. Dafermos, Asymptotic stability in viscoelasticity, Arch. Ration. Mech. Anal. 37 (1970), pp. 297-308.

[10] V. Danese, P. G. Geredeli and V. Pata, Exponential attractors for abstract equations with memory and applications to viscoelasticity, Discrete Contin. Dyn. Syst. 35 (7) (2015), pp. 2881-2904.

[11] R. Dautray and J.-L. Lions. Mathematical analysis and numerical methods for science and technology: volume 3 Spectral Theory and Applications, Springer Science & Business Media, 1999.
[12] R. Dautray and J.-L. Lions. Mathematical analysis and numerical methods for science and technology: volume 5 Evolution problems I, Springer Science & Business Media, 1999.

[13] A. Figotin and J. H. Schenker, Spectral theory of time dispersive and dissipative systems, J. Stat. Phys. 118 (1) (2005), pp. 199–263.

[14] B. Gralak, A. Tip, Macroscopic Maxwell's equations and negative index materials, J. Math. Phys. 51 (5) (2010), 052902.

[15] A. Hanyga and M. Seredyńska, On a mathematical framework for the constitutive equations of anisotropic dielectric relaxation, J. Stat. Phys. 131 (2) (2008), pp. 269–303.

[16] J. E. Muñoz Rivera, M. Grazia Naso and E. Vuk, Asymptotic behaviour of the energy for electromagnetic systems with memory, Math. Meth. Appl. Sci. 27 (7) (2004), pp. 819–841.

[17] S. Nicaise, Stabilization and asymptotic behavior of dispersive medium models, Syst. Control Lett. 61 (5) (2012), pp. 638–648.

[18] S. Nicaise, C. Pignotti, Asymptotic behavior of dispersive electromagnetic waves in bounded domains, Z. Angew. Math. Phys. 71 (3) (2020), pp. 71–76.

[19] S. Nicaise, Stability and asymptotic properties of dissipative evolution equations coupled with ordinary differential equations, Math. Control Relat. F. (2022).

[20] M. Fabrizio, F. Morro, Mathematical problems of linear viscoelasticity, *SIAM Studies in Applied Mathematics*, 1987.

[21] M. Fabrizio, F. Morro, Thermodynamics of Electromagnetic Isothermal Systems with Memory, J. Non-Equilib Thermodyn. 22 (1997), pp. 110–128.

[22] R. Schilling, R. Song and Z. Vondracek. Bernstein functions, Theory and applications, *De Gruyter*, 2010.

[23] A. T. Welters, Y. Avniel and S. G. Johnson, Speed-of-light limitations in passive linear media, Phys. Rev. A 90 (2) (2014), 023847.

[24] A. H. Zemanian, Realizability theory for continuous linear systems, *Courier Corporation*, 1972.