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

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

**Abstract** This paper is concerned with numerical modeling of the interaction between an electromagnetic field and a gain medium in the context of laser physics, with the goal of simulating the gain process that results in an increase in optical power. This phenomenon can be modelled by a four-level atomic differential system that couples Maxwell’s equations with a set of non-linear Ordinary Differential Equations (ODEs) to describe the electronic density evolution for each energy level. Most of the existing works dealing with this model consider the Finite Difference Time-Domain method (FDTD) as seen in [1]. In this article, we will present a novel numerical modeling leveraging a Discontinuous Galerkin Time-Domain method (DGTD) in 3D that we have formulate to solve this model and will propose an estimation for the continuous energy associated to the system. Based on the work done in [5], we used a second order Leap-Frog temporal scheme and made approximation for the nonlinear terms present in the ODEs. An energy estimate, for both continuous problem and the discrete scheme, inspired by the work done in [3] allows us to prove stability of the scheme. The method is validated in a 3D framework using a model problem with manufactured solution.

## 1.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 an electromagnetic radiation with a collection of atoms. This interaction has two effects. On one hand, the medium has an effect on the field, and 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. In

<footnote>Cédric Legrand · Stéphane Lanteri · Stéphane Descombes
Université Côte d'Azur, INRIA, CNRS, LJAD
e-mail: cedric.legrand@inria.fr</footnote>
fact, while the atoms lose energy through radiative or non-radiative mechanisms, the medium attenuates or amplifies, and phase-shifts the total field [1].

In this process, the electrons have to be pumped to different energy levels. The most common model of this process considers that the electrons reach four different levels of energy [2]. We denote by these four different levels of energy, $E_0$, $E_1$, $E_2$, $E_3$. First, electrons are pumped by an external mechanism to the third level $E_3$. They quickly relax to $E_2$, non radiatively. Then they transfer from $E_2$ to $E_1$ both radiatively and non radiatively when the density population at $E_2$ is higher than the one at $E_1$; this is the mechanism known as lasing. Finally they transfer quickly and non radiatively form $E_1$ to $E_0$ [2]. All the parameters presented in figure 1.1 will be explained in section 1.2.

<figure><img src="image_1.png" /><figcaption>Fig. 1.1: Four-level atomic system model of a gain medium.</figcaption></figure>

## 1.2 General model

In our problem we will consider two domains. The first one is a vacuum medium, surrounding the second one, the gain medium. We will study the interaction between the gain medium and the incident electromagnetic wave. Let us introduce $\Omega_1, \Omega_2 \subset \mathbb{R}^3$, which are bounded domains corresponding respectively to the gain medium and to the host medium. We denote by $\mathbf{E}$ and $\mathbf{H}$, respectively the electric field and the magnetic field vectors in $\mathbb{R}^3$. $\epsilon$ and $\mu$ are respectively the electric permittivity and the magnetic permeability of the material. They both depend on the position. We introduce $\mathbf{P}$ the electric polarization density of the gain material and $\mathbf{J}_p$ his derivative in the gain material. $\mathbf{J}_p$ will be considered equal to zero in $\Omega_2$. $\mathbf{E}$, $\mathbf{H}$ are functions of space and time and are solutions of the following time-dependent Maxwell's equation on the time domain $[0, T]$, with $T > 0$:

$$ \left\{ \begin{array}{l} \epsilon \frac{\partial \mathbf{E}}{\partial t} = \nabla \times \mathbf{H} - \mathbf{J}_p \quad \text{in } \Omega_1 \cup \Omega_2 \times [0, T], \\ \\ \mu \frac{\partial \mathbf{H}}{\partial t} = -\nabla \times \mathbf{E} \quad \text{in } \Omega_1 \cup \Omega_2 \times [0, T]. \end{array} \right. $$

At the boundary of the domain we impose perfect electric conditions as follows, $\mathbf{n} \times \mathbf{H} = J_s$, $\mathbf{n} \times \mathbf{E} = 0$, $\mathbf{n} \cdot \mathbf{H} = 0$, $\mathbf{n} \cdot \mathbf{E} = \rho_s / \epsilon$.

The atomic density at the $i^{th}$ energy level is denoted by $N_i$, and is a function of time and space. Introducing $\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$ the induced electric polarization density on the pumping transition between the ground state level ($E_0$) and the third level ($E_3$). The total electric polarization of the gain media is $\mathbf{P} = \mathbf{P}_a + \mathbf{P}_b$ and satisfies the equations [2] :

$$ \left\{ \begin{aligned} \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{aligned} \right. \quad (1) $$

The evolution of the electron 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$ 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 [2] we have to consider the following rate equations :

$$ \left\{ \begin{aligned} \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], \\ \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], \\ \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], \\ \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{aligned} \right. \quad (2) $$

