# Semiclassical numerical modeling of gain materials with a high-order Implicit-Explicit Discontinuous Galerkin time-domain solver

Cédric Legrand, Stéphane Lanteri, Stéphane Descombes

► To cite this version:

Cédric Legrand, Stéphane Lanteri, Stéphane Descombes. Semiclassical numerical modeling of gain materials with a high-order Implicit-Explicit Discontinuous Galerkin time-domain solver. RR-9579, INRIA Sophia Antipolis. 2025. hal-04987997

HAL Id: hal-04987997

https://inria.hal.science/hal-04987997v1

Submitted on 13 Mar 2025

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
Ínria

# Semiclassical numerical modeling of gain materials with a high-order Implicit-Explicit Discontinuous Galerkin time-domain solver

Cédric Legrand, Stéphane Lanteri, Stéphane Descombes

## RESEARCH REPORT

Nº 9579

March 2025

Project-Teams Atlantis
<figure><img src="image_2.png" /></figure>
Inria

Semiclassical numerical modeling of gain
materials with a high-order Implicit-Explicit
Discontinuous Galerkin time-domain solver

Cédric Legrand, Stéphane Lanteri, Stéphane Descombes

Project-Teams Atlantis

Research Report n° 9579 — March 2025 — 34 pages

**Abstract:** This study focuses on the numerical modeling of wave propagation in gain media, a critical element in laser physics that amplifies electromagnetic fields through electron interactions. The gain process is represented by a four-level atomic model combining Maxwell's equations with nonlinear Ordinary Differential Equations (ODEs) to describe the evolution of electronic populations. While the Finite Difference Time Domain (FDTD) method has traditionally dominated this domain, this work introduces a high-order Discontinuous Galerkin Time-Domain (DGTD) method, tailored for solving such complex systems in 3D. The DGTD method is implemented with a second-order Leap-Frog temporal scheme and includes approximations for nonlinear terms. Stability is proven via energy estimates, both for the continuous problem and its discrete counterpart. Numerical validation, performed on a 3D cubic cavity with manufactured solutions, confirms the method's accuracy and convergence. The work highlights the potential of DGTD methods to enhance the precision and flexibility of numerical simulations in laser physics.

**Key-words:** Laser physics, gain medium, time-domain Maxwell equations, rate equations, discontinuous Galerkin method, implicit-explicit temporal scheme

RESEARCH CENTRE
Centre Inria d'Université Côte d'Azur

2004 route des Lucioles - BP 93
06902 Sophia Antipolis Cedex
Modélisation numérique semi-classique des matériaux
amplificateurs avec un solveur implicite-explicite Galerkin
discontinu en domaine temporel.

**Résumé** : Cette étude porte sur la modélisation numérique de la propagation des ondes dans les milieux amplificateurs, un élément clé de la physique des lasers qui amplifie les champs électromagnétiques via des interactions avec les électrons. Le processus d'amplification est représenté par un modèle atomique à quatre niveaux, combinant les équations de Maxwell avec des équations différentielles ordinaires (EDO) non linéaires pour décrire l'évolution des populations électroniques. Bien que la méthode des différences finies dans le domaine temporel (FDTD) ait longtemps dominé ce domaine, ce travail propose une méthode de Galerkin Discontinue en domaine temporel (DGTD) d'ordre élevé, spécialement adaptée à la résolution de tels systèmes complexes en 3D. La méthode DGTD est implémentée avec un schéma temporel Leap-Frog d'ordre deux et comprend des approximations pour les termes non linéaires. La stabilité est démontrée grâce à des estimations d'énergie pour le problème continu et son équivalent discret. La validation numérique, réalisée sur une cavité cubique 3D avec des solutions manufacturées, confirme la précision et la convergence de la méthode. Ce travail met en évidence le potentiel des méthodes DGTD pour améliorer la précision et la flexibilité des simulations numériques en physique des lasers.

**Mots-clés** : Physique des lasers, milieu amplificateur, équations de Maxwell en domaine temporel, équations de taux, méthode Galerkin discontinue, schéma temporel implicite
# Contents

<table><tr><td>1</td><td>Introduction</td><td>4</td></tr><tr><td>2</td><td>Physical context</td><td>5</td></tr><tr><td>2.1</td><td>Presentation of the physical problem</td><td>5</td></tr><tr><td>2.2</td><td>Presentation of the model</td><td>6</td></tr><tr><td>2.2.1</td><td>Problem setting</td><td>6</td></tr><tr><td>2.2.2</td><td>General model</td><td>6</td></tr><tr><td>2.2.3</td><td>Adimensionalization of the system</td><td>7</td></tr><tr><td>3</td><td>Mathematical model and properties study</td><td>9</td></tr><tr><td>3.1</td><td>Continuous system</td><td>9</td></tr><tr><td>3.1.1</td><td>Analytical framework</td><td>9</td></tr><tr><td>3.1.2</td><td>Bounds for the continuous energy</td><td>11</td></tr><tr><td>3.2</td><td>Discontinuous Galerkin method and time discretization</td><td>13</td></tr><tr><td>3.2.1</td><td>Semi-discrete problem in time</td><td>13</td></tr><tr><td>3.2.2</td><td>Discretization in space</td><td>14</td></tr><tr><td>3.3</td><td>Fully discrete energy</td><td>17</td></tr><tr><td>4</td><td>Numerical results</td><td>24</td></tr><tr><td>4.1</td><td>Manufactured solutions</td><td>25</td></tr><tr><td>4.2</td><td>Numerical results</td><td>26</td></tr><tr><td>5</td><td>Conclusion</td><td>31</td></tr><tr><td>A</td><td>Appendix</td><td>31</td></tr></table>
# 1 Introduction

This study is concerned with the numerical modeling of the propagation of waves in gain media. Such propagation media are used in laser physics and are characterized by an amplification of the electromagnetic field due to its interaction with a collection of electrons. Gain materials are mostly studied in laser physics, as seen in [12]. It has the property to raise light amplification thanks to stimulated photon emission. Laser systems typically comprise two reflector mirrors and an active gain medium, which together generate stimulated photons, electron transitions, and light amplification. To enable this phenomenon, electrons must be pumped to a discrete set of energy levels, setting the stage for lasing through successive excitation, relaxation, and emission processes. The theoretical modeling of gain media often employs the Maxwell-Bloch equations or similar frameworks to describe the interaction between the electromagnetic field and the media (see [15, 16]). Among these, the two-level and four-level atomic model stand out for their comprehensive representation of the excitation and relaxation processes within gain materials. Electrons are pumped to a high energy state, non-radiatively relax to an intermediate level, and subsequently transition radiatively or non-radiatively to a lower level, facilitating lasing. These dynamics are governed by a coupled system of Maxwell's equations and nonlinear rate equations, which capture the evolution of electronic populations across energy levels. In this article, we will describe a numerical resolution of the four-level energy model described in [1, 5, 17]. The study of gain media, particularly the four-level energy model, has been extensively explored, with much of the numerical work historically focused on methodological developments based on the Finite Difference Time Domain (FDTD) method (see [2, 4, 9, 10, 14]). This work introduces a high-order Discontinuous Galerkin Time-Domain (DGTD) method, recognized for its flexibility and high accuracy in solving hyperbolic problems over complex geometries. This method allows for an efficient and robust approach to tackle the coupled equations governing gain media.

DGTD solvers have been widely employed for Maxwell's equations over the past decades, with various time-integration schemes being proposed. Explicit schemes, such as Runge-Kutta 4, have been used for Maxwell's problems, as detailed in [3]. The mathematical model in this article couples Maxwell's equations with a set of nonlinear ordinary differential equations (ODEs) that describe the dynamics of the electronic populations. We propose a numerical study of the four-level energy model using a DGTD method combined with a second-order time-domain Leap-Frog scheme featuring centered fluxes. Notably, this approach is not fully explicit, as it requires solving an implicit system for the ODE component, which is why we refer to it as an IMEX-DG solver. This temporal scheme, originally introduced for DGTD solvers of Maxwell's equations in [8, 7], has shown strong potential in ensuring stability and accuracy of the resulting DGTD method. Further inspiration is drawn from applications of DGTD methods to the Drude-Lorenz model, as seen in [18], given the mathematical similarities with the four-energy model for gain. One can also take benefit from [13, 11], particularly for treating the nonlinear terms. In these works, with the Gauss-Legendre quadrature points, the authors define a Lagrange basis and an interpolation operator. This operator satisfies proprieties that allow to project the non-linear part in the Lagrange basis, without approximation. This method requires a structured squared grid and, in this work, we want to use a tetrahedral mesh, so the projection has to be done in a different way which will be presented here. To ensure the stability of a DG scheme, it is common to define and estimate an energy associated with the problem. Previous works have demonstrated this approach for classical Maxwell's equations, as seen in [8], and extended it to coupled systems involving ODEs (see [6, 3]). In this study, we extend these principles to our four-level model, establishing stability through energy estimates for the continuous problem and adapting
these estimates to our Leap-Frog scheme. Finally, we implemented a solver and validated it in a three-dimensional cubic cavity using manufactured solutions inspired by the methodology in [18].

This paper is organized as follows. First, in section 2, we delve into the physical context by introducing the four-level atomic system model. We detail the governing equations, including Maxwell's equations and the rate equations that characterize the dynamics of gain media. Subsequently, in section 3, we establish the mathematical framework, where we derive energy estimates for both the continuous and discrete systems. Particular attention is given to the stability analysis, highlighting the impact of boundary conditions such as the Perfect Electric Conductor (PEC) and Silver-Müller (SM) conditions. Here, we also focus on the numerical methods and we describe the implementation of the Discontinuous Galerkin Time-Domain (DGTD) method coupled with the second-order Leap-Frog temporal scheme. In addition, we address the challenges of handling nonlinear terms and provide insights into spatial discretization using tetrahedral meshes. Then, in section 4, we present the numerical results, where the proposed methodology is validated through simulations in a three-dimensional PEC cubic cavity, for which we derive a manufactured solution. These results highlight the convergence properties and robustness of our approach.

# 2 Physical context

## 2.1 Presentation of the physical problem

Gain materials have been extensively studied during the last decades. In laser physics, gain or amplification is a process due to the interaction of electromagnetic radiation with a collection of atoms. This interaction has two effects. On the one hand, the medium affects the field; on the other hand, the incident field causes a change in the material parameters. As a result of these effects, one observes the total field, which is the sum of the incident field and the field radiated by atoms. While the atoms lose energy through radiative or non-radiative mechanisms, the medium attenuates or amplifies and phase shifts the total field [5].

In this study, we focus on the four-level atomic system model, which provides a detailed representation of energy level dynamics within gain materials as seen in [9]. As represented in Figure 1, the model includes four distinct energy levels: ($E_0$), ($E_1$), ($E_2$), ($E_3$). Electrons are pumped from the ground state ($E_0$) to the third energy level ($E_3$) by an external electromagnetic mechanism. Subsequently, these electrons non-radiatively relax to the second level ($E_2$). Lasing occurs when the population density at ($E_2$) exceed the one at ($E_1$) and transition both radiatively and non-radiatively. Finally, electrons return to the ground state ($E_0$) via rapid, non radiative decay. The lasing and pumping processes are modeled as electric dipole transitions, characterized by the frequencies $\omega_a$ and $\omega_b$, respectively. These processes are governed by the decay rates $(1/\tau_{ij})_{1\le i,j\le 4}$, where $\tau_{ij}$ represents the transition time between energy levels i and j. Electrons are pumped to higher energy states by an external electromagnetic field, and under certain conditions, this process can be simplified to assume homogeneous pumping, particularly in thin gain slabs with low laser gain.

<figure><img src="image_3.png" /><figcaption>Figure 1: Four-level atomic system model.</figcaption></figure>
In the next sections, we present the four-level model with no pumping rate, the goal is to solve the model in the most general case with a DGTD. In the first part, we present the system that combines Maxwell’s equations and the rate equations.

## 2.2 Presentation of the model

### 2.2.1 Problem setting

The geometry of our problem is split into two domains. The first one is the host medium, surrounding the second one, the gain medium. We study the interaction between the gain medium and an incident electromagnetic wave. The theoretical study and the numerical study are conducted in 3D on the temporal domain $[0, T]$ with $T > 0$ fixed.

We introduce $\Omega_1, \Omega_2 \subset \mathbb{R}^3$, which are bounded domains corresponding to the gain medium and the host medium. $\Omega$ is equal to $\bar{\Omega}_1 \cup \Omega_2$ and $\Gamma = \partial\Omega$. This domain is illustrated in Figure 2.

<figure><img src="image_4.png" /><figcaption>Figure 2: The physical domain $\Omega$.</figcaption></figure>

### 2.2.2 General model

We denote by $\mathbf{E}$ and $\mathbf{H}$, respectively the electric field and the magnetic field vectors in $\mathbb{R}^3$. The electric permittivity and the magnetic permeability of the material are respectively denoted $\varepsilon_1$ and $\mu_1$ in $\Omega_1$ and $\varepsilon_2$ and $\mu_2$ in $\Omega_2$; they both are constants. We introduce $\mathbf{P}$ the electric polarization density of the gain material and $\mathbf{J}_p$ its time derivative in the gain material. We introduce $\mathbf{J}_p$ the extension by 0 of $\mathbf{J}_p$ from $\Omega_1$ to $\Omega$ and $\mathbf{E}, \mathbf{H}$ solutions of the following time-dependent Maxwell’s equation on the time domain $[0, T]$