The initial conditions for the electronic densities are taken as $N_i(\mathbf{x}, 0) = 0$, for $i \in \{1, 2, 3\}$ and $N_0(\mathbf{x}, 0) = N$. The relaxation time of an electron from $E_i$ to $E_j$ is given by $\tau_{ij}$. We denote respectively by $\omega_a$ and $\omega_b$ the lasing and the pumping frequencies, such as $\omega_a = (E_2 - E_1)/\hbar$ and $\omega_b = (E_3 - E_0)/\hbar$. As defined in [2] the constants $\gamma_a$ and $\gamma_b$ are the line width of the atomic transition $\omega_a$ and $\omega_b$, and $\kappa_a$ and $\kappa_b$ are the coupling strength of $\mathbf{P}_a$ and $\mathbf{P}_b$ to the electric field.

We then performed several transformations of (1) and (2) in order to facilitate the theoretical study. We use the derivatives of $\mathbf{P}_a$ and $\mathbf{P}_b$, denoted as $\mathbf{J}_a$ and $\mathbf{J}_b$ respectively, to transform this system into a first-order system. With these new notations, we have $\mathbf{J}_p = \mathbf{J}_a + \mathbf{J}_b$. We also homogenized these systems of equations following the method presented in [4]. The new parameters of the problem are summarized in table 1.1.
<table><caption>Table 1.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>

## 1.3 Energy estimate and stability

To study the stability of the system we are going to rewrite the system with operators, by adapting the method used in [3] to our problem. Let us consider the following spaces :

$$ \mathbf{L}^2(\Omega_i) = \mathbf{L}^2(\Omega_i)^3, \mathbf{H}(\mathrm{curl}, \Omega_2) = \{\mathbf{v} \in \mathbf{L}^2(\Omega_2), \text{ such as } \nabla \times \mathbf{v} \in \mathbf{L}^2(\Omega_2)\}, \\ \text{and } \mathcal{E} = \mathbf{L}^2(\Omega_2)^2 \times \mathbf{L}^2(\Omega_1)^4 \times \mathbf{L}^2(\Omega_1) $$

Let $\mathcal{U}$ be the vector of unknowns, $\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$.

Let us also define by $\mathcal{A}$, the linear operator containing the linear part of the equations, and by $\mathcal{F}$ the non linear application that contain the non linear part of the equations :

$$ \begin{aligned} \mathcal{A} : \mathcal{D}(\mathcal{A}) &= \mathbf{H}(\mathrm{curl}, \Omega_2)^2 \times \mathbf{L}^2(\Omega_1)^4 \times \mathbf{L}^2(\Omega_1) \subset \mathcal{E} \rightarrow \mathcal{E}, \\ \mathcal{F} : \mathcal{D}(\mathcal{F}) &= \mathbf{H}(\mathrm{curl}, \Omega_2) \cap \mathbf{L}^3(\Omega_2) \times \mathbf{H}(\mathrm{curl}, \Omega_2) \times \mathbf{L}^3(\Omega_1)^2 \cap \mathbf{L}^2(\Omega_1)^2 \times \mathbf{L}^3(\Omega_1)^4 \rightarrow \mathcal{E}. \end{aligned} $$

Using the fact that $\mathcal{D}(\mathcal{F}) \subset \mathcal{D}(\mathcal{A})$, we can write the equations in a system form such as the problem becomes, with $\mathcal{U}_0 \in \mathcal{D}(\mathcal{F})$ an initial condition :

Find $\mathcal{U} \in \mathcal{D}(\mathcal{F})$ such as,

$$ \left\{ \begin{aligned} \frac{\partial \mathcal{U}}{\partial t} &= \mathcal{A}\mathcal{U} + \mathcal{F}(\mathcal{U}), \\ \mathcal{U}(t=0) &= \mathcal{U}_0. \end{aligned} \right. $$

We define a scalar product on $\mathcal{E}$ denoted by $\langle \cdot, \cdot \rangle_{\mathcal{E}}$ and $|| \cdot ||_{\mathcal{E}}$ its associated norm, such that, $\forall (\mathcal{U}, \mathcal{V}) \in \mathcal{E}^2$ :

$$ \begin{aligned} \langle \mathcal{U}, \mathcal{V} \rangle_{\mathcal{E}} ={}& (\epsilon_r \mathbf{u}_1, \mathbf{v}_1)_{\mathbf{L}^2(\Omega_2)} + (\mu_r \mathbf{u}_2, \mathbf{v}_2)_{\mathbf{L}^2(\Omega_2)} \\ & + \frac{1}{K_a} (\mathbf{u}_3, \mathbf{v}_3)_{\mathbf{L}^2(\Omega_1)} + \frac{1}{K_b} (\mathbf{u}_4, \mathbf{v}_4)_{\mathbf{L}^2(\Omega_1)} + \frac{\Omega_a^2}{K_a} (\mathbf{u}_5, \mathbf{v}_5)_{\mathbf{L}^2(\Omega_1)} + \frac{\Omega_b^2}{K_b} (\mathbf{u}_6, \mathbf{v}_6)_{\mathbf{L}^2(\Omega_1)} \\ & + \frac{1}{r_b} (u_7, v_7)_{\mathbf{L}^2(\Omega_1)} + \frac{1}{r_a} (u_8, v_8)_{\mathbf{L}^2(\Omega_1)} + \frac{1}{r_a} (u_9, v_9)_{\mathbf{L}^2(\Omega_1)} + \frac{1}{r_b} (u_{10}, v_{10})_{\mathbf{L}^2(\Omega_1)}. \end{aligned} $$
With this scalar product we find, for $\mathcal{U} \in \mathcal{D}(\mathcal{F})$ a solution of the problem with PEC boundary conditions :

$$ egin{cases} 
\langle \mathcal{F}(\mathcal{U}), \mathcal{U} \rangle_{\varepsilon} = 0, \\ 
\langle \mathcal{A}\mathcal{U}, \mathcal{U} \rangle_{\varepsilon} = -((\mathbf{J}_a + \mathbf{J}_b), \mathbf{E})_{L^2(\Omega_1)} - \frac{\Gamma_a}{K_a} \|\mathbf{J}_a\|_{L^2(\Omega_1)} - \frac{\Gamma_b}{K_b} \|\mathbf{J}_b\|_{L^2(\Omega_1)} \\ 
- \frac{1}{r_b l_{32}} \|N_3\|_{L^2(\Omega_1)} - \frac{1}{r_a l_{21}} \|N_2\|_{L^2(\Omega_1)} - \frac{1}{r_a l_{10}} \|N_1\|_{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)}. 
\end{cases} $$

For $\mathcal{U}$ a solution of the problem, we define $\xi$ the energy associated to the system such as, for $t > 0$ :

$$ \xi(t) = \frac{1}{2} \|\mathcal{U}(t)\|_{\varepsilon}^2. $$

For $t > 0$, we find by differentiating this expression :

$$ egin{aligned} \frac{\partial \xi}{\partial t}(t) &= \langle \frac{\partial \mathcal{U}}{\partial t}(t), \mathcal{U}(t) \rangle_{\varepsilon} = \langle \mathcal{A}\mathcal{U}(t), \mathcal{U}(t) \rangle_{\varepsilon} \\ &= -((\mathbf{J}_a + \mathbf{J}_b), \mathbf{E})_{L^2(\Omega_1)} - \frac{\Gamma_a}{K_a} \|\mathbf{J}_a\|_{L^2(\Omega_1)} - \frac{\Gamma_b}{K_b} \|\mathbf{J}_b\|_{L^2(\Omega_1)} \\ &\quad - \frac{1}{r_b l_{32}} \|N_3\|_{L^2(\Omega_1)} - \frac{1}{r_a l_{21}} \|N_2\|_{L^2(\Omega_1)} - \frac{1}{r_a l_{10}} \|N_1\|_{L^2(\Omega_1)} \\ &\quad + \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)}. \end{aligned} $$