$$ \forall i \in \{1, 2\} \quad \left\{ \begin{array}{ll} \varepsilon_i \frac{\partial \mathbf{E}}{\partial t} = \nabla \times \mathbf{H} - \bar{\mathbf{J}}_p & \text{in } \Omega_i \times [0, T], \\ \mu_i \frac{\partial \mathbf{H}}{\partial t} = -\nabla \times \mathbf{E} & \text{in } \Omega_i \times [0, T]. \end{array} \right. $$

These equations are completed with initial and boundary conditions, we consider two different types of boundary conditions for the electromagnetic fields on $\Gamma = \partial\Omega$ perfect electric conductor condition (PEC)

$$ \left\{ \begin{array}{l} \mathbf{n} \times \mathbf{H} = J_s, \\ \mathbf{n} \times \mathbf{E} = 0, \\ \mathbf{n} \cdot \mathbf{H} = 0, \\ \mathbf{n} \cdot \mathbf{E} = \frac{\rho_s}{\varepsilon}, \end{array} \right. $$

Inria
the surface current $J_s$ and the surface charge $\rho_s$ being unknowns, and Silver-Müller boundary condition (SM)

$$ \left\{ \begin{array}{l} \mathbf{n} \times (\mathbf{E} + \sqrt{\frac{\mu_2}{\epsilon_2}}(\mathbf{n} \times \mathbf{H})) = 0, \\ \mathbf{n} \times (\mathbf{H} - \sqrt{\frac{\epsilon_2}{\mu_2}}(\mathbf{n} \times \mathbf{E})) = 0. \end{array} \right. $$

We split the total electric polarization of the gain media $\mathbf{P}$ into two parts $\mathbf{P} = \mathbf{P}_a + \mathbf{P}_b$, $\mathbf{P}_a$ corresponding to the induced electric polarization density on the lasing transition between the upper ($E_2$) and lower ($E_1$) levels and $\mathbf{P}_b$ corresponding to the induced electric polarization density on the pumping transition between the ground state level ($E_0$) and the third level ($E_3$). For $1 \le i \le 4$, the electronic density at the $i^{th}$ energy level is denoted by $N_i$. The two polarization vectors $\mathbf{P}_a, \mathbf{P}_b$ and $N_i$, $1 \le i \le 4$, satisfy a system of ODEs in $\Omega_1$ that we now describe in detail. To do this, we reintroduce $\omega_b$ the pumping frequency, which is equal to $(E_3 - E_0)/\hbar$, where $\hbar$ is the Planck constant and $\omega_a$ the radiation frequency between the second level ($E_2$) and the first level ($E_1$) which is equal to $(E_2 - E_1)/\hbar$. As defined in [9] we introduce the constants $\gamma_a$ and $\gamma_b$ the line width of the atomic transition $\omega_a$ and $\omega_b$, and $\kappa_a$ and $\kappa_b$ the coupling strengths of $\mathbf{P}_a$ and $\mathbf{P}_b$ to the electric field. $\mathbf{P}_a$ and $\mathbf{P}_b$ satisfy the second-order system of differential equations

$$ \left\{ \begin{array}{l} \frac{\partial^2 \mathbf{P}_a}{\partial t^2} + \gamma_a \frac{\partial \mathbf{P}_a}{\partial t} + \omega_a^2 \mathbf{P}_a = -\kappa_a (N_2 - N_1) \mathbf{E} \quad \text{in } \Omega_1 \times [0, T], \\ \\ \frac{\partial^2 \mathbf{P}_b}{\partial t^2} + \gamma_b \frac{\partial \mathbf{P}_b}{\partial t} + \omega_b^2 \mathbf{P}_b = -\kappa_b (N_3 - N_0) \mathbf{E} \quad \text{in } \Omega_1 \times [0, T]. \end{array} \right. $$

The evolution of electronic densities is characterized by a system of ODEs called rate equations. The interaction between the electromagnetic field and the gain medium is characterized here by the presence of the field $\mathbf{E}$, which is coupled to the electric polarizations $\mathbf{P}_a$ and $\mathbf{P}_b$ induced by each transition in these equations. We also assume that at the initial time, all electrons are located at the lowest energy level with a density of $N_{tot}$ throughout the gain medium. If we use an external electromagnetic wave to pump the electrons from the ground state ($E_0$) to the third level ($E_3$), like in [9] we have to consider the following rate equations

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial N_3}{\partial t} = -\frac{N_3}{\tau_{32}} + \frac{1}{\hbar\omega_b} \mathbf{E} \cdot \frac{\partial \mathbf{P}_b}{\partial t} \quad \text{in } \Omega_1 \times [0, T], \\ \\ \displaystyle \frac{\partial N_2}{\partial t} = \frac{N_3}{\tau_{32}} - \frac{N_2}{\tau_{21}} + \frac{1}{\hbar\omega_a} \mathbf{E} \cdot \frac{\partial \mathbf{P}_a}{\partial t} \quad \text{in } \Omega_1 \times [0, T], \\ \\ \displaystyle \frac{\partial N_1}{\partial t} = \frac{N_2}{\tau_{21}} - \frac{N_1}{\tau_{10}} - \frac{1}{\hbar\omega_a} \mathbf{E} \cdot \frac{\partial \mathbf{P}_a}{\partial t} \quad \text{in } \Omega_1 \times [0, T], \\ \\ \displaystyle \frac{\partial N_0}{\partial t} = \frac{N_1}{\tau_{10}} - \frac{1}{\hbar\omega_b} \mathbf{E} \cdot \frac{\partial \mathbf{P}_b}{\partial t} \quad \text{in } \Omega_1 \times [0, T]. \end{array} \right. $$

The initial conditions for the electronic densities are taken as $N_i(\mathbf{x}, 0) = 0, 1 \le i \le 3$ and $N_0(\mathbf{x}, 0) = N_{tot}$.

### 2.2.3 Adimensionalization of the system

In this section we renormalize the system following an approach adapted from [18]. For this purpose, it is common to introduce $\epsilon_0$ the vacuum permittivity and $\mu_0$ the vacuum permeability,
and the relative permittivity and permeability $\epsilon_{r1}$, $\mu_{r1}$ in $\Omega_1$, and $\epsilon_{r2}$, $\mu_{r2}$ in $\Omega_2$; such as $\epsilon_i = \epsilon_0\epsilon_{ri}$ and $\mu_i = \mu_0\mu_{ri}$, for $i \in \{1, 2\}$. For greater clarity, we define

$$ \epsilon_r = \begin{cases} \epsilon_{r1} & \text{in } \Omega_1, \\ \epsilon_{r2} & \text{in } \Omega_2, \end{cases} \quad \text{and} \quad \mu_r = \begin{cases} \mu_{r1} & \text{in } \Omega_1, \\ \mu_{r2} & \text{in } \Omega_2. \end{cases} $$

In order to eliminate $\epsilon_0$ and $\mu_0$ from the system, we introduce $Z_0 = \sqrt{\mu_0/\epsilon_0}$ the vacuum impedance, $c_0 = 1/\sqrt{\epsilon_0\mu_0}$ the speed of light in the vacuum and the following new variables, $\tilde{H} = Z_0\mathbf{H}$, $\tilde{\mathbf{E}} = \mathbf{E}$, $\tilde{\mathbf{J}}_p = Z_0\mathbf{J}_p$, $\tilde{\mathbf{P}}_a = c_0Z_0\mathbf{P}_a$, $\tilde{\mathbf{P}}_b = c_0Z_0\mathbf{P}_b$, $\tilde{\mathbf{P}} = c_0Z_0\mathbf{P}$, $\tilde{N}_i = N_i/N_{\text{tot}}$, $1 \le i \le 4$ and $\tilde{t} = c_0t$. Then, the system becomes

$$ \left\{ \begin{aligned} \epsilon_r \frac{\partial \tilde{\mathbf{E}}}{\partial \tilde{t}} &= \nabla \times \tilde{\mathbf{H}} - \tilde{\mathbf{J}}_p & \text{in } \Omega \times [0, T], \\ \mu_r \frac{\partial \tilde{\mathbf{H}}}{\partial \tilde{t}} &= -\nabla \times \tilde{\mathbf{E}} & \text{in } \Omega \times [0, T], \\ \frac{\partial^2 \tilde{\mathbf{P}}_a}{\partial \tilde{t}^2} + \frac{\gamma_a}{c_0} \frac{\partial \tilde{\mathbf{P}}_a}{\partial \tilde{t}} + \frac{\omega_a^2}{c_0^2} \tilde{\mathbf{P}}_a &= -\frac{k_a Z_0 N_{\text{tot}}}{c_0} (\tilde{N}_2 - \tilde{N}_1) \tilde{\mathbf{E}} & \text{in } \Omega_1 \times [0, T], \\ \frac{\partial^2 \tilde{\mathbf{P}}_b}{\partial \tilde{t}^2} + \frac{\gamma_b}{c_0} \frac{\partial \tilde{\mathbf{P}}_b}{\partial \tilde{t}} + \frac{\omega_b^2}{c_0^2} \tilde{\mathbf{P}}_b &= -\frac{k_b Z_0 N_{\text{tot}}}{c_0} (\tilde{N}_3 - \tilde{N}_0) \tilde{\mathbf{E}} & \text{in } \Omega_1 \times [0, T], \\ \frac{\partial \tilde{N}_3}{\partial \tilde{t}} &= -\frac{\tilde{N}_3}{c_0 \tau_{32}} + \frac{1}{N_{\text{tot}} \hbar \omega_b c_0 Z_0} \tilde{\mathbf{E}} \cdot \frac{\partial \tilde{\mathbf{P}}_b}{\partial t} & \text{in } \Omega_1 \times [0, T], \\ \frac{\partial \tilde{N}_2}{\partial \tilde{t}} &= \frac{\tilde{N}_3}{c_0 \tau_{32}} - \frac{\tilde{N}_2}{c_0 \tau_{21}} + \frac{1}{N_{\text{tot}} \hbar \omega_a c_0 Z_0} \tilde{\mathbf{E}} \cdot \frac{\partial \tilde{\mathbf{P}}_a}{\partial t} & \text{in } \Omega_1 \times [0, T], \\ \frac{\partial \tilde{N}_1}{\partial \tilde{t}} &= \frac{\tilde{N}_2}{c_0 \tau_{21}} - \frac{\tilde{N}_1}{c_0 \tau_{10}} - \frac{1}{N_{\text{tot}} \hbar \omega_a Z_0 c_0} \tilde{\mathbf{E}} \cdot \frac{\partial \tilde{\mathbf{P}}_a}{\partial t} & \text{in } \Omega_1 \times [0, T], \\ \frac{\partial \tilde{N}_0}{\partial \tilde{t}} &= \frac{\tilde{N}_1}{c_0 \tau_{10}} - \frac{1}{N_{\text{tot}} \hbar \omega_b Z_0 c_0} \tilde{\mathbf{E}} \cdot \frac{\partial \tilde{\mathbf{P}}_b}{\partial t} & \text{in } \Omega_1 \times [0, T]. \end{aligned} \right. $$

For the rest of the paper, we introduce the following constants

$$ l_{i,j} = \tau_{i,j} c_0, \quad 1 \le i, j \le 4, \Gamma_{a,b} = \frac{\gamma_{a,b}}{c_0}, \Omega_{a,b} = \frac{\omega_{a,b}}{c_0}, $$

and

$$ r_{a,b} = \frac{1}{\hbar \omega_{a,b} c_0 Z_0 N_{\text{tot}}} \quad K_{a,b} = \frac{\kappa_{a,b} Z_0 N_{\text{tot}}}{c_0} = \kappa_{a,b} N_{\text{tot}} \mu_0. $$

To formulate a system of first-order ODEs, we introduce the temporal derivatives of the variables $\mathbf{P}_a$ and $\mathbf{P}_b$, respectively denoted as $\mathbf{J}_a$ and $\mathbf{J}_b$. By dropping the tilde notation and adding the initial condition, we obtain a new system of 10 equations and 10 variables

$$ \left\{ \begin{aligned} \epsilon_r \frac{\partial \mathbf{E}}{\partial t} &= \nabla \times \mathbf{H} - \overline{\mathbf{J}_a + \mathbf{J}_b} & \text{in } \Omega \times [0, T], \\ \mu_r \frac{\partial \mathbf{H}}{\partial t} &= -\nabla \times \mathbf{E} & \text{in } \Omega \times [0, T]. \end{aligned} \right. $$

Inria
$$ egin{cases} 
\frac{\partial \mathbf{J}_a}{\partial t} + \Gamma_a \mathbf{J}_a + \Omega_a^2 \mathbf{P}_a = -K_a (N_2 - N_1) \mathbf{E} & \text{in } \Omega_1 \times [0, T], \\
\frac{\partial \mathbf{J}_b}{\partial t} + \Gamma_b \mathbf{J}_b + \Omega_b^2 \mathbf{P}_b = -K_b (N_3 - N_0) \mathbf{E} & \text{in } \Omega_1 \times [0, T].
\end{cases} $$

$$ egin{cases} 
\frac{\partial N_3}{\partial t} = -\frac{N_3}{l_{32}} + r_b \mathbf{E} \cdot \frac{\partial \mathbf{P}_b}{\partial t} & \text{in } \Omega_1 \times [0, T], \\
\frac{\partial N_2}{\partial t} = \frac{N_3}{l_{32}} - \frac{N_2}{l_{21}} + r_a \mathbf{E} \cdot \frac{\partial \mathbf{P}_a}{\partial t} & \text{in } \Omega_1 \times [0, T], \\
\frac{\partial N_1}{\partial t} = \frac{N_2}{l_{21}} - \frac{N_1}{l_{10}} - r_a \mathbf{E} \cdot \frac{\partial \mathbf{P}_a}{\partial t} & \text{in } \Omega_1 \times [0, T], \\
\frac{\partial N_0}{\partial t} = \frac{N_1}{l_{10}} - r_b \mathbf{E} \cdot \frac{\partial \mathbf{P}_b}{\partial t} & \text{in } \Omega_1 \times [0, T].
\end{cases} $$

The initial conditions for the electronic densities are taken as $N_i(\mathbf{x}, 0) = 0, 1 \le i \le 4$ and $N_0(\mathbf{x}, 0) = 1$. The new parameters of the problem are summarized in table 1.

<table><caption>Table 1: Normalized parameters of the problem.</caption><thead><tr><th></th><th>E</th><th>H</th><th>P<sub>a,b</sub></th><th>J<sub>p,c</sub></th><th>N<sub>i</sub></th><th>t</th><th>τ<sub>i,j</sub></th><th>ω<sub>a,b</sub></th><th>γ<sub>a,b</sub></th><th>κ<sub>a,b</sub></th><th>r<sub>a,b</sub></th></tr></thead><tbody><tr><th>Original Unit</th><td>Vm<sup>-1</sup></td><td>Am<sup>-1</sup></td><td>sAm<sup>-2</sup></td><td>Am<sup>-2</sup></td><td>m<sup>-3</sup></td><td>s</td><td>s</td><td>s<sup>-1</sup></td><td>s<sup>-1</sup></td><td>Fm<sup>2</sup>s<sup>-2</sup></td><td>-</td></tr><tr><th>Normalized Unit</th><td>Vm<sup>-1</sup></td><td>Vm<sup>-1</sup></td><td>Vm<sup>-1</sup></td><td>Vm<sup>-2</sup></td><td>Ø</td><td>m</td><td>m</td><td>m<sup>-1</sup></td><td>m<sup>-1</sup></td><td>m<sup>-2</sup></td><td>V<sup>-2</sup>m<sup>2</sup></td></tr></tbody></table>

# 3 Mathematical model and properties study

This section outlines the mathematical framework used to derive an energy estimate for the continuous problem. To extend this energy estimate to the discrete problem, we introduce an IMEX time discretization paired with a Discontinuous Galerkin (DG) method and present the key results obtained for the discrete energy.

## 3.1 Continuous system

### 3.1.1 Analytical framework

In this section as in [8], we rewrite our problem as a Cauchy problem on a Hilbert space equipped with a weighted scalar product. These considerations are made to define and estimate the energy of this system. We introduce $L^2(\Omega)$ and $L^2(\Omega_1)$ the spaces of square-integrable vectorial functions on $\Omega$ and $\Omega_1$ and $L^2(\Omega_1) = L^2(\Omega_1)^3$, $L^2(\Omega) = L^2(\Omega)^3$. With these two different spaces, we can introduce

$$ \varepsilon = L^2(\Omega)^2 \times L^2(\Omega_1)^4 \times L^2(\Omega_1)^4. $$
We define a weighted scalar product on $\mathcal{E}$ denoted by $\langle \cdot, \cdot \rangle_{\mathcal{E}}$ and $\| \cdot \|_{\mathcal{E}}$ its associated norm, such that, for all $\mathcal{U} = (\mathbf{u}_i)_{1 \le i \le 10}$ and $\mathcal{V} = (\mathbf{v}_i)_{1 \le i \le 10}$ belonging to $\mathcal{E}^2$

$$ \langle \mathcal{U}, \mathcal{V} \rangle_{\mathcal{E}} = (\varepsilon_r \mathbf{u}_1, \mathbf{v}_1)_{\mathcal{L}^2(\Omega)} + (\mu_r \mathbf{u}_2, \mathbf{v}_2)_{\mathcal{L}^2(\Omega)} + \\ + \frac{1}{K_a} (\mathbf{u}_3, \mathbf{v}_3)_{\mathcal{L}^2(\Omega_1)} + \frac{1}{K_b} (\mathbf{u}_4, \mathbf{v}_4)_{\mathcal{L}^2(\Omega_1)} + \frac{\Omega_a^2}{K_a} (\mathbf{u}_5, \mathbf{v}_5)_{\mathcal{L}^2(\Omega_1)} + \frac{\Omega_b^2}{K_b} (\mathbf{u}_6, \mathbf{v}_6)_{\mathcal{L}^2(\Omega_1)} \\ + \frac{1}{r_b} (u_7, v_7)_{\mathcal{L}^2(\Omega_1)} + \frac{1}{r_a} (u_8, v_8)_{\mathcal{L}^2(\Omega_1)} + \frac{1}{r_a} (u_9, v_9)_{\mathcal{L}^2(\Omega_1)} + \frac{1}{r_b} (u_{10}, v_{10})_{\mathcal{L}^2(\Omega_1)}. $$

We now recall the following Sobolev spaces

$$ \begin{aligned} \mathbf{H}(\mathrm{curl}, \Omega) &= \{ \mathbf{v} \in \mathcal{L}^2(\Omega), \text{ such that } \nabla \times \mathbf{v} \in \mathcal{L}^2(\Omega) \}, \\ \mathbf{H}_0(\mathrm{curl}, \Omega) &= \{ \mathbf{v} \in \mathbf{H}(\mathrm{curl}, \Omega), \text{ such that } \mathbf{v} \times \mathbf{n}|_{\Gamma} = 0 \}, \\ \mathbf{H}(\mathrm{div}, \Omega) &= \{ \mathbf{v} \in \mathcal{L}^2(\Omega), \text{ such that } \mathrm{div} \, \mathbf{v} \in \mathcal{L}^2(\Omega) \}, \\ \mathbf{H}_0(\mathrm{div}, \Omega) &= \{ \mathbf{v} \in \mathbf{H}(\mathrm{div}, \Omega), \text{ such that } \mathbf{v} \cdot \mathbf{n}|_{\Gamma} = 0 \}. \end{aligned} $$

We introduce $\mathcal{H}$ the following subspace of $\mathcal{E}$

$$ \mathcal{H} = \mathbf{H}(\mathrm{div}, \Omega) \cap \mathbf{H}_0(\mathrm{curl}, \Omega) \times \mathbf{H}(\mathrm{curl}, \Omega) \cap \mathbf{H}_0(\mathrm{div}, \Omega) \times \mathcal{L}^2(\Omega_1)^4 \times \mathcal{L}^2(\Omega_1)^4. $$

We now define $\mathcal{A}$, $\mathcal{A}_{\mathrm{ME}}$ the operators containing the linear part of the ODEs and the Maxwell's equations and by $\mathcal{F}$, the nonlinear operator containing the nonlinear part of the equations. Let be $\mathcal{U}$ the vector of unknown, $\mathcal{U} = (\mathbf{E}, \mathbf{H}, \mathbf{J}_a, \mathbf{J}_b, \mathbf{P}_a, \mathbf{P}_b, N_3, N_2, N_1, N_0)^T$, $\mathcal{D}(\mathcal{A}) = \mathcal{D}(\mathcal{A}_{\mathrm{ME}}) = \mathcal{H}$ and for $\mathcal{U}$ in $\mathcal{H}$

$$ \begin{aligned} \mathcal{A}_{\mathrm{ME}} \mathcal{U} &= \left( \frac{1}{\varepsilon_r} \nabla \times \mathbf{H} - \frac{1}{\varepsilon_r} (\overline{\mathbf{J}_a + \mathbf{J}_b}), \quad \frac{1}{\mu_r} \nabla \times \mathbf{E}, \quad \mathbf{0}_{10,1} \right)^T, \\ \mathcal{A} \mathcal{U} &= \left( \mathbf{0}_{2,1}, -\Gamma_a \mathbf{J}_a - \Omega_a^2 \mathbf{P}_a, -\Gamma_b \mathbf{J}_b - \Omega_b^2 \mathbf{P}_b, \mathbf{J}_a, \mathbf{J}_b, -\frac{N_3}{l_{32}}, \frac{N_3}{l_{32}} - \frac{N_2}{l_{21}}, \frac{N_2}{l_{21}} - \frac{N_1}{l_{10}}, \frac{N_1}{l_{10}} \right)^T. \end{aligned} $$

To define the nonlinear operator $\mathcal{F}$ we assume that $\mathcal{D}(\mathcal{F})$ is a subspace of $\mathcal{H}$ such that the products between $N_i$ and $\mathbf{E}$, $1 \le i \le 4$ and the scalar products $\mathbf{E} \cdot \mathbf{J}_a$, $\mathbf{E} \cdot \mathbf{J}_b$ are well defined. For $\mathcal{U}$ in $\mathcal{D}(\mathcal{F})$, we define

$$ \mathcal{F}(\mathcal{U}) = (\mathbf{0}_{2,1}, -K_a(N_2-N_1)\mathbf{E}, -K_b(N_3-N_0)\mathbf{E}, \mathbf{0}_{2,1}, r_b\mathbf{E} \cdot \mathbf{J}_b, r_a\mathbf{E} \cdot \mathbf{J}_a, -r_a\mathbf{E} \cdot \mathbf{J}_a, -r_b\mathbf{E} \cdot \mathbf{J}_b)^T. $$

With these considerations, we can write the equations in a system form such as the problem with perfect electric boundary conditions is equivalent to

$$ \left\{ \begin{array}{l} \text{Find } \mathcal{U} \in \mathcal{D}(\mathcal{F}) \subset \mathcal{H} \text{ such that,} \\ \frac{\partial \mathcal{U}}{\partial t} = (\mathcal{A}_{\mathrm{ME}} + \mathcal{A}) \mathcal{U} + \mathcal{F}(\mathcal{U}), \\ \mathcal{U}(t=0) = \mathcal{U}_0, \end{array} \right. \qquad (1) $$

where $\mathcal{U}_0 \in \mathcal{D}(\mathcal{F})$ is an initial condition. On the other hand, by slightly changing $\mathcal{H}$ for $\mathcal{H}_1$ made out of functions $(\mathbf{u}, \mathbf{v}, w, z) \in \mathbf{H}^1(\Omega) \times \mathbf{H}^1(\Omega) \times \mathcal{L}^2(\Omega_1)^4 \times \mathcal{L}^2(\Omega_1)^4$ such that

$$ \mathbf{n} \times \mathbf{u} + \sqrt{\frac{\mu_r}{\epsilon_r}} \mathbf{n} \times \mathbf{n} \times \mathbf{v} = 0, $$

Inria
the problem with SM boundary conditions can also be written in the same form

$$ egin{cases} 
\text{Find } \mathcal{U} \in \mathcal{D}(\mathcal{F}) \subset \mathcal{H}_1 \text{ such that,} \\ 
\frac{\partial \mathcal{U}}{\partial t} = (\mathcal{A}_{\text{ME}} + \mathcal{A})\mathcal{U} + \mathcal{F}(\mathcal{U}), \\ 
\mathcal{U}(t=0) = \mathcal{U}_0, 
\end{cases} 
\quad (2) $$

where $\mathcal{U}_0 \in \mathcal{D}(\mathcal{F})$ is an initial condition. For the next section, we assume that the solutions of (1) and (2) exist and belong to $C^1([0, T], \mathcal{E})$ to carry out all necessary computations.

### 3.1.2 Bounds for the continuous energy

In this part, we prove some energy estimates with the two boundary conditions. Let $T > 0$ and $\mathcal{U}$ the solution of (1) or (2), we define $\xi$ the energy associated to the system such as, for $t \in [0, T]$, $\xi(t) = \|\mathcal{U}(t)\|_{\mathcal{E}}^2/2$. We have the following result:

**Theorem 3.1.1.** Let

$$ C = \sqrt{\frac{(K_a + K_b)}{\varepsilon_{r1}}} + \frac{1}{l_{32}} \sqrt{\frac{r_b}{r_a}} + \frac{1}{l_{21}} + \frac{1}{l_{10}} \sqrt{\frac{r_a}{r_b}}, $$

the solution $\mathcal{U}$ of (1) or (2) satisfies

$$ \forall t \in [0, T], \xi(t) \le \xi(0) \exp(2CT). \quad (3) $$

*Proof.* By differentiating the energy on $[0, T]$ we obtain, with (1) and (2) that for $t \in [0, T]$

$$ 
\begin{aligned}
\frac{\partial \xi}{\partial t}(t) &= \left\langle \frac{\partial \mathcal{U}}{\partial t}(t), \mathcal{U}(t) \right\rangle_{\mathcal{E}} \\
&= \langle (\mathcal{A}_{\text{ME}} + \mathcal{A})\mathcal{U}(t), \mathcal{U}(t) \rangle_{\mathcal{E}} + \langle \mathcal{F}(\mathcal{U})(t), \mathcal{U}(t) \rangle_{\mathcal{E}}.
\end{aligned}
 $$

We can notice that $\langle \mathcal{F}(\mathcal{U}), \mathcal{U} \rangle_{\mathcal{E}} = 0$ indeed

$$ 
\begin{aligned}
\langle \mathcal{F}(\mathcal{U}), \mathcal{U} \rangle_{\mathcal{E}} &= \frac{1}{K_a}(-K_a(N_2 - N_1)\mathbf{E}, \mathbf{J}_a)_{\mathbf{L}^2(\Omega_1)} + \frac{1}{K_b}(-K_b(N_3 - N_0)\mathbf{E}, \mathbf{J}_b)_{\mathbf{L}^2(\Omega_1)} \\
&\quad + \frac{1}{r_b}(r_b\mathbf{E} \cdot \mathbf{J}_b, N_3)_{\mathbf{L}^2(\Omega_1)} + \frac{1}{r_a}(r_a\mathbf{E} \cdot \mathbf{J}_a, N_2)_{\mathbf{L}^2(\Omega_1)} \\
&\quad + \frac{1}{r_a}(-r_a\mathbf{E} \cdot \mathbf{J}_a, N_1)_{\mathbf{L}^2(\Omega_1)} + \frac{1}{r_b}(-r_b\mathbf{E} \cdot \mathbf{J}_b, N_0)_{\mathbf{L}^2(\Omega_1)} \\
&= -((N_2 - N_1)\mathbf{E}, \mathbf{J}_a)_{\mathbf{L}^2(\Omega_1)} - ((N_3 - N_0)\mathbf{E}, \mathbf{J}_b)_{\mathbf{L}^2(\Omega_1)} \\
&\quad + (\mathbf{E} \cdot \mathbf{J}_b, (N_3 - N_0))_{\mathbf{L}^2(\Omega_1)} + (\mathbf{E} \cdot \mathbf{J}_a, (N_2 - N_1))_{\mathbf{L}^2(\Omega_1)} \\
&= 0.
\end{aligned}
 $$
For the linear part, we obtain the following expression

$$ \begin{aligned} \langle (\mathcal{A}_{\text{ME}} + \mathcal{A})\mathcal{U}, \mathcal{U} \rangle_{\mathcal{E}} = & \left( \varepsilon_r \frac{1}{\varepsilon_r} \nabla \times \mathbf{H}, \mathbf{E} \right)_{\mathbf{L}^2(\Omega)} + \varepsilon_{r1} \frac{1}{\varepsilon_{r1}} (-\mathbf{J}_a - \mathbf{J}_b, \mathbf{E})_{\mathbf{L}^2(\Omega_1)} \\ & + \left( \mu_r \left( -\frac{1}{\mu_r} \nabla \times \mathbf{E} \right), \mathbf{H} \right)_{\mathbf{L}^2(\Omega)} \\ & + \frac{1}{K_a} (-\Gamma_a \mathbf{J}_a - \Omega_a^2 \mathbf{P}_a, \mathbf{J}_a)_{\mathbf{L}^2(\Omega_1)} + \frac{1}{K_b} (-\Gamma_b \mathbf{J}_b - \Omega_b^2 \mathbf{P}_b, \mathbf{J}_b)_{\mathbf{L}^2(\Omega_1)} \\ & + \frac{\Omega_a^2}{K_a} (\mathbf{J}_a, \mathbf{P}_a)_{\mathbf{L}^2(\Omega_1)} + \frac{\Omega_b^2}{K_b} (\mathbf{J}_b, \mathbf{P}_b)_{\mathbf{L}^2(\Omega_1)} \\ & + \frac{1}{r_b} \left( -\frac{N_3}{l_{32}}, N_3 \right)_{\mathbf{L}^2(\Omega_1)} + \frac{1}{r_a} \left( \frac{N_3}{l_{32}} - \frac{N_2}{l_{21}}, N_2 \right)_{\mathbf{L}^2(\Omega_1)} \\ & + \frac{1}{r_a} \left( \frac{N_2}{l_{21}} - \frac{N_1}{l_{10}}, N_1 \right)_{\mathbf{L}^2(\Omega_1)} + \frac{1}{r_b} \left( \frac{N_1}{l_{10}}, N_0 \right)_{\mathbf{L}^2(\Omega_1)}. \end{aligned} $$

First, we have

$$ \begin{aligned} L(\mathbf{E}, \mathbf{H}) &= \left( \mu_r \left( -\frac{1}{\mu_r} \nabla \times \mathbf{E} \right), \mathbf{H} \right)_{\mathbf{L}^2(\Omega)} + \left( \varepsilon_r \frac{1}{\varepsilon_r} \nabla \times \mathbf{H}, \mathbf{E} \right)_{\mathbf{L}^2(\Omega)} \\ &= - \int_{\Omega} (\nabla \times \mathbf{E} \cdot \mathbf{H}) (x) dx + \int_{\Omega} (\nabla \times \mathbf{H} \cdot \mathbf{E}) (x) dx \\ &= - \left( \int_{\Omega} (\nabla \times \mathbf{H} \cdot \mathbf{E}) (x) dx - \int_{\Gamma} (\mathbf{E} \times \mathbf{n} \cdot \mathbf{H}) (x) dx \right) \\ &\quad + \int_{\Omega} (\nabla \times \mathbf{H} \cdot \mathbf{E}) (x) dx \\ &= \int_{\Gamma} (\mathbf{E} \times \mathbf{n} \cdot \mathbf{H}) (x) dx. \end{aligned} $$

The polarization vectors can be simplified in the equation so that we obtain

$$ \begin{aligned} \langle (\mathcal{A}_{\text{ME}} + \mathcal{A})\mathcal{U}, \mathcal{U} \rangle_{\mathcal{E}} = & - ((\mathbf{J}_a + \mathbf{J}_b), \mathbf{E})_{\mathbf{L}^2(\Omega_1)} + \int_{\Gamma} (\mathbf{E} \times \mathbf{n} \cdot \mathbf{H}) (x) dx \\ & - \frac{\Gamma_a}{K_a} (\mathbf{J}_a, \mathbf{J}_a)_{\mathbf{L}^2(\Omega_1)} - \frac{\Gamma_b}{K_b} (\mathbf{J}_b, \mathbf{J}_b)_{\mathbf{L}^2(\Omega_1)} \\ & - \frac{1}{r_b l_{32}} (N_3, N_3)_{\mathbf{L}^2(\Omega_1)} - \frac{1}{r_a l_{21}} (N_2, N_2)_{\mathbf{L}^2(\Omega_1)} - \frac{1}{r_a l_{10}} (N_1, N_1)_{\mathbf{L}^2(\Omega_1)} \\ & + \frac{1}{r_a l_{32}} (N_3, N_2)_{\mathbf{L}^2(\Omega_1)} + \frac{1}{r_a l_{21}} (N_2, N_1)_{\mathbf{L}^2(\Omega_1)} + \frac{1}{r_b l_{10}} (N_1, N_0)_{\mathbf{L}^2(\Omega_1)}. \end{aligned} $$

We have to estimate the integral on the boundary $\Gamma$, $\int_{\Gamma} (\mathbf{E} \times \mathbf{n} \cdot \mathbf{H})(x)dx$, for both boundaries conditions. We notice that for perfect electric conductor condition $\int_{\Gamma} (\mathbf{E} \times \mathbf{n} \cdot \mathbf{H})(x)dx = 0$ and for Silver Müller boundary condition we have

$$ \int_{\Gamma} (\mathbf{E} \times \mathbf{n} \cdot \mathbf{H}) (x) dx = - \sqrt{\frac{\mu_r}{\varepsilon_r}} \int_{\Gamma} \| (\mathbf{H} \times \mathbf{n}) (x) \|^2_{\mathbb{R}^3} dx \le 0, $$

then in both case, $\int_{\Gamma} (\mathbf{E} \times \mathbf{n} \cdot \mathbf{H})(x)dx \le 0$.
With these results, we obtain that for $t \in [0, T]$

$$ 
\begin{aligned}
\frac{\partial \xi}{\partial t}(t) = & -((\mathbf{J}_a + \mathbf{J}_b), \mathbf{E})_{\mathbf{L}^2(\Omega_1)} + \int_{\Gamma} (\mathbf{E} \times \mathbf{n} \cdot \mathbf{H}) (x) dx - \frac{\Gamma_a}{K_a} \|\mathbf{J}_a\|_{\mathbf{L}^2(\Omega_1)}^2 - \frac{\Gamma_b}{K_b} \|\mathbf{J}_b\|_{\mathbf{L}^2(\Omega_1)}^2 \\
& - \frac{1}{r_b l_{32}} \|N_3\|_{L^2(\Omega_1)}^2 - \frac{1}{r_a l_{21}} \|N_2\|_{L^2(\Omega_1)}^2 - \frac{1}{r_a l_{10}} \|N_1\|_{L^2(\Omega_1)}^2 \\
& + \frac{1}{r_a l_{32}} (N_3, N_2)_{L^2(\Omega_1)} + \frac{1}{r_a l_{21}} (N_2, N_1)_{L^2(\Omega_1)} + \frac{1}{r_b l_{10}} (N_1, N_0)_{L^2(\Omega_1)} \\
\le & -((\mathbf{J}_a + \mathbf{J}_b), \mathbf{E})_{\mathbf{L}^2(\Omega_1)} + \frac{1}{r_a l_{32}} (N_3, N_2)_{L^2(\Omega_1)} \\
& + \frac{1}{r_a l_{21}} (N_2, N_1)_{L^2(\Omega_1)} + \frac{1}{r_b l_{10}} (N_1, N_0)_{L^2(\Omega_1)} \\
\le & \|\mathbf{J}_a + \mathbf{J}_b\|_{\mathbf{L}^2(\Omega_1)} \|\mathbf{E}\|_{\mathbf{L}^2(\Omega_1)} + \frac{1}{r_a l_{32}} \|N_3\|_{L^2(\Omega_1)} \|N_2\|_{L^2(\Omega_1)} \\
& + \frac{1}{r_a l_{21}} \|N_2\|_{L^2(\Omega_1)} \|N_1\|_{L^2(\Omega_1)} + \frac{1}{r_b l_{10}} \|N_1\|_{L^2(\Omega_1)} \|N_0\|_{L^2(\Omega_1)} \\
\le & \frac{\sqrt{K_a + K_b}}{\sqrt{\varepsilon r_1}} \|\mathbf{U}(t)\|_{\mathcal{E}}^2 + \frac{\sqrt{r_b r_a}}{r_a l_{32}} \|\mathbf{U}(t)\|_{\mathcal{E}}^2 + \frac{\sqrt{r_a r_a}}{r_a l_{21}} \|\mathbf{U}(t)\|_{\mathcal{E}}^2 + \frac{\sqrt{r_a r_b}}{r_a l_{10}} \|\mathbf{U}(t)\|_{\mathcal{E}}^2 \\
\le & \left( \sqrt{\frac{(K_a + K_b)}{\varepsilon r_1}} + \frac{1}{l_{32}} \sqrt{\frac{r_b}{r_a}} + \frac{1}{l_{21}} + \frac{1}{l_{10}} \sqrt{\frac{r_a}{r_b}} \right) \|\mathbf{U}(t)\|_{\mathcal{E}}^2.
\end{aligned}
$$

Denoting by

$$ C = \sqrt{\frac{(K_a + K_b)}{\varepsilon r_1}} + \frac{1}{l_{32}} \sqrt{\frac{r_b}{r_a}} + \frac{1}{l_{21}} + \frac{1}{l_{10}} \sqrt{\frac{r_a}{r_b}}, $$

we obtain that for $t \in [0, T]$

$$ \frac{\partial \xi}{\partial t}(t) \le 2C\xi(t). $$

By integrating the previous inequality we obtain that for $t \in [0, T]$

$$ \xi(t) \le \xi(0) \exp(2Ct) \le \xi(0) \exp(2CT). $$

This concludes the proof of Theorem 3.1.1. □

## 3.2 Discontinuous Galerkin method and time discretization

### 3.2.1 Semi-discrete problem in time

For the time integration we adopt a second-order leap-frog scheme as in [18]. A common choice is to split every time iteration between $t$ and $t + \Delta t$ in two steps. $\mathbf{E}$ is approximated at even time stations $t_n = n\Delta t$, while $\mathbf{H}$, $\mathbf{J}_a$, $\mathbf{J}_b$, $\mathbf{P}_a$, $\mathbf{P}_b$ and the electronic densities $(N_i)_{i \in \{0,3\}}$ are approximated at odd time stations $t_{n+\frac{1}{2}} = (n + \Delta t/2)$. The temporal scheme consists in seeking the value of $\mathbf{E}^{n+1}$, $\mathbf{H}^{n+3/2}$, $\mathbf{J}_a^{n+3/2}$, $\mathbf{J}_b^{n+3/2}$, $\mathbf{P}_a^{n+3/2}$, $\mathbf{P}_b^{n+3/2}$, and $N_i^{n+3/2}$, knowing $\mathbf{E}^n$, $\mathbf{H}^{n+1/2}$, $\mathbf{J}_a^{n+1/2}$, $\mathbf{J}_b^{n+1/2}$, $\mathbf{P}_a^{n+1/2}$, $\mathbf{P}_b^{n+1/2}$, and $N_i^{n+1/2}$. For the approximation of the ODEs, we drew inspiration from [13] and [11] to ensure the stability of the method. For a vector variable $\mathbf{U}$ or a scalar variable $\mathbf{V}$, we define $\mathbf{U}^{[n+\frac{1}{2}]}$ and $\mathbf{V}^{[n+\frac{1}{2}]}$ such as

$$ \mathbf{U}^{[n+\frac{1}{2}]} = \frac{\mathbf{U}^{n+\frac{1}{2}} + \mathbf{U}^{n+3/2}}{2} \quad \text{and} \quad \mathbf{V}^{[n+\frac{1}{2}]} = \frac{\mathbf{V}^{n+\frac{1}{2}} + \mathbf{V}^{n+3/2}}{2}. $$
The second-order leap-frog scheme adapted to our problem is presented as follows

$$
\left\{
\begin{aligned}
& \epsilon_r \frac{\mathbf{E}^{n+1} - \mathbf{E}^n}{\Delta t} = \nabla \times \mathbf{H}^{n+\frac{1}{2}} - \mathbf{J}_a^{n+\frac{1}{2}} - \mathbf{J}_b^{n+\frac{1}{2}}, \\
& \mu_r \frac{\mathbf{H}^{n+3/2} - \mathbf{H}^{n+\frac{1}{2}}}{\Delta t} = -\nabla \times \mathbf{E}^{n+1}, \\
& \left\{
\begin{aligned}
& \frac{\mathbf{J}_a^{n+3/2} - \mathbf{J}_a^{n+\frac{1}{2}}}{\Delta t} + \Gamma_a \mathbf{J}_a^{[n+\frac{1}{2}]} + \Omega_a^2 \mathbf{P}_a^{[n+\frac{1}{2}]} = -K_a (N_2^{[n+\frac{1}{2}]} - N_1^{[n+\frac{1}{2}]}) \mathbf{E}^{n+1}, \\
& \frac{\mathbf{J}_b^{n+3/2} - \mathbf{J}_b^{n+\frac{1}{2}}}{\Delta t} + \Gamma_b \mathbf{J}_b^{[n+\frac{1}{2}]} + \Omega_b^2 \mathbf{P}_b^{[n+\frac{1}{2}]} = -K_b (N_3^{[n+\frac{1}{2}]} - N_0^{[n+\frac{1}{2}]}) \mathbf{E}^{n+1}, \\
& \frac{\mathbf{P}_a^{n+3/2} - \mathbf{P}_a^{n+\frac{1}{2}}}{\Delta t} = \mathbf{J}_a^{[n+\frac{1}{2}]}, \\
& \frac{\mathbf{P}_b^{n+3/2} - \mathbf{P}_b^{n+\frac{1}{2}}}{\Delta t} = \mathbf{J}_b^{[n+\frac{1}{2}]},
\end{aligned}
\right.
\end{aligned}
\right.
\left\{
\begin{aligned}
& \frac{N_3^{n+3/2} - N_3^{n+\frac{1}{2}}}{\Delta t} = -\frac{N_3^{[n+\frac{1}{2}]}}{l_{32}} + r_b \mathbf{E}^{n+1} \cdot \mathbf{J}_b^{[n+\frac{1}{2}]}, \\
& \frac{N_2^{n+3/2} - N_2^{n+\frac{1}{2}}}{\Delta t} = \frac{N_3^{[n+\frac{1}{2}]}}{l_{32}} - \frac{N_2^{[n+\frac{1}{2}]}}{l_{21}} + r_a \mathbf{E}^{n+1} \cdot \mathbf{J}_a^{[n+\frac{1}{2}]}, \\
& \frac{N_1^{n+3/2} - N_1^{n+\frac{1}{2}}}{\Delta t} = \frac{N_2^{[n+\frac{1}{2}]}}{l_{21}} - \frac{N_1^{[n+\frac{1}{2}]}}{l_{10}} - r_a \mathbf{E}^{n+1} \cdot \mathbf{J}_a^{[n+\frac{1}{2}]}, \\
& \frac{N_0^{n+3/2} - N_0^{n+\frac{1}{2}}}{\Delta t} = \frac{N_1^{[n+\frac{1}{2}]}}{l_{10}} - r_b \mathbf{E}^{n+1} \cdot \mathbf{J}_b^{[n+\frac{1}{2}]}. \\
& \frac{N_3^{n+3/2} - N_3^{n+\frac{1}{2}}}{\Delta t} = -\frac{N_3^{[n+\frac{1}{2}]}}{l_{32}} + r_b \mathbf{E}^{n+1} \cdot \mathbf{J}_b^{[n+\frac{1}{2}]},
\end{aligned}
\right.
$$

As we can observe, the temporal scheme for Maxwell's equations is fully explicit, as the current vectors are injected at the time iteration $n + \frac{1}{2}$. In contrast, the ODE system requires solving an implicit system where $J_{a,b}^{[n+3/2]}$, $P_{a,b}^{[n+3/2]}$, and $N_i^{[n+3/2]}$ are the unknowns. This involves inverting a matrix of size $16 \times 16$ and we then have a Implicit-Explicit temporal scheme.

### 3.2.2 Discretization in space

In this section we are going to set up the space discretization done with the discontinuous Galerkin method. We will not provide a detailed explanation of the implementation of this method, but we will emphasize the approximations used for handling nonlinear terms. Here we will directly introduce the matrix formulation for the fully discrete scheme.

Let $\Omega_{1h}$ and $\Omega_{2h}$ discretizations of, respectively, $\Omega_1$ and $\Omega_2$, relying on a quasi-uniform tri-angulation $\mathcal{T}_h = (T_i)_{[1,N]}$, with $N \in \mathbb{N}^*$. Let $\Omega_h = \Omega_{1h} \cup \Omega_{2h}$, a discretization of $\Omega$. We denote $a_{ij} = T_i \cap T_j$, an internal face of the discretization for $(i, j) \in [1, N]$ such as $T_i$ and $T_j$ are adjacent. Let $\mathbf{n}_{ij}$ be the unit normal vector of $a_{ik}$, oriented from $T_i$ to $T_j$. For $i \in [1, N]$, we denote $\mathcal{V}_i = \{j \in [1, N], T_i \cap T_j \neq \emptyset\}$. We introduce the discrete space $\mathcal{V}_h^k = \{v \in L^2(\Omega_1) \text{ such as for all } i, v|_{T_i} \in P^p(T_i)\}$, $P^p(T_i)$ being the set of polynomials on $T_i$ with degree $p$, for $p \in \mathbb{N}^*$. For each cell $T_i$ we define a set of scalar basis functions $(\phi_{ij})_{1 \le j \le d_i}$
where $d_i$ is the number of degrees of freedom per dimension.

For a vectorial semi-discrete variable by $\mathbf{U}^n$ and a scalar semi-discrete variable $V^n$, we respectively denote by $\mathbf{U}_h^n$ and $V_h^n$ the discrete in space and times variables. Their restrictions on each cell $T_i$ are denoted respectively $\mathbf{U}_i^n$ and $V_i^n$. For a vectorial vector $\mathbf{U}_i^n$, we introduce its coordinates such as

$$ \mathbf{U}_i^n = \begin{pmatrix} U_i^{xn} \\ U_i^{yn} \\ U_i^{zn} \end{pmatrix}. $$

Each of these variables is locally expended on the chosen set of basis functions

$$ V_i^n = \sum_{j=1}^{d_i} V_{ij}^n \phi_{ij}, \text{ and } \forall v \in \{x, y, z\}, U_i^{vn} = \sum_{j=1}^{d_i} U_{ij}^{vn} \phi_{ij}. $$

Therefore, in practice, we define vectors of $d_i$ components and $3d_i$ components

$$ \bar{V}_i^n = \begin{pmatrix} V_{i1}^n \\ \vdots \\ V_{id_i}^n \end{pmatrix} \quad \text{and} \quad \bar{\mathbf{U}}_i^n = \begin{pmatrix} (U_{ij}^{xn})_{1 \le j \le d_i} \\ (U_{ij}^{yn})_{1 \le j \le d_i} \\ (U_{ij}^{zn})_{1 \le j \le d_i} \end{pmatrix}. $$

The approximate fields $(E_h, H_h)$ inside the finite element $T_i$. We denote by $P_i = \text{Span}\{\phi_h\}_i$ the approximation space for the fields. The approximate fields, defined for all $i$ by $(E_h|_{T_i} = E_i, H_h|_{T_i} = H_i)$, are allowed to be discontinuous across element boundaries. For such a discontinuous field, we use centered flux and denoted them by $\{\mathbf{E}^{n+1}\}_{il}$ and $\{\mathbf{H}^{n+\frac{1}{2}}\}_{il}$ through any internal face $a_{il}$ as

$$ \{\mathbf{H}^{n+\frac{1}{2}}\}_{il} = \frac{\mathbf{H}_i^{n+\frac{1}{2}} + \mathbf{H}_l^{n+\frac{1}{2}}}{2}, \text{ and } \{\mathbf{E}^{n+1}\}_{il} = \frac{\mathbf{E}_i^{n+1} + \mathbf{E}_l^{n+1}}{2}. $$

We introduce the three matrices defined in [18], the mass matrix $\mathbf{M}_i$, the stiffness matrices, for $v \in \{x, y, z\}$, $\mathbb{K}_i^v$, the flux matrix $\mathbb{S}_i$ such as

$$ \mathbf{M}_i = \left( \int_{T_i} \phi_{ij} \phi_{ik} \right)_{(k,j) \in [1, d_i]^2}, \quad \mathbb{K}_i^v = \left( \int_{T_i} \phi_{ij} \frac{\partial \phi_{ik}}{\partial v} \right)_{(k,j) \in [1, d_i]^2} \quad \text{and} \quad \mathbb{S}_i = \left( \int_{a_{il}} \phi_{ij} \phi_{ik} \right)_{(k,j) \in [1, d_i]^2}. $$

To ensure clarity, we omit the variable of integration in the integrals for the remainder of this section.

With these three matrices we can define three other matrices that will be used to obtain a generalized formula for Maxwell's equation in the fully discrete scheme

$$ \bar{\mathbf{M}}_i = \begin{pmatrix} \mathbf{M}_i & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{M}_i & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{M}_i \end{pmatrix}, \quad \bar{\mathbf{K}} = \begin{pmatrix} \mathbb{K}_i^x \\ \mathbb{K}_i^y \\ \mathbb{K}_i^z \end{pmatrix} \quad \text{and} \quad \bar{\mathbf{S}}_{il} = \begin{pmatrix} \mathbb{S}_{il} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbb{S}_{il} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbb{S}_{il} \end{pmatrix}. $$

With the previous notations, we obtain directly the following matrix formulation for the Maxwell's equations

$$ \left\{ \begin{array}{l} \varepsilon_r \bar{\mathbf{M}}_i \frac{\bar{\mathbf{E}}_i^{n+1} - \bar{\mathbf{E}}_i^n}{\Delta t} = \bar{\mathbb{K}}_i \times \bar{\mathbf{H}}_i^{n+\frac{1}{2}} - \bar{\mathbf{M}}_i (\bar{\mathbf{J}}_{a_i}^{[n+\frac{1}{2}]} + \bar{\mathbf{J}}_{b_i}^{[n+\frac{1}{2}]} ) - \sum_{l \in \mathcal{V}_i} \bar{\mathbf{S}}_{il} (\{\mathbf{H}^{n+\frac{1}{2}}\}_{il} \times \mathbf{n}_{il}), \\ \\ \mu_r \bar{\mathbf{M}}_i \frac{\bar{\mathbf{H}}_i^{n+3/2} - \bar{\mathbf{H}}_i^{n+\frac{1}{2}}}{\Delta t} = -\bar{\mathbb{K}}_i \times \bar{\mathbf{E}}_i^{n+\frac{1}{2}} + \sum_{l \in \mathcal{V}_i} \bar{\mathbf{S}}_{il} (\{\mathbf{E}^{n+1}\}_{il} \times \mathbf{n}_{il}). \end{array} \right. \qquad (\text{ME}) $$
The presence of nonlinear terms complicates the task. Indeed, if we project the variables onto the polynomial basis as usual, we end up having to compute integrals involving a triple product of functions. For example, for two scalar variables $V_i$ and $W_i$ we would have in the weak formulation applied on a basis function $\phi_{ik}$, $k \in [1, d_i]$,

$$ \int_{T_i} V_i W_i \phi_{ik} = \sum_{j=1}^{d_i} \sum_{l=1}^{d_i} \left( V_{ij} W_{il} \int_{T_i} \phi_{ij} \phi_{il} \phi_{ik} \right). $$

In the matrix formulation, this would result in the presence of a tensor of size $d_i^3$. However, this is too complex and computationally demanding, we project directly each product of variables onto the function basis. With the notations from the previous example, we would thus obtain

$$ \int_{T_i} V_i W_i \phi_{ik} \approx \int_{T_i} \left( \sum_{r=1}^{d_i} \left[ \sum_{j=1}^{d_i} \sum_{l=1}^{d_i} V_{ij} W_{il} \phi_{ij} \phi_{il} \right] \middle|_{x_r} \phi_{ir} \right) \phi_{ik} \approx \sum_{r=1}^{d_i} \left( V_{ir} W_{ir} \int_{T_i} \phi_{ir} \phi_{ik} \right). $$

We can see that calculating the integral leads to the inclusion of the mass matrix in the matrix formulation, which ultimately simplifies when combined with the matrix in front of the other terms. We can directly apply this method of approximation to the non-linear terms of the rate equations and polarization equations. We obtain the following system of discrete equations

$$ \left\{ \begin{aligned} \frac{\overline{\mathbf{P}}_{a_i}^{n+3/2} - \overline{\mathbf{P}}_{a_i}^{n+1/2}}{\Delta t} &= \overline{\mathbf{J}}_{a_i}^{[n+1/2]}, \\ \frac{\overline{\mathbf{P}}_{b_i}^{n+3/2} - \overline{\mathbf{P}}_{b_i}^{n+1/2}}{\Delta t} &= \overline{\mathbf{J}}_{b_i}^{[n+1/2]}, \\ \frac{\overline{\mathbf{J}}_{a_i}^{n+3/2} - \overline{\mathbf{J}}_{a_i}^{n+1/2}}{\Delta t} &= -\Gamma_a \overline{\mathbf{J}}_{a_i}^{[n+1/2]} - \Omega_a^2 \overline{\mathbf{P}}_{a_i}^{[n+1/2]} - K_a \overline{\mathbf{E}} \overline{\mathbf{N}}_{21i}^{n+1}, \\ \frac{\overline{\mathbf{J}}_{b_i}^{n+3/2} - \overline{\mathbf{J}}_{b_i}^{n+1/2}}{\Delta t} &= -\Gamma_b \overline{\mathbf{J}}_{b_i}^{[n+1/2]} - \Omega_b^2 \overline{\mathbf{P}}_{b_i}^{[n+1/2]} - K_b \overline{\mathbf{E}} \overline{\mathbf{N}}_{30i}^{n+1}, \end{aligned} \right. \quad (\text{PE}) $$

$$ \begin{aligned} \text{with } \forall v \in \{x, y, z\}, \quad \overline{EN}_{21i}^{v,n+1} &= \left( \left( N_{2ij}^{[n+1/2]} - N_{1ij}^{[n+1/2]} \right) E_{ij}^{v,n+1} \right)_{1 \le j \le d_i} \\ \text{and } \overline{EN}_{30i}^{v,n+1} &= \left( \left( N_{3ij}^{[n+1/2]} - N_{0ij}^{[n+1/2]} \right) E_{ij}^{v,n+1} \right)_{1 \le j \le d_i}. \end{aligned} $$

$$ \left\{ \begin{aligned} \frac{\overline{N}_{3i}^{n+3/2} - \overline{N}_{3i}^{n+1/2}}{\Delta t} &= -\frac{\overline{N}_{3i}^{[n+1/2]}}{l_{32}} + r_b \sum_{v \in \{x,y,z\}} \overline{E}_{b_i}^{v,n+1}, \\ \frac{\overline{N}_{2i}^{n+3/2} - \overline{N}_{2i}^{n+1/2}}{\Delta t} &= \frac{\overline{N}_{3i}^{[n+1/2]}}{l_{32}} - \frac{\overline{N}_{2i}^{[n+1/2]}}{l_{21}} + r_a \sum_{v \in \{x,y,z\}} \overline{E}_{a_i}^{v,n+1}, \\ \frac{\overline{N}_{1i}^{n+3/2} - \overline{N}_{1i}^{n+1/2}}{\Delta t} &= \frac{\overline{N}_{2i}^{[n+1/2]}}{l_{21}} - \frac{\overline{N}_{1i}^{[n+1/2]}}{l_{10}} - r_a \sum_{v \in \{x,y,z\}} \overline{E}_{a_i}^{v,n+1}, \\ \frac{\overline{N}_{0i}^{n+3/2} - \overline{N}_{0i}^{n+1/2}}{\Delta t} &= \frac{\overline{N}_{1i}^{[n+1/2]}}{l_{10}} - r_b \sum_{v \in \{x,y,z\}} \overline{E}_{b_i}^{v,n+1}, \end{aligned} \right. \quad (\text{RE}) $$

$$ \text{with } \overline{E}_{a_i}^{v,n} = \left( \overline{E}_{ij}^{v,n+1} \overline{J}_{a_{ij}}^{v,[n+1/2]} \right)_{1 \le j \le d_i} \quad \text{and} \quad \overline{E}_{b_i}^{v,n} = \left( \overline{E}_{ij}^{v,n+1} \overline{J}_{b_{ij}}^{v,[n+1/2]} \right)_{1 \le j \le d_i}. $$

Therefore, it is necessary to define the boundary conditions in the discrete case, for the restrictions
of the electromagnetic field components. As it was done in [8], for $a_{il} \subset \partial\Omega_h$ we consider the discrete PEC boundary condition

$$ \left\{ \begin{array}{l} \mathbf{E}_i^n|_{a_{il}} = -\mathbf{E}_l^n|_{a_{il}}, \\ \mathbf{H}_i^{n+\frac{1}{2}}|_{a_{il}} = \mathbf{H}_l^{n+\frac{1}{2}}|_{a_{il}}, \end{array} \right. \qquad (4) $$

and the discrete SM boundary condition

$$ \left\{ \begin{array}{l} \mathbf{H}_l^{n+\frac{1}{2}}|_{a_{il}} = -\sqrt{\frac{\varepsilon_{r2}}{\mu_{r2}}} \mathbf{E}_i^{[n+\frac{1}{2}]}|_{a_{il}} \times \mathbf{n}_{il}, \\ \mathbf{E}_l^{n+1}|_{a_{il}} = \sqrt{\frac{\mu_{r2}}{\varepsilon_{r2}}} \mathbf{H}_i^{[n+\frac{1}{2}]}|_{a_{il}} \times \mathbf{n}_{il}. \end{array} \right. \qquad (5) $$

## 3.3 Fully discrete energy

Here, we aim to define energy of the discrete scheme in a way that is consistent with the continuous setting. The goal is to achieve an approximation similar to that derived in Section 3.1.2.

First, for a $T_i$ we denote $(\cdot, \cdot)_{T_i}$ the $L_2$ scalar product on $T_i$ and $\|\cdot\|_{T_i}$ its associated norm, for $n \in \mathbb{N}$, we define the fully discrete energy on $T_i$ at the time step $n\Delta t$ such as, if $T_i \subset \Omega_{1h}$,

$$ \mathcal{E}_i^n = \frac{1}{2} \left( \varepsilon_{r1} \|\mathbf{E}_i^n\|_{T_i}^2 + \mu_{r1} (\mathbf{H}_i^{n+\frac{1}{2}}, \mathbf{H}^{n-1/2})_{T_i} \right. \\ + \frac{1}{K_a} \|\mathbf{J}_{a_i}^{n+\frac{1}{2}}\|_{T_i}^2 + \frac{1}{K_b} \|\mathbf{J}_{b_i}^{n+\frac{1}{2}}\|_{T_i}^2 + \frac{\Omega_a^2}{K_a} \|\mathbf{P}_{a_i}^{n+\frac{1}{2}}\|_{T_i}^2 + \frac{\Omega_b^2}{K_b} + \|\mathbf{P}_{b_i}^{n+\frac{1}{2}}\|_{T_i}^2 \\ \left. + \frac{1}{r_b} \|N_{3_i}^{n+\frac{1}{2}}\|_{T_i}^2 + \frac{1}{r_a} \|N_{2_i}^{n+\frac{1}{2}}\|_{T_i}^2 + \frac{1}{r_a} \|N_{1_i}^{n+\frac{1}{2}}\|_{T_i}^2 + \frac{1}{r_b} \|N_{0_i}^{n+\frac{1}{2}}\|_{T_i}^2 \right), $$

if $T_i \subset \Omega_{2h}$,

$$ \mathcal{E}_i^n = \frac{1}{2} \left( \varepsilon_{r2} \|\mathbf{E}_i^n\|_{T_i}^2 + \mu_{r2} (\mathbf{H}_i^{n+\frac{1}{2}}, \mathbf{H}^{n-1/2})_{T_i} \right). $$

We also define the total discrete energy on $\Omega_h$ such as $\mathcal{E}_{\Omega_h}^n = \sum_{i=1}^N \mathcal{E}_i^n$.

In [7] and [8], the authors show that under a condition on $\Delta t$ depending of the size of each $T_i$, the electromagnetic energy of the system is positive. This relation is known as the CFL condition and is necessary to ensure the consistency of the system. In our work, the added terms of our energy, depending on the unknowns of the polarization equations and rate equations, are positive and do not affect the positivity of the energy. In the following proposition, we present the CFL condition of our problem, extended to $\Omega_h$.

**Proposition 3.3.1 (Positivity of the energy and CFL condition).** Under PEC boundary condition defined at (4), there exists a generic constant $K > 0$ such as, the discrete energy $\mathcal{E}_{\Omega_h}^n$ defines a quadratic form on the unknown, if

$$ \frac{\Delta t}{\sqrt{\max(\varepsilon_{r1}\mu_{r1}, \varepsilon_{r2}\mu_{r2})}} \leq Kh. $$
**Theorem 3.3.1.** Let $n \in \mathbb{N}$, C the constant defined in Theorem 3.1.1 and $\rho$ defined by

$$ \rho = \frac{1}{l_{32}} \sqrt{\frac{r_b}{r_a}} + \frac{1}{l_{21}} + \frac{1}{l_{10}} \sqrt{\frac{r_a}{r_b}}, $$

under PEC boundary condition and the CFL condition of Proposition 3.3.1, if we further assume that

$$ \Delta t \le \frac{1}{2 \left( \sqrt{\frac{K_a + K_b}{2\varepsilon_{r1}}} + \rho \right)}, $$

the total discrete energy at stage n satisfies the following estimate

$$ \mathcal{E}_{\Omega_h}^{n} \leq (1 + 2^{5/2} C \Delta t)^n \mathcal{E}_{\Omega_h}^{0}. \qquad (6) $$

*Proof.* First, we are going to evaluate the variation of the energy at level $n \in \mathbb{N}$, in $T_i \subset \Omega_{1h}$,

$$
\begin{aligned}
\mathcal{E}_i^{n+1} - \mathcal{E}_i^n = & \frac{1}{2}\varepsilon_{r1} (\|\mathbf{E}_i^{n+1}\|_{T_i} - \|\mathbf{E}_i^n\|_{T_i}^2) + \mu_{r1} \left( (\mathbf{H}_i^{n+3/2}, \mathbf{H}_i^{n+1/2})_{T_i} - (\mathbf{H}_i^{n+1/2}, \mathbf{H}_i^{n-1/2})_{T_i} \right) \\
& + \frac{1}{2K_a} \left( \|\mathbf{J}_{a_i}^{n+3/2}\|_{T_i}^2 - \|\mathbf{J}_{a_i}^{n+1/2}\|_{T_i}^2 \right) + \frac{1}{2K_b} \left( \|\mathbf{J}_{b_i}^{n+3/2}\|_{T_i}^2 - \|\mathbf{J}_{b_i}^{n+1/2}\|_{T_i}^2 \right) \\
& + \frac{\Omega_a^2}{2K_a} \left( \|\mathbf{P}_{a_i}^{n+3/2}\|_{T_i}^2 - \|\mathbf{P}_{a_i}^{n+1/2}\|_{T_i}^2 \right) + \frac{\Omega_b^2}{2K_b} \left( \|\mathbf{P}_{b_i}^{n+3/2}\|_{T_i}^2 - \|\mathbf{P}_{b_i}^{n+1/2}\|_{T_i}^2 \right) \\
& + \frac{1}{2r_b} \left( \left( \|N_{3_i}^{n+3/2}\|_{T_i}^2 - \|N_{3_i}^{n+1/2}\|_{T_i}^2 \right) + \left( \|N_{0_i}^{n+3/2}\|_{T_i}^2 - \|N_{0_i}^{n+1/2}\|_{T_i}^2 \right) \right) \\
& + \frac{1}{2r_a} \left( \left( \|N_{2_i}^{n+3/2}\|_{T_i}^2 - \|N_{2_i}^{n+1/2}\|_{T_i}^2 \right) + \left( \|N_{1_i}^{n+3/2}\|_{T_i}^2 - \|N_{1_i}^{n+1/2}\|_{T_i}^2 \right) \right) \\
= & \varepsilon_{r1} (\mathbf{E}_i^{n+1} - \mathbf{E}_i^n, \mathbf{E}_i^{[n+1/2]})_{T_i} + \frac{1}{2}\mu_{r1} (\mathbf{H}_i^{n+1/2}, \mathbf{H}_i^{n+3/2} - \mathbf{H}_i^{n-1/2})_{T_i} \\
& + \frac{1}{K_a} (\mathbf{J}_{a_i}^{n+3/2} - \mathbf{J}_{a_i}^{n+1/2}, \mathbf{J}_{a_i}^{[n+1/2]})_{T_i} + \frac{1}{K_b} (\mathbf{J}_{b_i}^{n+3/2} - \mathbf{J}_{b_i}^{n+1/2}, \mathbf{J}_{b_i}^{[n+1/2]})_{T_i} \\
& + \frac{\Omega_a^2}{K_a} (\mathbf{P}_{a_i}^{n+3/2} - \mathbf{P}_{a_i}^{n+1/2}, \mathbf{P}_{a_i}^{[n+1/2]})_{T_i} + \frac{\Omega_b^2}{K_b} (\mathbf{P}_{b_i}^{n+3/2} - \mathbf{P}_{b_i}^{n+1/2}, \mathbf{P}_{b_i}^{[n+1/2]})_{T_i} \\
& + \frac{1}{r_b} (N_{3_i}^{n+3/2} - N_{3_i}^{n+1/2}, N_{3_i}^{[n+1/2]})_{T_i} + (N_{0_i}^{n+3/2} - N_{0_i}^{n+1/2}, N_{0_i}^{[n+1/2]})_{T_i} \\
& + \frac{1}{r_a} (N_{2_i}^{n+3/2} - N_{2_i}^{n+1/2}, N_{2_i}^{[n+1/2]})_{T_i} + (N_{1_i}^{n+3/2} - N_{1_i}^{n+1/2}, N_{1_i}^{[n+1/2]})_{T_i}.
\end{aligned}
$$

Then, we are going to evaluate

$$ L^n(\mathbf{E}_i, \mathbf{H}_i) = \varepsilon_{r1} (\mathbf{E}_i^{n+1} - \mathbf{E}_i^n, \mathbf{E}_i^{[n+1/2]})_{T_i} + \frac{1}{2}\mu_{r1} (\mathbf{H}_i^{n+1/2}, \mathbf{H}_i^{n+3/2} - \mathbf{H}_i^{n-1/2})_{T_i}, $$

from (ME) at level n, we find

$$ \left( \varepsilon_{r1} \frac{\mathbf{E}_i^{n+1} - \mathbf{E}_i^n}{\Delta t}, \mathbf{E}_i^{[n+1/2]} \right)_{T_i} = \int_{T_i} \varepsilon_{r1} \frac{\mathbf{E}_i^{n+1} - \mathbf{E}_i^n}{\Delta t} \cdot \mathbf{E}_i^{[n+1/2]} $$
$$ = \int_{T_i} \mathbf{H}_i^{n+\frac{1}{2}} \cdot \nabla \times \mathbf{E}_i^{[n+\frac{1}{2}]} - \sum_{l \in \mathcal{V}_i} \int_{a_{il}} \{\mathbf{H}^{n+\frac{1}{2}}\}_{il} \times \mathbf{n}_{il} \cdot \mathbf{E}_i^{[n+\frac{1}{2}]} \\ - \int_{T_i} (\mathbf{J}_{a_i}^{n+\frac{1}{2}} + \mathbf{J}_{b_i}^{n+\frac{1}{2}}) \cdot \mathbf{E}_i^{[n+\frac{1}{2}]}, $$

and

$$ \frac{1}{2} \left( \mu_{r1} \frac{\mathbf{H}_i^{n+3/2} - \mathbf{H}_i^{n-1/2}}{\Delta t}, \mathbf{H}_i^{n+\frac{1}{2}} \right)_{T_i} = \frac{1}{2} \int_{T_i} \mu_r \frac{\mathbf{H}_i^{n+3/2} - \mathbf{H}_i^{n-1/2}}{\Delta t} \cdot \mathbf{H}_i^{n+\frac{1}{2}} \\ = - \int_{T_i} \mathbf{E}_i^{[n+\frac{1}{2}]} \cdot \nabla \times \mathbf{H}_i^{n+\frac{1}{2}} \\ + \sum_{l \in \mathcal{V}_i} \int_{a_{il}} \{\mathbf{E}^{[n+\frac{1}{2}]}\}_{il} \times \mathbf{n}_{il} \cdot \mathbf{H}_i^{n+\frac{1}{2}}. $$

As it is shown in [8], by integrating by part and reorganizing the terms, we have that

$$ L^n(\mathbf{E}_i, \mathbf{H}_i) = - \frac{\Delta t}{2} \sum_{l \in \mathcal{V}_i} \int_{a_{il}} \left( \mathbf{E}_i^{[n+\frac{1}{2}]} \times \mathbf{H}_i^{n+\frac{1}{2}} + \mathbf{E}_l^{[n+\frac{1}{2}]} \times \mathbf{H}_i^{n+\frac{1}{2}} \right) \cdot \mathbf{n}_{il} \\ - \Delta t \int_{T_i} (\mathbf{J}_{a_i}^{n+\frac{1}{2}} + \mathbf{J}_{b_i}^{n+\frac{1}{2}}) \cdot \mathbf{E}_i^{[n+\frac{1}{2}]} $$

Finally

$$ \begin{aligned} \mathcal{E}_i^{n+1} - \mathcal{E}_i^n = & - \frac{\Delta t}{2} \sum_{l \in \mathcal{V}_i} \int_{a_{il}} \left( \mathbf{E}_i^{[n+\frac{1}{2}]} \times \mathbf{H}_i^{n+\frac{1}{2}} + \mathbf{E}_l^{[n+\frac{1}{2}]} \times \mathbf{H}_i^{n+\frac{1}{2}} \right) \cdot \mathbf{n}_{il} \\ & - \Delta t \left( \mathbf{J}_{a_i}^{n+\frac{1}{2}} + \mathbf{J}_{b_i}^{n+\frac{1}{2}}, \mathbf{E}_i^{[n+\frac{1}{2}]} \right)_{T_i} \\ & + \frac{1}{K_a} \left( \mathbf{J}_{a_i}^{n+\frac{3}{2}} - \mathbf{J}_{a_i}^{n+\frac{1}{2}}, \mathbf{J}_{a_i}^{[n+\frac{1}{2}]} \right)_{T_i} + \frac{1}{K_b} \left( \mathbf{J}_{b_i}^{n+\frac{3}{2}} - \mathbf{J}_{b_i}^{n+\frac{1}{2}}, \mathbf{J}_{b_i}^{[n+\frac{1}{2}]} \right)_{T_i} \\ & + \frac{\Omega_a^2}{K_a} \left( \mathbf{P}_{a_i}^{n+\frac{3}{2}} - \mathbf{P}_{a_i}^{n+\frac{1}{2}}, \mathbf{P}_{a_i}^{[n+\frac{1}{2}]} \right)_{T_i} + \frac{\Omega_b^2}{K_b} \left( \mathbf{P}_{b_i}^{n+\frac{3}{2}} - \mathbf{P}_{b_i}^{n+\frac{1}{2}}, \mathbf{P}_{b_i}^{[n+\frac{1}{2}]} \right)_{T_i} \\ & + \frac{1}{r_b} \left( N_{3_i}^{n+\frac{3}{2}} - N_{3_i}^{n+\frac{1}{2}}, N_{3_i}^{[n+\frac{1}{2}]} \right)_{T_i} + \left( N_{0_i}^{n+\frac{3}{2}} - N_{0_i}^{n+\frac{1}{2}}, N_{0_i}^{[n+\frac{1}{2}]} \right)_{T_i} \\ & + \frac{1}{r_a} \left( N_{2_i}^{n+\frac{3}{2}} - N_{2_i}^{n+\frac{1}{2}}, N_{2_i}^{[n+\frac{1}{2}]} \right)_{T_i} + \left( N_{1_i}^{n+\frac{3}{2}} - N_{1_i}^{n+\frac{1}{2}}, N_{1_i}^{[n+\frac{1}{2}]} \right)_{T_i} \end{aligned} $$

With (PE) and (RE) at level n and by noticing the same simplification than in the proof of
Theorem 3.1.1, we obtain for $T_i \subset \Omega_{1h}$

$$
\begin{aligned}
\mathcal{E}_i^{n+1} - \mathcal{E}_i^n = & -\frac{\Delta t}{2} \sum_{l \in \mathcal{V}_i} \int_{a_{il}} (\mathbf{E}_i^{[n+\frac{1}{2}]} \times \mathbf{H}_l^{n+\frac{1}{2}} + \mathbf{E}_l^{[n+\frac{1}{2}]} \times \mathbf{H}_i^{n+\frac{1}{2}}) \cdot \mathbf{n}_{il} \\
& - \Delta t (\mathbf{J}_{a_i}^{n+\frac{1}{2}} + \mathbf{J}_{b_i}^{n+\frac{1}{2}}, \mathbf{E}_i^{[n+\frac{1}{2}]})_{T_i} \\
& - \Delta t \left(\frac{\Gamma_a}{K_a} \|\mathbf{J}_{a_i}^{[n+\frac{1}{2}]}\|_{T_i}^2 + \frac{\Gamma_b}{K_b} \|\mathbf{J}_{b_i}^{[n+\frac{1}{2}]}\|_{T_i}^2\right) \\
& - \Delta t \left(\frac{1}{r_b l_{32}} \|N_{3_i}^{[n+\frac{1}{2}]}\|_{T_i}^2 + \frac{1}{r_a l_{21}} \|N_{2_i}^{[n+\frac{1}{2}]}\|_{T_i}^2 + \frac{1}{r_a l_{10}} \|N_{1_i}^{[n+\frac{1}{2}]}\|_{T_i}^2\right) \\
& + \Delta t \left(\frac{1}{r_a l_{32}} (N_{3_i}^{[n+\frac{1}{2}]}, N_{2_i}^{[n+\frac{1}{2}]})_{T_i} + \frac{1}{r_a l_{21}} (N_{2_i}^{[n+\frac{1}{2}]}, N_{1_i}^{[n+\frac{1}{2}]})_{T_i}\right. \\
& \left. + \frac{1}{r_b l_{10}} (N_{1_i}^{[n+\frac{1}{2}]}, N_{0_i}^{[n+\frac{1}{2}]})_{T_i}\right).
\end{aligned}
$$

For $T_i \subset \Omega_{2h}$, by the same considerations for the Maxwell part of the discrete energy, we obtain

$$
\mathcal{E}_i^{n+1} - \mathcal{E}_i^n = - \frac{\Delta t}{2} \sum_{l \in \mathcal{V}_i} \int_{a_{il}} \left( \mathbf{E}_i^{[n+\frac{1}{2}]} \times \mathbf{H}_l^{n+\frac{1}{2}} + \mathbf{E}_l^{[n+\frac{1}{2}]} \times \mathbf{H}_i^{n+\frac{1}{2}} \right) \cdot \mathbf{n}_{il}.
$$

Now we are going to study the energy defined in all the domain $\Omega_h = \Omega_{1h} \cup \Omega_{2h}$ such as
$\mathcal{E}_{\Omega_h}^n = \sum_{i=1}^N \mathcal{E}_i^n$.

By summing the previous equations $\forall i \in [1, N]$, all the internal boundary terms are simplified
and we obtain for the variation of the discrete energy

$$
\begin{aligned}
\mathcal{E}_{\Omega_h}^{n+1} - \mathcal{E}_{\Omega_h}^{n} = & -\frac{\Delta t}{2} \sum_{a_{il} \subset \partial \Omega_h} \int_{a_{il}} (\mathbf{E}_i^{[n+\frac{1}{2}]} \times \mathbf{H}_l^{n+\frac{1}{2}} + \mathbf{E}_l^{[n+\frac{1}{2}]} \times \mathbf{H}_i^{n+\frac{1}{2}}) \cdot \mathbf{n}_{il} \\
& - \Delta t (\mathbf{J}_{a_h}^{n+\frac{1}{2}} + \mathbf{J}_{b_h}^{n+\frac{1}{2}}, \mathbf{E}_h^{[n+\frac{1}{2}]})_{\Omega_{1h}} \\
& - \Delta t \left(\frac{\Gamma_a}{K_a} \|\mathbf{J}_{a_h}^{[n+\frac{1}{2}]}\|_{\Omega_{1h}}^2 + \frac{\Gamma_b}{K_b} \|\mathbf{J}_{b_h}^{[n+\frac{1}{2}]}\|_{\Omega_{1h}}^2\right) \\
& - \Delta t \left(\frac{1}{r_b l_{32}} \|N_{3_h}^{[n+\frac{1}{2}]}\|_{\Omega_{1h}}^2 + \frac{1}{r_a l_{21}} \|N_{2_h}^{[n+\frac{1}{2}]}\|_{\Omega_{1h}}^2 + \frac{1}{r_a l_{10}} \|N_{1_h}^{[n+\frac{1}{2}]}\|_{\Omega_{1h}}^2\right) \\
& + \Delta t \left(\frac{1}{r_a l_{32}} (N_{3_h}^{[n+\frac{1}{2}]}, N_{2_h}^{[n+\frac{1}{2}]})_{\Omega_{1h}}\right. \\
& \left. + \frac{1}{r_a l_{21}} (N_{2_h}^{[n+\frac{1}{2}]}, N_{1_h}^{[n+\frac{1}{2}]})_{\Omega_{1h}} + \frac{1}{r_b l_{10}} (N_{1_h}^{[n+\frac{1}{2}]}, N_{0_h}^{[n+\frac{1}{2}]})_{\Omega_{1h}}\right),
\end{aligned}
$$

that yields

$$
\begin{aligned}
\mathcal{E}_{\Omega_h}^{n+1} - \mathcal{E}_{\Omega_h}^{n} \leq & -\frac{\Delta t}{2} \sum_{a_{il} \subset \partial \Omega_h} \int_{a_{il}} (\mathbf{E}_i^{[n+\frac{1}{2}]} \times \mathbf{H}_l^{n+\frac{1}{2}} + \mathbf{E}_l^{[n+\frac{1}{2}]} \times \mathbf{H}_i^{n+\frac{1}{2}}) \cdot \mathbf{n}_{il} \\
& - \Delta t (\mathbf{J}_{a_h}^{n+\frac{1}{2}} + \mathbf{J}_{b_h}^{n+\frac{1}{2}}, \mathbf{E}_h^{[n+\frac{1}{2}]})_{\Omega_{1h}} \\
& + \Delta t \left(\frac{1}{r_a l_{32}} (N_{3_h}^{[n+\frac{1}{2}]}, N_{2_h}^{[n+\frac{1}{2}]})_{\Omega_{1h}}\right. \\
& \left. + \frac{1}{r_a l_{21}} (N_{2_h}^{[n+\frac{1}{2}]}, N_{1_h}^{[n+\frac{1}{2}]})_{\Omega_{1h}} + \frac{1}{r_b l_{10}} (N_{1_h}^{[n+\frac{1}{2}]}, N_{0_h}^{[n+\frac{1}{2}]})_{\Omega_{1h}}\right).
\end{aligned}
$$
We can estimate the boundary terms for both boundary conditions. For PEC boundary condition, we have as defined at (4), for $a_{il} \subset \partial\Omega_h$,

$$ \mathbf{E}_{i|a_{il}}^{n} = -\mathbf{E}_{l|a_{il}}^{n}, \quad \mathbf{E}_{i}^{n+1|a_{il}} = -\mathbf{E}_{l}^{n+1|a_{il}} \quad \text{and} \quad \mathbf{H}_{l}^{n+\frac{1}{2}|a_{il}} = \mathbf{H}_{i}^{n+\frac{1}{2}|a_{il}} \quad \text{so} \\ \mathbf{E}_{i}^{[n+\frac{1}{2}]}|_{a_{il}} \times \mathbf{H}_{l}^{n+\frac{1}{2}|a_{il}} = -\mathbf{E}_{i}^{[n+\frac{1}{2}]}|_{a_{il}} \times \mathbf{H}_{i}^{n+\frac{1}{2}|a_{il}} $$

and finally,

$$ \int_{a_{il}} (\mathbf{E}_{i}^{[n+\frac{1}{2}]} \times \mathbf{H}_{l}^{n+\frac{1}{2}} + \mathbf{E}_{l}^{[n+\frac{1}{2}]} \times \mathbf{H}_{i}^{n+\frac{1}{2}}) \cdot \mathbf{n}_{il} = 0, \forall a_{il} \subset \partial\Omega_{h}. $$

We have now that

$$ \mathcal{E}_{\Omega_h}^{n+1} - \mathcal{E}_{\Omega_h}^{n} \leq -\Delta t (\mathbf{J}_{a_h}^{n+\frac{1}{2}} + \mathbf{J}_{b_h}^{n+\frac{1}{2}}, \mathbf{E}_{h}^{[n+\frac{1}{2}]})_{\Omega_{1h}} \\ + \Delta t \left( \frac{1}{r_a l_{32}} (N_{3h}^{[n+\frac{1}{2}]}, N_{2h}^{[n+\frac{1}{2}]})_{\Omega_{1h}} \right. \\ \left. + \frac{1}{r_a l_{21}} (N_{2h}^{[n+\frac{1}{2}]}, N_{1h}^{[n+\frac{1}{2}]})_{\Omega_{1h}} + \frac{1}{r_b l_{10}} (N_{1h}^{[n+\frac{1}{2}]}, N_{0h}^{[n+\frac{1}{2}]})_{\Omega_{1h}} \right). $$

With Cauchy-Schwarz inequality we have

$$ \mathcal{E}_{\Omega_h}^{n+1} - \mathcal{E}_{\Omega_h}^{n} \leq \frac{\Delta t}{2} \left\| \mathbf{J}_{a_h}^{n+\frac{1}{2}} + \mathbf{J}_{b_h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}} \left\| \mathbf{E}_h^{n+1} + \mathbf{E}_h^n \right\|_{\Omega_{1h}} \\ + \frac{\Delta t}{4} \left( \frac{1}{r_a l_{32}} \left\| N_{3h}^{n+\frac{3}{2}} + N_{3h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}} \left\| N_{2h}^{n+\frac{3}{2}} + N_{2h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}} \right. \\ \left. + \frac{1}{r_a l_{21}} \left\| N_{2h}^{n+\frac{3}{2}} + N_{2h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}} \left\| N_{1h}^{n+\frac{3}{2}} + N_{1h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}} \right. \\ \left. + \frac{1}{r_b l_{10}} \left\| N_{1h}^{n+\frac{3}{2}} + N_{1h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}} \left\| N_{0h}^{n+\frac{3}{2}} + N_{0h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}} \right). $$

With Young inequality, we get that for $\omega, \omega_1, \omega_2, \omega_3 > 0$

$$ \left\| \mathbf{J}_{a_h}^{n+\frac{1}{2}} + \mathbf{J}_{b_h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}} \left\| \mathbf{E}_h^{n+1} + \mathbf{E}_h^n \right\|_{\Omega_{1h}} \leq \frac{1}{2\omega} \left\| \mathbf{J}_{a_h}^{n+\frac{1}{2}} + \mathbf{J}_{b_h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}}^2 + \frac{\omega}{2} \left\| \mathbf{E}_h^{n+1} + \mathbf{E}_h^n \right\|_{\Omega_{1h}}^2, $$

$$ \left\| N_{3h}^{n+\frac{3}{2}} + N_{3h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}} \left\| N_{2h}^{n+\frac{3}{2}} + N_{2h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}} \leq \frac{1}{2\omega_1} \left\| N_{3h}^{n+\frac{3}{2}} + N_{3h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}}^2 + \frac{\omega_1}{2} \left\| N_{2h}^{n+\frac{3}{2}} + N_{2h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}}^2, $$

$$ \left\| N_{1h}^{n+\frac{3}{2}} + N_{1h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}} \left\| N_{2h}^{n+\frac{3}{2}} + N_{2h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}} \leq \frac{1}{2\omega_2} \left\| N_{1h}^{n+\frac{3}{2}} + N_{1h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}}^2 + \frac{\omega_2}{2} \left\| N_{2h}^{n+\frac{3}{2}} + N_{2h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}}^2, $$

and

$$ \left\| N_{1h}^{n+\frac{3}{2}} + N_{1h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}} \left\| N_{0h}^{n+\frac{3}{2}} + N_{0h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}} \leq \frac{1}{2\omega_3} \left\| N_{1h}^{n+\frac{3}{2}} + N_{1h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}}^2 + \frac{\omega_3}{2} \left\| N_{0h}^{n+\frac{3}{2}} + N_{0h}^{n+\frac{1}{2}} \right\|_{\Omega_{1h}}^2. $$
So we can have, for $\omega, \omega_1, \omega_2, \omega_3 > 0$

$$
\begin{align*}
\mathcal{E}_{\Omega_h}^{n+1} - \mathcal{E}_{\Omega_h}^n &\le \frac{\Delta t}{2} \left( \frac{1}{2\omega} \|\mathbf{J}_{a_h}^{n+\frac{1}{2}} + \mathbf{J}_{b_h}^{n+\frac{1}{2}}\|_{\Omega_{1h}}^2 + \frac{\omega}{2} \|\mathbf{E}_h^{n+1} + \mathbf{E}_h^n\|_{\Omega_{1h}}^2 \right) \\
&\quad + \frac{\Delta t}{4} \left( \frac{1}{r_a l_{32}} \left( \frac{1}{2\omega_1} \|N_{3h}^{n+\frac{3}{2}} + N_{3h}^{n+\frac{1}{2}}\|_{\Omega_{1h}}^2 + \frac{\omega_1}{2} \|N_{2h}^{n+\frac{3}{2}} + N_{3h}^{n+\frac{1}{2}}\|_{\Omega_{1h}}^2 \right) \right. \\
&\qquad \left. + \frac{1}{r_a l_{21}} \left( \frac{1}{2\omega_2} \|N_{1h}^{n+\frac{3}{2}} + N_{1h}^{n+\frac{1}{2}}\|_{\Omega_{1h}}^2 + \frac{\omega_2}{2} \|N_{2h}^{n+\frac{3}{2}} + N_{3h}^{n+\frac{1}{2}}\|_{\Omega_{1h}}^2 \right) \right. \\
&\qquad \left. + \frac{1}{r_b l_{10}} \left( \frac{1}{2\omega_3} \|N_{1h}^{n+\frac{3}{2}} + N_{1h}^{n+\frac{1}{2}}\|_{\Omega_{1h}}^2 + \frac{\omega_3}{2} \|N_{0h}^{n+\frac{3}{2}} + N_{0h}^{n+\frac{1}{2}}\|_{\Omega_{1h}}^2 \right) \right) \\
&\le \frac{\Delta t}{2} \left( \frac{1}{\omega} \|\mathbf{J}_{a_h}^{n+\frac{1}{2}}\|_{\Omega_{1h}}^2 + \frac{1}{\omega} \|\mathbf{J}_{b_h}^{n+\frac{1}{2}}\|_{\Omega_{1h}}^2 + \omega \|\mathbf{E}_h^{n+1}\|_{\Omega_{1h}}^2 + \omega \|\mathbf{E}_h^n\|_{\Omega_{1h}}^2 \right) \\
&\quad + \frac{\Delta t}{4} \left( \frac{1}{r_a l_{32}} \left( \frac{1}{\omega_1} \left( \|N_{3h}^{n+\frac{3}{2}}\|_{\Omega_{1h}}^2 + \|N_{3h}^{n+\frac{1}{2}}\|_{\Omega_{1h}}^2 \right) + \omega_1 \left( \|N_{2h}^{n+\frac{3}{2}}\|_{\Omega_{1h}}^2 + \|N_{2h}^{n+\frac{1}{2}}\|_{\Omega_{1h}}^2 \right) \right) \right. \\
&\qquad \left. + \frac{1}{r_a l_{21}} \left( \frac{1}{\omega_2} \left( \|N_{1h}^{n+\frac{3}{2}}\|_{\Omega_{1h}}^2 + \|N_{1h}^{n+\frac{1}{2}}\|_{\Omega_{1h}}^2 \right) + \omega_2 \left( \|N_{2h}^{n+\frac{3}{2}}\|_{\Omega_{1h}}^2 + \|N_{2h}^{n+\frac{1}{2}}\|_{\Omega_{1h}}^2 \right) \right) \right. \\
&\qquad \left. + \frac{1}{r_b l_{10}} \left( \frac{1}{\omega_3} \left( \|N_{1h}^{n+\frac{3}{2}}\|_{\Omega_{1h}}^2 + \|N_{1h}^{n+\frac{1}{2}}\|_{\Omega_{1h}}^2 \right) + \omega_3 \left( \|N_{0h}^{n+\frac{3}{2}}\|_{\Omega_{1h}}^2 + \|N_{0h}^{n+\frac{1}{2}}\|_{\Omega_{1h}}^2 \right) \right) \right) \\
&\le \Delta t \left( \frac{1}{\omega} (K_a + K_b) \mathcal{E}_{\Omega_h}^n + \frac{\omega}{\varepsilon_{r1}} (\mathcal{E}_{\Omega_h}^{n+1} + \mathcal{E}_{\Omega_h}^n) \right) \\
&\quad + \frac{\Delta t}{2} \left( \frac{1}{r_a l_{32}} \left( \frac{r_b}{\omega_1} (\mathcal{E}_{\Omega_h}^{n+1} + \mathcal{E}_{\Omega_h}^n) + r_a \omega_1 (\mathcal{E}_{\Omega_h}^{n+1} + \mathcal{E}_{\Omega_h}^n) \right) \right. \\
&\qquad \left. + \frac{1}{r_a l_{21}} \left( \frac{r_a}{\omega_2} (\mathcal{E}_{\Omega_h}^n + \mathcal{E}_{\Omega_h}^{n+1}) + r_a \omega_2 (\mathcal{E}_{\Omega_h}^n + \mathcal{E}_{\Omega_h}^{n+1}) \right) \right. \\
&\qquad \left. + \frac{1}{r_b l_{10}} \left( \frac{r_a}{\omega_3} (\mathcal{E}_{\Omega_h}^n + \mathcal{E}_{\Omega_h}^{n+1}) + r_b \omega_3 (\mathcal{E}_{\Omega_h}^n + \mathcal{E}_{\Omega_h}^{n+1}) \right) \right).
\end{align*}
$$