In order to evaluate the derivative of the energy over time, we also assume that $\epsilon_r \ge \epsilon_-$ > 0. After some computations, we obtain the following estimate for the variation of the energy :

$$ \frac{\partial \xi}{\partial t}(t) \le 2C\xi(t), \text{ with } C = \sqrt{\frac{(K_a + K_b)}{\epsilon^{-}}} + \frac{1}{l_{32}} \sqrt{\frac{r_b}{r_a}} + \frac{1}{l_{21}} + \frac{1}{l_{10}} \sqrt{\frac{r_a}{r_b}}. $$

Let $t \in [0, T]$, by integrating the previous inequality over $[0, t]$, we obtain :

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

Under the assumption that the solution of the problem exists, we have thus shown that the $L^2$ norm of the solution is bounded on $[0, T]$.
# 1.4 Discontinuous Galerkin method and time discretization

## 1.4.1 Semi-discrete problem in time

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

$$ U^{[n+1/2]} = \frac{U^{n+1/2} + U^{n+3/2}}{2} \quad \text{and} \quad V^{[n+1/2]} = \frac{V^{n+1/2} + V^{n+3/2}}{2}. $$

## 1.4.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 set up the matrix formulation for the fully discrete scheme.

Let $\Omega_h$ be a discretization of $\Omega$, relying on a quasi-uniform triangulation $\mathcal{T}_h = (T_i)_{\ll 1, N\rr}$, where $N \in \mathbb{N}^*$. We denote $a_{ij} = T_i \cap T_j$, an internal face of the discretization for $(i, j) \in \ll 1, N\rr$ such as $T_i$ and $T_j$ are adjacent. Let $n_{ij}$ be the unit normal vector of $a_{ij}$, oriented from $T_i$ to $T_j$. For $i \in \ll 1, N\rr$, we denote $V_i = \{j \in \ll 1, N\rr, T_i \cap T_j \neq \emptyset\}$. We introduce the discrete space $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 freedom per dimension.

For a vectorial semi-discrete variable by $U^n$ and a scalar semi-discrete variable $V^n$ we respectively denote by $U_h^n$ and $V_h^n$ the discrete in space and times variables. Their restriction on each cell $T_i$ are denoted respectively $U_i^n$ and $V_i^n$. We will denote their coordinates, for a cell $T_i$, $i \in \ll 1, d_i\rr$, in the polynomial basis as follows, respectively :

$$ \overline{V}_i^n = (V_{i1}, \dots, V_{id_i})^T \quad \text{and} \quad \overline{U}_i^n = (\overline{U}_i^{x^n}, \overline{U}_i^{y^n}, \overline{U}_i^{z^n})^T. $$
For Maxwell's equation, we use the same method as in [4]. We introduce the three matrices defined in [4], the mass matrix $\overline{M}_i$, the stiffness matrix $\overline{K}_i$, the flux matrix $\overline{S}_i$. With the previous notations, we obtain directly the following matrix formulation for the Maxwell's equations:

$$ (\text{ME}) \left\{ \begin{array}{l} \epsilon_r \overline{M}_i \frac{\overline{E}_i^{n+1} - \overline{E}_i^n}{\Delta t} = \overline{K}_i \times \overline{H}_i^{n+1/2} - \overline{M}_i (\overline{J}_{ai}^{[n+1/2]} + \overline{J}_{bi}^{[n+1/2]}) \\ - \sum_{l \in V_i} \overline{S}_{il} (\overline{H}_i^{n+1/2} \times \mathbf{n}_{il}), \\ \\ \mu_r \overline{M}_i \frac{\overline{H}_i^{n+3/2} - \overline{H}_i^{n+1/2}}{\Delta t} = -\overline{K}_i \times \overline{E}_i^{n+1/2} + \sum_{l \in V_i} \overline{S}_{il} (\overline{E}_i^{n+1} \times \mathbf{n}_{il}). \end{array} \right. $$

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 elements. For example, the integral on $T_i$ of the product of two scalar variables $V_i$ and $W_i$ and a basis function $\phi_{ik}$, $k \in [1, d_i]$ is:

$$ \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 will make an approximation. We have decided to simply project the nonlinear term 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] \right) \Big|_{x_r} \phi_{ir} \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 it will then be a matter of calculating the integral, which results in the presence of the mass matrix in the matrix formulation, which ultimately simplifies along 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:

$$ (\text{PE}) \left\{ \begin{array}{l} \frac{\overline{\mathbf{P}}_{ai}^{n+3/2} - \overline{\mathbf{P}}_{ai}^{n+1/2}}{\overline{\mathbf{P}}_{bi}^{n+3/2} - \overline{\mathbf{P}}_{bi}^{n+1/2}} = \overline{\mathbf{J}}_{ai}^{[n+1/2]}, \\ \\ \frac{\overline{\mathbf{P}}_{bi}^{n+3/2} - \overline{\mathbf{P}}_{bi}^{n+1/2}}{\overline{\mathbf{J}}_{ai}^{n+3/2} - \overline{\mathbf{J}}_{ai}^{n+1/2}} = \overline{\mathbf{J}}_{bi}^{[n+1/2]}, \\ \\ \frac{\overline{\mathbf{J}}_{ai}^{n+3/2} - \overline{\mathbf{J}}_{ai}^{n+1/2}}{\overline{\mathbf{J}}_{bi}^{n+3/2} - \overline{\mathbf{J}}_{bi}^{n+1/2}} = -\Gamma_a \overline{\mathbf{J}}_{ai}^{[n+1/2]} - \Omega_a^2 \overline{\mathbf{P}}_{ai}^{[n+1/2]} - K_a \overline{\mathbf{E}}_{21i}^{n+1}, \\ \\ \frac{\overline{\mathbf{J}}_{bi}^{n+3/2} - \overline{\mathbf{J}}_{bi}^{n+1/2}}{\Delta t} = -\Gamma_b \overline{\mathbf{J}}_{bi}^{[n+1/2]} - \Omega_b^2 \overline{\mathbf{P}}_{ai}^{[n+1/2]} - K_b \overline{\mathbf{E}}_{30i}^{n+1}, \end{array} \right. $$
with $\forall v \in \{x, y, z\}$, $\overline{EN}_{21i}^{v,n+1} = ((N_{2ij}^{[n+1/2]} - N_{1ij}^{[n+1/2]}) E_{ij}^{v,n+1})_{1 \le j \le d_i}$
and $\overline{EN}_{30i}^{v,n+1} = ((N_{3ij}^{[n+1/2]} - N_{0ij}^{[n+1/2]}) E_{ij}^{v,n+1})_{1 \le j \le d_i}$.

$$ (\text{RE}) \quad \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}_{bi}^{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}_{ai}^{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}_{ai}^{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}_{bi}^{v,n+1}, \end{aligned} \right. $$

with $E_{a_i}^{v,n} = (E_{ij}^{v,n+1} J_{a_{ij}}^{v,[n+1/2]})_{1 \le j \le d_i}$ and $E_{b_i}^{v,n} = (E_{ij}^{v,n+1} J_{b_{ij}}^{v,[n+1/2]})_{1 \le j \le d_i}$.

### 1.4.3 Fully discrete energy

In order to study the stability of the scheme we define the fully discrete energy :

$$ \begin{aligned} \mathcal{E}_{\Omega_h}^n ={}& \frac{1}{2} \int_{\Omega_h} \left( \epsilon_r |\mathbf{E}_h^n|^2 + \mu_r \mathbf{H}_h^{n+1/2} \cdot \mathbf{H}^{n-1/2} \right) \\ & + \frac{1}{2} \int_{\Omega_h} \left( \frac{1}{K_a} |\mathbf{J}_a^{n+1/2}|^2 + \frac{1}{K_b} |\mathbf{J}_b^{n+1/2}|^2 + \frac{\Omega_a^2}{K_a} |\mathbf{P}_a^{n+1/2}|^2 + \frac{\Omega_b^2}{K_b} |\mathbf{P}_b^{n+1/2}|^2 \right) \\ & + \frac{1}{2} \int_{\Omega_h} \left( \frac{1}{r_b} |N_{3h}^{n+1/2}|^2 + \frac{1}{r_a} |N_{2h}^{n+1/2}|^2 + \frac{1}{r_a} |N_{1h}^{n+1/2}|^2 + \frac{1}{r_b} |N_{0h}^{n+1/2}|^2 \right). \end{aligned} $$