Finally, we find

$$
\mathcal{E}_{\Omega_h}^{n+1} - \mathcal{E}_{\Omega_h}^n \le \Delta t \left( \frac{1}{\omega} (K_a + K_b) + \frac{\omega}{\varepsilon_{r1}} + \rho(\omega_1, \omega_2, \omega_3) \right) \mathcal{E}_{\Omega_h}^n + \Delta t \left( \frac{\omega}{\varepsilon_{r1}} + \rho(\omega_1, \omega_2, \omega_3) \right) \mathcal{E}_{\Omega_h}^{n+1},
$$

$$
\text{with } \rho(\omega_1, \omega_2, \omega_3) = \frac{1}{2} \left( \frac{r_b}{r_a l_{32} \omega_1} + \frac{\omega_1}{l_{32}} + \frac{1}{l_{21} \omega_2} + \frac{\omega_2}{l_{21}} + \frac{r_a}{r_b l_{10} \omega_3} + \frac{\omega_3}{l_{32}} \right).
$$

We now suppose that

$$
\Delta t < \frac{1}{\frac{\omega}{\varepsilon_{r1}} + \rho(\omega_1, \omega_2, \omega_3)}
$$

we have under this condition, the following inequality for $\mathcal{E}_{\Omega_h}^{n+1}$, with $\omega > 0$

$$
\mathcal{E}_{\Omega_h}^{n+1} \le \frac{1 + \Delta t \left( \frac{1}{\omega} (K_a + K_b) + \frac{\omega}{\varepsilon_{r1}} + \rho(\omega_1, \omega_2, \omega_3) \right)}{1 - \Delta t \left( \frac{\omega}{\varepsilon_{r1}} + \rho(\omega_1, \omega_2, \omega_3) \right)} \mathcal{E}_{\Omega_h}^n.
$$
Let us denote $C(\Delta t) = \frac{1 + \Delta t \left( \frac{1}{\omega}(K_a + K_b) + \frac{\omega}{\varepsilon_{r1}} + \rho(\omega_1, \omega_2, \omega_3) \right)}{1 - \Delta t \left( \frac{\omega}{\varepsilon_{r1}} + \rho(\omega_1, \omega_2, \omega_3) \right)}$.

The previous inequality gives us that

$$ \mathcal{E}_{\Omega_h}^n \le \mathcal{E}_{\Omega_h}^0 C(\Delta t)^n. $$

Under the previous condition on $\Delta t$, we have that $C(\Delta t) > 0$, so we can assume that we reach energy stability.

We can show that, $\forall \eta > 1$, if $\Delta t \le \frac{1 - \frac{1}{\eta}}{\left( \frac{\omega}{\varepsilon_{r1}} + \rho(\omega_1, \omega_2, \omega_3) \right)}$ we have