To complete the definition of the fully discrete energy we discretize the boundary conditions as in [3], for $a_{il} \subset \partial\Omega_h$ we use, $E_i^n|_{a_{il}} = -E_l^n|_{a_{il}}$, $H_i^{n+1/2}|_{a_{il}} = H_l^{n+1/2}|_{a_{il}}$.

Then, by introducing $\rho = (\frac{1}{r_a l_{32}} + \frac{1}{r_b l_{10}})(r_b + r_a) + \frac{2}{l_{21}}$, we can show that under the assumption $\Delta t < \frac{1}{\sqrt{2\frac{K_a + K_b}{\epsilon_r} + \rho}}$,

the fully discrete energy satisfies the following bound :

$$ \mathcal{E}_{\Omega_h}^n \le \left( 1 + \left( 4\sqrt{2\frac{K_a + K_b}{\epsilon_r}} + 2\rho \right) \Delta t \right)^n \mathcal{E}_{\Omega_h}^0. $$
Thus the fully discrete energy is bounded for a finite time.

## 1.5 Results

To validate our solvers we built an artificial case with manufactured solutions. 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 condition. The simulations were conducted using elements ranging from $P_1$ to $P_2$. The set of parameters for the gain medium model is presented in table 1.2.

<table><caption>**Table 1.2:** Set of parameters for the following test.</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><th>Values</th><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></tbody></table>

We present in figure 1.2 the variation of the total error with respect to the mesh sizes $N \in \{11, 21, 31, 41\}$ at the logarithmic scale for $P_1$ and $P_2$ Lagrange elements. In figure 1.3 we can observe the evolution of the energy for both elements. As excepted, we see in table 1.3 an order of convergence of 1 for $P_1$ and 2 for $P_2$.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.2:** Logarithmic errors for DGTD $P_1$ and DGTD $P_2$ methods.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 1.3:** Evolution of total energy for $T = 1$ and $N = 41$ for DGTD $P_1$ and DGTD $P_2$ methods.</figcaption></figure>
<table><caption>**Table 1.3:** Evolution of errors and convergence orders for DGTD P<sub>1</sub> and DGTD P<sub>2</sub> methods.</caption><thead><tr><th rowspan="2">N</th><th colspan="2">P<sub>1</sub></th><th colspan="2">P<sub>2</sub></th></tr><tr><th>Error</th><th>Order</th><th>Error</th><th>Order</th></tr></thead><tbody><tr><td>11</td><td>2.522 × 10<sup>-2</sup></td><td>-</td><td>1.993 × 10<sup>-3</sup></td><td>-</td></tr><tr><td>21</td><td>1.048 × 10<sup>-2</sup></td><td>1.36</td><td>4.542 × 10<sup>-4</sup></td><td>2.29</td></tr><tr><td>31</td><td>6.765 × 10<sup>-3</sup></td><td>1.28</td><td>1.984 × 10<sup>-4</sup></td><td>1.28</td></tr><tr><td>41</td><td>5.030 × 10<sup>-3</sup></td><td>1.23</td><td>1.110 × 10<sup>-4</sup></td><td>2.20</td></tr></tbody></table>

## 1.6 Conclusion

The work carried out up to this point and presented in this article includes a theoretical study and the development of the DGTD method of a 4-level model of a gain medium. The next objective is to adapt this scheme to an upwind flux DG scheme in space coupled to a Low Storage Runge-Kutta scheme for the time integration, with the goal of achieving higher numerical convergence order.

## References

1. Amit S. Nagra, A. York. in *FDTD Analysis of Wave Propagation in Nonlinear Absorbing and gain Media*, IEEE Transactions on antennas and propagation (1998),

2. A. Fang, T. Koschny, C. M. Soukoulis. in *Lasing in metamaterials nanostrucutres*, IOP Science (2010),

3. L. Fezoui, S. Lanteri, S. Lohrengel, S. Piperno, *Convergence and stability of a discontinous Galerkin Time-Domain method for the 3D heterogeneous Maxwell equations on unstructured meshes*, ESAIM (2005),

4. J. Viquerat, *Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method*, INRIA d'Université Côte d'Azur (2018),

5. M. L. Vrushali, A. Bokil, Y. Cheng, F. Li, *Energy stable discontinuous Galerkin methods for Maxwell's equations in nonlinear optical media*, Rensselaer Polytechnic Institute (2017).