$$ 0 < C(\Delta t) \le 1 + \eta \left( \frac{1}{\omega}(K_a + K_b) + \frac{2\omega}{\varepsilon_{r1}} + 2\rho(\omega_1, \omega_2, \omega_3) \right) \Delta t. $$

We denote $K(\omega, \omega_1, \omega_2, \omega_3) = \frac{1}{\omega}(K_a + K_b) + \frac{2\omega}{\varepsilon_r} + \rho(\omega_1, \omega_2, \omega_3)$.

$K$ is minimal for

$$ \tilde{\omega} = \sqrt{\frac{\varepsilon_{r1}(K_a + K_b)}{2}}, \quad \tilde{\omega}_1 = \sqrt{\frac{r_b}{r_a}}, \quad \tilde{\omega}_2 = 1, \quad \tilde{\omega}_3 = \sqrt{\frac{r_a}{r_b}}. $$

In this case

$$ \rho = \rho(\tilde{\omega}_1, \tilde{\omega}_2, \tilde{\omega}_3) = \frac{1}{l_{32}} \sqrt{\frac{r_b}{r_a}} + \frac{1}{l_{21}} + \frac{1}{l_{10}} \sqrt{\frac{r_a}{r_b}}. $$

So, $\forall \eta > 1$, under the assumption

$$ \Delta t \le \frac{1 - \frac{1}{\eta}}{\sqrt{\frac{K_a + K_b}{2\varepsilon_{r1}} + \rho}}, $$

we have

$$ \mathcal{E}_{\Omega_h}^n \le \left( 1 + 2\eta \left( \sqrt{2\frac{K_a + K_b}{\varepsilon_{r1}}} + \rho \right) \Delta t \right)^n \mathcal{E}_{\Omega_h}^0. $$

In Theorem 3.1.1, we have defined the following constant, $C = \sqrt{\frac{K_a + K_b}{\varepsilon_{r1}}} + \rho$.

We can see here that

$$ \sqrt{2\frac{K_a + K_b}{\varepsilon_{r1}}} + \rho < \sqrt{2} \left( \sqrt{\frac{K_a + K_b}{\varepsilon_{r1}}} + \rho \right) = \sqrt{2C}, $$

So we have, $\forall \eta > 1$, under the assumption

$$ \Delta t \le \frac{1 - \frac{1}{\eta}}{\sqrt{\frac{K_a + K_b}{2\varepsilon_{r1}} + \rho}}, \text{ we have } \mathcal{E}_{\Omega_h}^n \le \left( 1 + 2^{3/2} C \eta \Delta t \right)^n \mathcal{E}_{\Omega_h}^0. $$
In particular if we take $\eta = 2$, under the assumption

$$ \Delta t \le \frac{1}{2\left(\sqrt{\frac{K_a + K_b}{2\varepsilon_{r1}} + \rho}\right)}, \text{ we have } \mathcal{E}_{\Omega_h}^n \le \left(1 + 2^{5/2}C\Delta t\right)^n \mathcal{E}_{\Omega_h}^0. $$

This concludes the proof of Theorem 3.3.1. □

To estimate the discrete energy of the system for SM boundary condition, we have to correct the electromagnetic part of the energy, by defining a new quantity, as it has been done in [8]

$$ \mathcal{F}_{\Omega_h}^n = \mathcal{E}_{\Omega_h}^n + \frac{\Delta t}{4} \sqrt{\frac{\mu_{r2}}{\varepsilon_{r2}}} \sum_{a_{il} \subset \Omega_h} \int_{a_{il}} \left( \mathbf{n}_{il} \times \mathbf{H}_i^{n-1/2} \right) \cdot \left( \mathbf{n}_{il} \times \frac{\mathbf{H}_i^{n+1/2} + \mathbf{H}_i^{n-1/2}}{2} \right). $$

*Remark.* We can also see in [8], that under SM boundary condition the CFL condition is exactly the same than in Proposition 3.3.1.

Under these considerations, we have the following result for the energy under SM boundary condition on $\partial\Omega_h$.

**Theorem 3.3.2.** Let $n \in \mathbb{N}$, under SM boundary condition and the CFL condition of Proposition 3.3.1, if we further assume that

$$ \Delta t \le \frac{1}{2\left(\sqrt{\frac{K_a + K_b}{2\varepsilon_{r1}} + \rho}\right)} $$

the total discrete energy at stage n satisfies the following estimate

$$ \mathcal{F}_{\Omega_h}^n \le \left(1 + 2^{5/2}C\Delta t\right)^n \mathcal{F}_{\Omega_h}^0. \qquad (7) $$

*Proof.* As it has been shown in [8], the electromagnetic part of the discrete energy with correction is decreasing under Silver Müller boundary condition. The rest of the energy can be treated exactly as for the proof of theorem 3.1.1 for PEC boundary condition. □

# 4 Numerical results

The proposed DGTD method in 3D has been implemented in Fortran 95. Here, we limit ourselves to a validation of this implementation using an artificial test case, physical studies will be presented in another paper.

For our validation tests, we have considered a cubic cavity of length L, such as $\Omega = \Omega_1 = [0, L]^3$, with PEC condition at the boundary. We cannot have simple exact solutions for this problem so we have to build artificial solutions. Here we present how this manufactured solution has been constructed by drawing inspiration from the solution used for the pure 3D Maxwell's equations problem. For the expression of this analytic solution we need to define some functions and parameters, the definitions can be found in appendix A.
## 4.1 Manufactured solutions

In this part, we build an artificial analytic solution for the problem. We start by imposing a solution for the electric and magnetic fields. We chose to take a classical artificial solution for the electromagnetic field for Maxwell problem in a PEC cubic cavity of size L. We can find a similar approach for Drude-Lorenz model in [18]. We respectively denote $\overline{\mathbf{E}}$ and $\overline{\mathbf{H}}$ the artificial electric and magnetic fields

$$ \overline{\mathbf{E}}(\mathbf{x}, t) = \cos(\omega t) \begin{pmatrix} -f(\mathbf{x}) \\ 0 \\ g(\mathbf{x}) \end{pmatrix} \quad \text{and} \quad \overline{\mathbf{H}}(\mathbf{x}, t) = \frac{x_m}{\omega \mu_r} \sin(\omega t) \begin{pmatrix} -h(\mathbf{x}) \\ 2i(\mathbf{x}) \\ -j(\mathbf{x}) \end{pmatrix}, $$

where $x_m = n_m\pi/L$, with $n_m \in \mathbb{N}^*$ the number of the mode and f, g, h, i and j are spatial functions, product of sines and cosines dependent on $x_m$ and defined in appendix A.

$\overline{\mathbf{E}}$ and $\overline{\mathbf{H}}$ defined as followed are solution of

$$ \left\{ \begin{array}{ll} \displaystyle \frac{\partial \overline{\mathbf{E}}}{\partial t} = \mu_r \frac{\omega^2}{3x_m^2} \nabla \times \overline{\mathbf{H}} & \text{in } \Omega_1 \times [0, T], \\ \\ \displaystyle \mu_r \frac{\partial \overline{\mathbf{H}}}{\partial t} = -\nabla \times \overline{\mathbf{E}} & \text{in } \Omega_1 \times [0, T]. \end{array} \right. $$

At this stage, by knowing $\overline{\mathbf{E}}$, we could obtain an explicit expression for the other parameters in the form $t \to \exp(\int_0^t B(s)ds)U_0$. However, this calculation is complex as it requires evaluating a non-trivial matrix exponential. Therefore, we will also consider artificial solutions for $\mathbf{P}_a$, $\mathbf{P}_b$ and their respective derivatives $\mathbf{J}_a$, $\mathbf{J}_b$, which we will note respectively $\overline{\mathbf{P}}_b$, $\overline{\mathbf{P}}_a$, $\overline{\mathbf{J}}_a$ and $\overline{\mathbf{J}}_b$. These two artificial solutions will be the solution of the Drude-Lorentz following problem

$$ \left\{ \begin{array}{ll} \displaystyle \frac{\partial^2 \overline{\mathbf{P}}_a}{\partial t^2} + \Gamma_a \frac{\partial \overline{\mathbf{P}}_a}{\partial t} + \Omega_a^2 \overline{\mathbf{P}}_b = -K_a \overline{\mathbf{E}} & \text{in } \Omega_1 \times [0, T], \\ \\ \displaystyle \frac{\partial^2 \overline{\mathbf{P}}_b}{\partial t^2} + \Gamma_b \frac{\partial \overline{\mathbf{P}}_b}{\partial t} + \Omega_b^2 \overline{\mathbf{P}}_b = -K_b \overline{\mathbf{E}} & \text{in } \Omega_1 \times [0, T]. \end{array} \right. $$

We do not have any particular initial conditions imposed here on $\overline{\mathbf{P}}_b$, $\overline{\mathbf{P}}_a$, and their derivatives. Therefore, we have chosen the initial conditions such that the time solutions of these equations are a linear combination of sines and cosines. This consideration facilitates the resolution of the rate equations. Then, we find the following analytic expression for the polarization vectors, depending on parameters defined in appendix A

$$ \begin{aligned} \overline{\mathbf{P}}_a(\mathbf{x}, t) &= -Z_a ((\Omega_a^2 - \omega^2) \sin(\omega t) + \omega \Gamma_a \cos(\omega t)) \begin{pmatrix} -f(\mathbf{x}) \\ 0 \\ g(\mathbf{x}) \end{pmatrix}, \\ \overline{\mathbf{P}}_b(\mathbf{x}, t) &= -Z_b ((\Omega_b^2 - \omega^2) \sin(\omega t) + \omega \Gamma_b \cos(\omega t)) \begin{pmatrix} -f(\mathbf{x}) \\ 0 \\ g(\mathbf{x}) \end{pmatrix}. \end{aligned} $$

To incorporate these fields in our equations, we need to add an explicit source current $\mathbf{J}_s$ in Maxwell's equations and two other sources $\mathbf{P}_{as}$ and $\mathbf{P}_{bs}$ in polarization equations. $\overline{\mathbf{E}}$, $\overline{\mathbf{H}}$, $\overline{\mathbf{J}}_a$, and $\overline{\mathbf{J}}_b$ being explicitly known at this point, we can deduce an expression for these sources from
the following equations

$$
\left\{
\begin{aligned}
\epsilon_r \frac{\partial \bar{\mathbf{E}}}{\partial t} &= \nabla \times \bar{\mathbf{H}} - (\bar{\mathbf{J}}_a + \bar{\mathbf{J}}_b) + \mathbf{J}_s & \text{in } \Omega_1 \times [0, T], \\
\frac{\partial^2 \bar{\mathbf{P}}_b}{\partial t^2} + \Gamma_a \frac{\partial \bar{\mathbf{P}}_b}{\partial t} + \Omega_a^2 (\bar{\mathbf{P}}_b - \mathbf{P}_{as}) &= -K_a (N_2 - N_1) \bar{\mathbf{E}} & \text{in } \Omega_1 \times [0, T], \\
\frac{\partial^2 \bar{\mathbf{P}}_b}{\partial t^2} + \Gamma_b \frac{\partial \bar{\mathbf{P}}_b}{\partial t} + \Omega_b^2 (\bar{\mathbf{P}}_b - \mathbf{P}_{bs}) &= -K_b (N_3 - N_0) \bar{\mathbf{E}} & \text{in } \Omega_1 \times [0, T].
\end{aligned}
\right.
$$

We obtain the following analytical expression for the sources current **J**<sub>s</sub> and the polarization sources **P**<sub>as</sub> and **P**<sub>bs</sub>

$$
\begin{align*}
\mathbf{J}_s(\mathbf{x}, t) &= \left( \omega \left( 3 \frac{x_m^2}{\mu_r \omega^2} - \epsilon_r + Z_a (\Omega_a^2 - \omega^2) + Z_b (\Omega_b^2 - \omega^2) \right) \sin(\omega t) \right. \\
&\qquad \left. - \omega^2 (Z_a \Gamma_a + Z_b \Gamma_b) \cos(\omega t) \right) \begin{pmatrix} -f(\mathbf{x}) \\ 0 \\ g(\mathbf{x}) \end{pmatrix}, \\
\mathbf{P}_{as}(\mathbf{x}, t) &= -\frac{K_a}{\Omega_a^2} (1 - (N_2(\mathbf{x}, t) - N_1(\mathbf{x}, t))) \cos(\omega t) \begin{pmatrix} -f(\mathbf{x}) \\ 0 \\ g(\mathbf{x}) \end{pmatrix}, \\
\mathbf{P}_{bs}(\mathbf{x}, t) &= -\frac{K_b}{\Omega_b^2} (1 - (N_3(\mathbf{x}, t) - N_0(\mathbf{x}, t))) \cos(\omega t) \begin{pmatrix} -f(\mathbf{x}) \\ 0 \\ g(\mathbf{x}) \end{pmatrix}.
\end{align*}
$$

Then, knowing $\overline{\mathbf{E}}$, $\overline{\mathbf{J}}_a$ and $\overline{\mathbf{J}}_b$, the rate equations system becomes a first-order linear ODE system which can be solved considering that, $\forall \mathbf{x} \in \Omega_1$, $N_0(\mathbf{x}, 0) = 1$, and $N_1(\mathbf{x}, 0) = N_2(\mathbf{x}, 0) = N_3(\mathbf{x}, 0) = 0$. The solutions of these equations are a linear combination of sines, cosines, and exponential, depending on parameters defined in appendix A

$$
\left\{
\begin{aligned}
N_3(\mathbf{x}, t) &= \left( \alpha_3 e^{-t/l_{32}} + \beta_3 \cos(\omega t) \sin(\omega t) + \gamma_3 \cos^2(\omega t) + \delta_3 \right) (f^2(\mathbf{x}) + g^2(\mathbf{x})), \\
N_2(\mathbf{x}, t) &= \left( \alpha_2 e^{-t/l_{21}} + \frac{\alpha_3 l_{21}}{l_{32}-l_{21}} e^{-t/l_{32}} + \beta_2 \cos(\omega t) \sin(\omega t) \right. \\
&\qquad \left. + \gamma_2 \cos^2(\omega t) + \delta_2 \right) (f^2(\mathbf{x}) + g^2(\mathbf{x})), \\
N_1(\mathbf{x}, t) &= \left( \alpha_1 e^{-t/l_{10}} + \frac{\alpha_2 l_{10}}{l_{21}-l_{10}} e^{-t/l_{21}} + \frac{\alpha_3 l_{32} l_{10}}{(l_{32}-l_{21})(l_{32}-l_{10})} e^{-t/l_{32}} \right. \\
&\qquad \left. + \beta_1 \cos(\omega t) \sin(\omega t) + \gamma_1 \cos^2(\omega t) + \delta_1 \right) (f^2(\mathbf{x}) + g^2(\mathbf{x})), \\
N_0(\mathbf{x}, t) &= \left( \alpha_0 t - \alpha_1 e^{-t/l_{10}} - \frac{\alpha_2 l_{21}}{l_{21}-l_{10}} e^{-t/l_{21}} - \frac{\alpha_3 l_{32}^2}{(l_{32}-l_{21})(l_{32}-l_{10})} e^{-t/l_{32}} \right. \\
&\qquad \left. + \beta_0 \cos(\omega t) \sin(\omega t) + \gamma_0 \cos^2(\omega t) + \delta_0 \right) (f^2(\mathbf{x}) + g^2(\mathbf{x})) + 1.
\end{aligned}
\right.
$$

## 4.2 Numerical results

We present here the convergence order results obtained for simulations carried out with $T = 1$ in a cubic cavity of length $L = 1$, with PEC boundary conditions. The simulations were carried out using elements from $\mathbb{P}_1$ to $\mathbb{P}_2$ and various tetrahedral mesh sizes. These meshes are generated from uniform cubic grids, where $h$ represents the size of each cube, and each cube is
subdivided into six uniform tetrahedrons. Let $N$ be the number of cubes along each dimension leading to $h = 1/(N - 1)$. The simulations were conducted for four different mesh sizes, with $N \in \{11, 21, 31, 41\}$. Let $n_s$ and $n_t$ denote the number of vertices and tetrahedrons, respectively. Table 2 summarizes the key parameters for each mesh configuration.

<table><caption>Table 2: Meshes characteristics for the cubic PEC cavity case.</caption><thead><tr><th></th><th>M<sub>1</sub></th><th>M<sub>2</sub></th><th>M<sub>3</sub></th><th>M<sub>4</sub></th></tr></thead><tbody><tr><td>N</td><td>11</td><td>21</td><td>31</td><td>41</td></tr><tr><td>h</td><td>0.1</td><td>0.05</td><td>0.034</td><td>0.025</td></tr><tr><td>n<sub>s</sub></td><td>1331</td><td>9261</td><td>29791</td><td>68921</td></tr><tr><td>n<sub>t</sub></td><td>1600</td><td>48000</td><td>162000</td><td>384000</td></tr></tbody></table>

For a mesh with parameter $h$, we define the system's error by comparing the approximated solution $\mathcal{U}_h^n$ with the exact solution $\mathcal{U}_{ex}^n$ at each time step $n$, using the norm defined in Section 3.1.1. Let $e_{tot}^h$ be the total error for a mesh with parameter $h$, given by

$$e_{tot}^h = \max_n (||\mathcal{U}_h^n - \mathcal{U}_{ex}^n||_\varepsilon).$$

To analyze the influence of every parameters on the total error, we introduce three additional error metrics. Let $e_{E,H}^h$, $e_{pol}^h$, and $e_{elec}^h$ denote, respectively, the electromagnetic error, the polarization error, and the electronic density error. These are defined as follows, based on the parameters of the equation

$$e_{E,H}^h = \max_n \left( \varepsilon_r ||\mathbf{E}_h^n - \mathbf{E}_{ex}^n||_{\mathbf{L}^2(\Omega)}^2 + ||\mathbf{H}_h^{n+1/2} - \mathbf{H}_{ex}^{n+1/2}||_{\mathbf{L}^2(\Omega)}^2 \right)^{\frac{1}{2}},$$

$$e_{pol}^h = \max_n \left( \frac{\Omega_a^2}{K_a} ||\mathbf{P}_{a_h}^{n+1/2} - \mathbf{P}_{a_{ex}}^{n+1/2}||_{\mathbf{L}^2(\Omega)}^2 + \frac{\Omega_b^2}{K_b} ||\mathbf{P}_{b_h}^{n+1/2} - \mathbf{P}_{b_{ex}}^{n+1/2}||_{\mathbf{L}^2(\Omega)}^2 \right)^{\frac{1}{2}},$$

and

$$e_{elec}^h = \max_n \left( \frac{1}{r_b} \left( ||N_{3h}^{n+1/2} - N_{3ex}^{n+1/2}||_{L^2(\Omega)}^2 + ||N_{0h}^{n+1/2} - N_{0ex}^{n+1/2}||_{L^2(\Omega)}^2 \right) \right. \\ \left. + \frac{1}{r_a} \left( ||N_{2h}^{n+1/2} - N_{2ex}^{n+1/2}||_{L^2(\Omega)}^2 + ||N_{1h}^{n+1/2} - N_{1ex}^{n+1/2}||_{L^2(\Omega)}^2 \right) \right)^{\frac{1}{2}}.$$

We used three different sets of parameters for the gain medium model, they are presented in Table 3.

<table><caption>Table 3: Values of the homogenized parameters of the gain media.</caption><thead><tr><th>Constants</th><th>Ω<sub>a</sub></th><th>Ω<sub>b</sub></th><th>Γ<sub>a</sub></th><th>Γ<sub>b</sub></th><th>K<sub>a</sub></th><th>K<sub>b</sub></th><th>r<sub>a</sub></th><th>r<sub>b</sub></th><th>l<sub>10</sub></th><th>l<sub>21</sub></th><th>l<sub>32</sub></th></tr></thead><tbody><tr><td>Set 1</td><td>0.5</td><td>4</td><td>4.</td><td>1.</td><td>1.</td><td>3.</td><td>4.</td><td>1.5</td><td>0.6</td><td>0.2</td><td>0.4</td></tr><tr><td>Set 2</td><td>7.</td><td>9.</td><td>20.</td><td>15.</td><td>0.1</td><td>0.01</td><td>0.4</td><td>0.1</td><td>0.1</td><td>0.2</td><td>0.15</td></tr><tr><td>Set 3</td><td>0.5</td><td>0.3</td><td>0.1</td><td>0.2</td><td>6.</td><td>8.</td><td>9.</td><td>6.</td><td>8.</td><td>9.</td><td>11.</td></tr></tbody></table>
Figure 3 presents the main results for set 1, with a particular focus on convergence analysis. In Figure 3a, a logarithmic-scale graph illustrates the evolution of the total error $e_{tot}^h$ across the four meshes. Figure 3b presents the variation of the total energy at each time step, as defined in 3.1.1, for $M_4$. Finally, Figures 3c and 3d show the error evolution at each time step for each mesh, using the DGTD $P_1$ and $P_2$ methods, respectively.

<figure><img src="image_5.png" /><figcaption>Figure 3: Comparison of error evolution and energy evolution for DGTD $P_1$ and DGTD $P_2$ methods.</figcaption></figure>

We see in Table 4 the different values presented in Figure 3a and the order of convergences for DGTD $P_1$ and $P_2$ methods. As expected we have an order of convergence of 1 for $P_1$ and 2 for $P_2$.
<table><caption>Table 4: Evolution of errors and convergence orders for DGTD $P_1$ and DGTD $P_2$ methods for second order leap frog scheme for set 1.</caption><thead><tr><th rowspan="2">Mesh</th><th colspan="2">P<sub>1</sub></th><th colspan="2">P<sub>2</sub></th></tr><tr><th>c<sup>h</sup><sub>tot</sub></th><th>Order</th><th>c<sup>h</sup><sub>tot</sub></th><th>Order</th></tr></thead><tbody><tr><td>M<sub>1</sub></td><td>2.522 × 10<sup>-2</sup></td><td>-</td><td>1.992 × 10<sup>-3</sup></td><td>-</td></tr><tr><td>M<sub>2</sub></td><td>1.048 × 10<sup>-2</sup></td><td>1.36</td><td>4.539 × 10<sup>-4</sup></td><td>2.29</td></tr><tr><td>M<sub>3</sub></td><td>6.765 × 10<sup>-3</sup></td><td>1.28</td><td>1.982 × 10<sup>-4</sup></td><td>2.23</td></tr><tr><td>M<sub>4</sub></td><td>5.030 × 10<sup>-3</sup></td><td>1.23</td><td>1.111 × 10<sup>-4</sup></td><td>2.20</td></tr></tbody></table>

We now present a brief analysis of the influence of the parameters on the error, particularly in relation to the nonlinear terms of the equation. In set 2, we chose higher values for $\Gamma_{a,b}$ and $\Omega_{a,b}$ compared to $K_{a,b}$ and $r_{a,b}$. Conversely, in set 3, we assigned lower values to $\Gamma_{a,b}$ and $\Omega_{a,b}$ while increasing $K_{a,b}$ and $r_{a,b}$. The objective is to examine how the significance of the nonlinear terms impacts each error component. Figure 4 shows the evolution of each errors at a logarithmic scale, for set 2 and set 3.

<figure><img src="image_6.png" /><figcaption>Figure 4: Evolution of the logarithmic error for each component, for set 2 and set 3 for DGTD $P_1$ and DGTD $P_2$ methods.</figcaption></figure>
We observe that all errors tend to exhibit a convergence order of 1 for $\mathbb{P}_1$ and 2 for $\mathbb{P}_2$. Additionally, we notice that $e_{pol}^h$ and $e_{elec}^h$ are significantly higher for set 3 compared to set 2. This highlights the impact of the nonlinear terms on the approximation of the polarization vectors, and more notably on the electronic density. However, when examining the behavior of $e_{tot}^h$, we do not observe a significant difference between sets 2 and 3. This is because the approximation error of the electromagnetic fields remains larger than that of the ODEs components, making it the dominant contributor to the total error.

We present in Table 5 the values used for the evolution of the error in Figure 4 and other associated order of convergence.

<table><caption>Table 5: Evolution of errors and convergence orders for DGTD $\mathbb{P}_1$ and DGTD $\mathbb{P}_2$ methods for second order leap frog scheme for set 3 and 4.</caption><thead><tr><th>Mesh</th><th>Set</th><th>e<sup>h</sup><sub>E,H</sub></th><th>e<sup>h</sup><sub>pol</sub></th><th>e<sup>h</sup><sub>elec</sub></th><th>e<sup>h</sup><sub>tot</sub></th><th>Order</th></tr></thead><tbody><tr><td rowspan="6">P<sub>1</sub></td><td>M<sub>1</sub></td><td>2</td><td>2.41 × 10<sup>-2</sup></td><td>1.999 × 10<sup>-4</sup></td><td>4.991 × 10<sup>-6</sup></td><td>2.411 × 10<sup>-2</sup></td><td>-</td></tr><tr><td>M<sub>1</sub></td><td>3</td><td>2.395 × 10<sup>-2</sup></td><td>2.149 × 10<sup>-3</sup></td><td>9.673 × 10<sup>-3</sup></td><td>2.592 × 10<sup>-2</sup></td><td>-</td></tr><tr><td>M<sub>2</sub></td><td>2</td><td>2.411 × 10<sup>-2</sup></td><td>8.912 × 10<sup>-5</sup></td><td>2.150 × 10<sup>-6</sup></td><td>9.862 × 10<sup>-3</sup></td><td>1.38</td></tr><tr><td>M<sub>2</sub></td><td>3</td><td>9.938 × 10<sup>-3</sup></td><td>1.051 × 10<sup>-3</sup></td><td>4.540 × 10<sup>-3</sup></td><td>1.010 × 10<sup>-2</sup></td><td>1.33</td></tr><tr><td>M<sub>3</sub></td><td>2</td><td>6.320 × 10<sup>-3</sup></td><td>5.848 × 10<sup>-5</sup></td><td>1.401 × 10<sup>-6</sup></td><td>6.321 × 10<sup>-3</sup></td><td>1.26</td></tr><tr><td>M<sub>3</sub></td><td>3</td><td>6.400 × 10<sup>-3</sup></td><td>6.989 × 10<sup>-4</sup></td><td>2.970 × 10<sup>-3</sup></td><td>7.090 × 10<sup>-3</sup></td><td>1.30</td></tr><tr><td rowspan="6">P<sub>2</sub></td><td>M<sub>4</sub></td><td>2</td><td>4.658 × 10<sup>-3</sup></td><td>4.370 × 10<sup>-5</sup></td><td>1.044 × 10<sup>-6</sup></td><td>4.658 × 10<sup>-3</sup></td><td>1.22</td></tr><tr><td>M<sub>4</sub></td><td>3</td><td>4.718 × 10<sup>-3</sup></td><td>5.253 × 10<sup>-4</sup></td><td>2.221 × 10<sup>-3</sup></td><td>5.241 × 10<sup>-3</sup></td><td>1.25</td></tr><tr><td>M<sub>1</sub></td><td>2</td><td>1.900 × 10<sup>-3</sup></td><td>1.767 × 10<sup>-5</sup></td><td>4.689 × 10<sup>-7</sup></td><td>1.897 × 10<sup>-3</sup></td><td>-</td></tr><tr><td>M<sub>1</sub></td><td>3</td><td>1.906 × 10<sup>-2</sup></td><td>32.275 × 10<sup>-4</sup></td><td>9.086 × 10<sup>-4</sup></td><td>2.123 × 10<sup>-3</sup></td><td>-</td></tr><tr><td>M<sub>2</sub></td><td>2</td><td>4.298 × 10<sup>-4</sup></td><td>4.136 × 10<sup>-6</sup></td><td>1.096 × 10<sup>-7</sup></td><td>4.300 × 10<sup>-4</sup></td><td>2.30</td></tr><tr><td>M<sub>2</sub></td><td>3</td><td>4.370 × 10<sup>-4</sup></td><td>5.448 × 10<sup>-5</sup></td><td>2.107 × 10<sup>-4</sup></td><td>4.882 × 10<sup>-4</sup></td><td>2.27</td></tr><tr><td>M<sub>3</sub></td><td>2</td><td>1.868 × 10<sup>-4</sup></td><td>1.809 × 10<sup>-6</sup></td><td>4.810 × 10<sup>-8</sup></td><td>1.869 × 10<sup>-4</sup></td><td>2.24</td></tr><tr><td>M<sub>3</sub></td><td>3</td><td>1.904 × 10<sup>-4</sup></td><td>2.396 × 10<sup>-4</sup></td><td>9.184 × 10<sup>-5</sup></td><td>2.123 × 10<sup>-4</sup></td><td>2.23</td></tr><tr><td>M<sub>4</sub></td><td>2</td><td>1.041 × 10<sup>-4</sup></td><td>9.956 × 10<sup>-7</sup></td><td>2.713 × 10<sup>-8</sup></td><td>1.041 × 10<sup>-4</sup></td><td>2.21</td></tr><tr><td>M<sub>4</sub></td><td>3</td><td>1.062 × 10<sup>-4</sup></td><td>3.413 × 10<sup>-5</sup></td><td>5.1209 × 10<sup>-5</sup></td><td>1.187 × 10<sup>-4</sup></td><td>2.20</td></tr></tbody></table>
# 5 Conclusion

This study demonstrated the effectiveness of a high-order Discontinuous Galerkin Time-Domain (DGTD) method for modeling wave propagation in gain media. The proposed approach, coupled with a second-order Leap-Frog scheme, ensured stability through energy estimates and accurately captured the nonlinear dynamics of the four-level atomic model. Numerical validation in a 3D cubic cavity confirmed the method's accuracy and convergence, highlighting its potential for advanced laser physics simulations. However, further improvements can be made, particularly in the matrix inversion method used for solving the implicit system, to enhance computational efficiency and scalability.

# A Appendix

## Construction of the artificial analytic solution for a cubic PEC problem

### Spatial functions for artificial solutions

$$ 
\left\{ 
\begin{aligned}
f(x, y, z) &= \cos(x_m x) \sin(x_m y) \sin(x_m z), \\
g(x, y, z) &= \sin(x_m x) \sin(x_m y) \cos(x_m z), \\
h(x, y, z) &= \sin(x_m x) \cos(x_m y) \cos(x_m z), \\
i(x, y, z) &= \cos(x_m x) \sin(x_m y) \cos(x_m z), \\
j(x, y, z) &= \cos(x_m x) \cos(x_m y) \sin(x_m z).
\end{aligned}
\right. 
$$

### Parameters for artificial polarization vectors

$$ Z_a = \frac{K_a}{(\Gamma_a \omega)^2 + (\Omega_a^2 - \omega^2)^2} \quad \text{and} \quad Z_b = \frac{K_b}{(\Gamma_b \omega)^2 + (\Omega_b^2 - \omega^2)^2}. $$
Electronic densities constants

The constants $\alpha_i$, $\beta_i$, $\gamma_i$ and $\delta_i$ for $i \in \{0, 1, 2, 3\}$ are defined as follows:

$$
\left\{
\begin{aligned}
\alpha_3 &= \frac{r_b}{1 + 4\omega^2 l_{32}^2} (B_b l_{32}^2 \omega - A_b l_{32} (1 + 2l_{32}^2 \omega^2)) \\
\beta_3 &= \frac{r_b}{1 + 4\omega^2 l_{32}^2} (B_b l_{32} + 2A_b l_{32}^2 \omega) \\
\gamma_3 &= \frac{r_b}{1 + 4\omega^2 l_{32}^2} (A_b l_{32} - 2B_b l_{32}^2 \omega) \\
\delta_3 &= \beta_3 l_{32} \omega
\end{aligned}
\right.
$$

$$
\left\{
\begin{aligned}
\alpha_2 &= -\frac{\alpha_3 l_{21}}{l_{32} - l_{21}} - \gamma_2 - \delta_2 \\
\beta_2 &= \frac{1}{1 + 4\omega^2 l_{21}^2} \left( \left( \frac{\beta_3}{l_{32}} + r_a B_a \right) l_{21} + 2 \left( \frac{\gamma_3}{l_{32}} + r_a A_a \right) l_{21}^2 \omega \right) \\
\gamma_2 &= \frac{1}{1 + 4\omega^2 l_{21}^2} \left( \left( \frac{\gamma_3}{l_{32}} + r_a A_a \right) l_{21} - 2 \left( \frac{\beta_3}{l_{32}} + r_a B_a \right) l_{21}^2 \omega \right) \\
\delta_2 &= \beta_2 l_{21} \omega + \frac{\delta_3 l_{21}}{l_{32}}
\end{aligned}
\right.
$$

$$
\left\{
\begin{aligned}
\alpha_1 &= -\frac{\alpha_2 l_{10}}{l_{21} - l_{10}} - \frac{\alpha_3 l_{32} l_{10}}{(l_{32} - l_{21})(l_{32} - l_{10})} - \gamma_1 - \delta_1 \\
\beta_1 &= \frac{1}{1 + 4\omega^2 l_{10}^2} \left( \left( \frac{\beta_2}{l_{21}} - r_a B_a \right) l_{10} + 2 \left( \frac{\gamma_2}{l_{21}} - r_a A_a \right) l_{10}^2 \omega \right) \\
\gamma_1 &= \frac{1}{1 + 4\omega^2 l_{10}^2} \left( \left( \frac{\gamma_2}{l_{21}} - r_a A_a \right) l_{10} - 2 \left( \frac{\beta_2}{l_{21}} - r_a B_a \right) l_{10}^2 \omega \right) \\
\delta_1 &= \beta_1 l_{10} \omega + \frac{\delta_2 l_{10}}{l_{21}}
\end{aligned}
\right.
$$

$$
\left\{
\begin{aligned}
\alpha_0 &= \beta_0 \omega + \frac{\delta_1}{l_{10}} = 0 \\
\beta_0 &= \frac{1}{2\omega} \left( \frac{\gamma_1}{l_{10}} - r_b A_b \right) \\
\gamma_0 &= -\frac{1}{2\omega} \left( \frac{\beta_1}{l_{10}} - r_b B_b \right) \\
\delta_0 &= -\gamma_0 + \alpha_1 + \frac{\alpha_2 l_{21}}{l_{21} - l_{10}} + \frac{\alpha_3 l_{32}^2}{(l_{32} - l_{21})(l_{32} - l_{10})}
\end{aligned}
\right.
$$
References

[1] Shih-Hui Chang, Allen Taflove, *Finite-difference time-domain model of lasing action in a four-level two-electron atomic system*, Optics Express **12** (2004), no. 16, 3827–3833.

[2] Bo Wu, Bingbing Sun, Hui Xue, Feng Xiao, Zhixiang Huang, Xianliang Wu, *Study on the gain media with four energy level model using two dimensional fdtd method*, Physica B: Condensed Matter **423** (2013), 21–25.

[3] Stéphane Lanteri, Claire Scheid, Jonathan Viquerat, *Analysis of a generalized dispersive model coupled to a dgtd method with application to nanophotonics*, SIAM Journal on Scientific Computing **39** (2017), no. 3, A831–A859.

[4] Qing Xia, Jeffrey W. Banks, William D. Henshaw, Alexander V. Kildishev, Gregor Kovačič, Ludmila J. Prokopeva, Donald W. Schwendeman, *High-order accurate schemes for maxwell's equations with nonlinear active media and material interfaces*, Journal of computational physics **456** (2022), 111051.

[5] Amit S. Nagra, Robert A. York, *Fdtd analysis of wave propagation in nonlinear absorbing and gain media*, IEEE Transactions on antennas and propagation **46** (1998), no. 3, 334–340.

[6] Claire Scheid, Stéphane Lanteri, *Convergence of a discontinuous galerkin scheme for the mixed time domain maxwell's equations in dispersive media*, IMA Journal of Numerical Analysis **33** (2013), no. 2, 432–459.

[7] Hassan Fahs, Stéphane Lanteri, *A high-order non-conforming discontinuous galerkin method for time-domain electromagnetics*, Journal of computational and applied mathematics **234** (2010), no. 4, 1088–1096.

[8] Loula Fezoui, Stéphane Lanteri, Stéphane Lohrengel, Serge Piperno, *Convergence and stability of a discontinuous galerkin time-domain method for the 3d heterogeneous maxwell equations on unstructured meshes*, ESAIM: Mathematical Modelling and Numerical Analysis **39** (2005), no. 6, 1149–1176.

[9] Anan Fang, Thomas Koschny, Costas M Soukoulis, *Lasing in metamaterials nanostrucutres*, Journal of optics **12** (2010), no. 2, 024013.

[10] ——, *Self-consistent calculations of loss-compensated fishnet metamaterials*, Physical Review B—Condensed Matter and Materials Physics **82** (2010), no. 12, 121102.

[11] Maohui Lyu, Vrushali A. Bokil, Yingda Cheng, Fengyan Li, *Energy stable nodal discontinuous galerkin methods for nonlinear maxwell's equations in multi-dimensions*, Journal of Scientific Computing **89** (2021), no. 2, 45.

[12] Coleman, William F, *Lasers-an introduction*, 1982.

[13] Vrushali A. Bokil, Yingda Cheng, Yan Jiang, Fengyan Li, *Energy stable discontinuous galerkin methods for maxwell's equations in nonlinear optical media*, Journal of Computational Physics **350** (2017), 450–452.

[14] Hui Xue, Zhixiang Huang, Xianliang Wu, *Study on the gain material with four energy level model using fdtd method*, Optics and photonics Journal **3** (2013), 15–20.

[15] Germán J. de Valcárcel, Eugenio Roldan, Franco Prato, *Semiclassical theory of amplification and lasing*, Revista mexicana de física **52** (2006), no. 2, 198–214.
[16] BR Mollow, *Power spectrum of light scattered by two-level systems*, Physical Review **188** (1969), no. 5.

[17] Shouyuan Shi, Brandon Redding, Tim Creazzo, Elton Marchena, Dennis W Prather, et al., *Quantum electrodynamic modeling of silicon-based active devices*, Advances in Optical Technologies **2008** (2008).

[18] Jonathan Viquerat, *Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous galerkin time-domain method*, Ph.D. thesis, INRIA d'Université Côte d'Azur, 2018.
Inria

RESEARCH CENTRE
Centre Inria d'Université Côte d'Azur

2004 route des Lucioles - BP 93
06902 Sophia Antipolis Cedex

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

ISSN 0249-6399