# Numerical modeling of wave propagation in passive and time-modulated media

Roman Gelly, Stéphane Lanteri, Mahmoud Elsawy

► To cite this version:

Roman Gelly, Stéphane Lanteri, Mahmoud Elsawy. Numerical modeling of wave propagation in passive and time-modulated media. RR-9569, Centre Inria d'Université Côte d'Azur. 2025, pp.49. hal-04873417

HAL Id: hal-04873417

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

Submitted on 8 Jan 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
Inria

Numerical modeling of
wave propagation in
passive and
time-modulated media

Roman Gelly, Stéphane Lanteri, Mahmoud Elsawy

RESEARCH
REPORT

Nº 9569

December 2024

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

# Numerical modeling of wave propagation in passive and time-modulated media

Roman Gelly*, Stéphane Lanteri *, Mahmoud Elsawy *

Project-Team Atlantis

Research Report n° 9569 — December 2024 — 49 pages

**Abstract:** This report summarizes the work achieved during my 12-month apprenticeship at the Inria research center at Université Côte d'Azur, as part of the Atlantis† project-team. The research activity of the team focuses on the numerical modeling of problems involving wave-matter interaction at the nanoscale (nanophotonics). The team mainly contributes to the development of computational methods to solve Maxwell's equations in heterogeneous media, implemented in the DIOGENeS‡ software suite.

This report focuses on the numerical study of *optical metasurfaces*, which are ultra-thin artificial materials composed of nanometric light-scattering elements, engineered to achieve a given optical function. Metasurfaces have already been proven to be useful to fields such as imaging, beam forming, holography and polarimetry. The first part of this report contains the numerical study of a passive phase-gradient metasurface through inverse design and simulation with a high-order Discontinuous Galerkin Time-Domain (DGTD) solver. The second part is motivated by the study of another type of metasurfaces called *active* or *time-modulated* metasurfaces. We first extend the DGTD method to propagation in time-varying media, then we validate its implementation through the comparison of numerical results with artificial and physical analytical solutions.

**Key-words:** nanophotonics, Discontinuous Galerkin, metasurfaces, active metasurfaces

* Project-Team Atlantis, Centre Inria d'Université Côte d'Azur

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

2004 route des Lucioles - BP 93
06902 Sophia Antipolis Cedex
# Modélisation numérique de la propagation d'ondes dans des milieux passifs et modulés temporellement

**Résumé** : Ce rapport contient le travail effectué durant mon année d'apprentissage au Centre Inria d'Université Côte d'Azur, au sein de l'équipe-projet Atlantis<fnref n="4" />. L'activité de recherche de l'équipe est centrée sur la modélisation numérique de problèmes d'interaction onde-matière à l'échelle nanométrique (nanophotonique). L'équipe contribue principalement au développement de méthodes numériques appliquées à la résolution des équations de Maxwell dans des milieux hétérogènes, implémentées dans la suite DIOGENES<fnref n="5" />.

Ce rapport est centré sur l'étude numérique de métasurfaces optiques, qui sont des matériaux artificiels ultra-fins constitués d'éléments diffractants nanométriques, conçus pour accomplir une fonction optique donnée. Les métasurfaces ont déjà prouvé leur utilité dans des domaines tels que l'imagerie, la formation de faisceaux, l'holographie et la polarimétrie. La première partie de ce rapport contient l'étude numérique d'une métasurface passive à gradient de phase par *design* inverse et des simulations avec un solveur temporel Galerkin Discontinu (DGTD) d'ordre élevé. La deuxième partie est motivée par l'étude d'un autre type de métasurfaces appelées métasurfaces *actives* ou *modulées temporellement*. Nous étendons d'abord la méthode DGTD à la propagation dans des milieux dépendants du temps, puis nous validons son implémentation par la comparaison de la solution numérique avec des solutions analytiques artificielles et physiques.

**Mots-clés** : nanophotonique, Galerkin Discontinu, metasurfaces, metasurfaces actives

<footnote n="4">⁴www-sop.inria.fr/atlantis/</footnote>

<footnote n="5">⁵diogenes.inria.fr</footnote>
CENTRE INRIA D'UNIVERSITÉ CÔTE D'AZUR

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

UNIVERSITÉ
CÔTE D'AZUR

Înria

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

Master 2 Apprenticeship:

Project-team ATLANTIS

sept. 2023 - sept. 2024

**Part I:** Inverse design of metasurfaces through numerical modeling

**Part II:** Numerical modeling of wave propagation in time-varying media

Roman GELLY

Supervised by:

Stéphane LANTERI
& Mahmoud ELSAWY
A Derivation of analytical solutions for a stationary wave in a resonant cavity filled with a time-varying medium

41

A.1 Derivation for a general permittivity modulation profile $\epsilon_r(t)$

41

A.2 Linear permittivity modulation

42

A.3 Quadratic permittivity modulation

45
# Part I

# Inverse design of metasurfaces through numerical modeling

## 1 Study of a diffraction grating case

The goal of this first section is to study the diffraction of an incident polychromatic plane wave by a periodic grating deposited on a insulator substrate in order to investigate the influence of numerical parameters on the accuracy of the solution and the computational cost of the simulation. The simulations are conducted using the DGTD method (Discontinuous Galerkin Time Domain) implemented in diogenes-dgtd of the DIOGENES suite.

### 1.1 Design of a numerical model for an infinite surface grating

The first step is to design the numerical model of the grating which we will use to run the DGTD solver. This includes defining the dimensions of the different components, the physical properties of the materials (ε, μ), as well as additional non-physical layers and surfaces to ensure the numerical accuracy of the solution and to impose the desired boundary conditions on the borders of the domain.

In our case, we seek to model an optical diffraction surface grating, composed of a repeating pattern of dielectric segments, deposited on an insulating substrate. To simplify the problem, we assume that the grating is doubly periodic along the lateral directions, which implies that the segments will be of infinite width with an infinitely repeating pattern.

In the following, we will use the structure dimensions and the materials of the diffraction grating provided as example with the diogenes-dgtd solver. Since the provided mesh is hexahedral, we can also rebuild the model using the diogenes-gfactory Python library, allowing for a more versatile numerical model. The structure of the model and the different materials are illustrated in Fig. 1.

The periodic pattern is composed of six 1 µm high gallium nitride (GaN) segments, deposited on a 300 nm thick titanium dioxide (TiO₂) substrate and surrounded by air above the grating. We also need to define several additional layers and surfaces to properly set up the numerical model.

First, we need to define the boundary conditions for the 6 domain borders. As mentioned before, we want to impose double periodicity along the lateral axes, hence we impose Periodic Boundary Condition (PBC) to the four lateral borders of the domain. For the remaining top and bottom borders, since we want to simulate a wave propagating through an infinite domain so as to prevent spurious reflections in the domain, we impose an absorbing boundary condition (ABC) on these surfaces, in the form of the Silver-$M\overline{A}^4$ller condition (SMC). This ABC allows waves to leave the domain, while preventing waves to enter. While SMC is very efficient in absorbing waves normally incident to the boundary, its performance decreases when the angle of incidence increases, creating spurious reflections. However, in the general case, we expect the scatterer to diffract the incident wave in multiple directions. To overcome this problem, we introduce an additional region around the physical domain called Perfectly Matched Layer (PML)[1], where the waves are artificially damped by several order of magnitude, by modifying Maxwell's equations in this region.

Since the plane wave will be injected from below, we need to define an injecting surface below
<figure><img src="image_5.png" /><figcaption>**Figure 1:** Diagram of the diffraction grating model and two types of meshing: (b) hexahedral and (c) tetrahedral.</figcaption></figure>

the substrate which will also be used to monitor the reflection of the wave by the grating. This surface is called a TF/SF interface and allows to separate the domain into a Total Field (TF) region, where the scatterer is located and a Scattered Field (SF) region where the observables will be computed. This allows to compute the observables on the scattered field without taking into account the contribution from the incident field.

We also need to define another surface between the grating and the top PML layer that will be used to capture the part of the wave transmitted upwards.

## 1.2 Simulation with the diogenes-dgtd solver

### 1.2.1 Evolution of energy and EM field

As a first simulation run, we start by running the DIOGENES DGTD solver using the numerical model and the default settings provided in DIOGENES. In this case, the mesh is hexahedral (Fig. 1b), which is well suited to the rectangular geometry of the model, allowing for a reduced number of elements and of similar sizes.

Now, it is necessary to choose a criterion to stop the simulation. The most straightforward criterion would be to impose a maximum physical time, but it our case, it would be difficult to predict the time taken by the wave to leave the domain. We thus use a criterion based on the energy decrease in the domain, meaning that the simulation is stopped when the energy decreases below a certain percentage of the maximum recorded energy. In our case, we first fix the energy criterion to 0.001 % of the maximum energy in order to fully capture the propagation
of the wave through the domain.

We must also choose the polynomial order for the spatial scheme (Discontinuous Galerkin). Since the example provided is discretized by an hexahedral mesh, we must configure the simulation for quadratic elements and we choose a polynomial order of 2 for the moment.

In Fig. 2a, we plot the evolution of the total energy contained in the domain as a function of physical time. For better visibility, we limit the time axis to period before reaching an energy lower than 1% of the maximum peak. We observe that the bulk of energy rapidly increases and

<figure><img src="image_6.png" /><figcaption>(a) Evolution of the total energy computed in the domain.
(b) Evolution of the $E_x$ field measured below (Probe 1) and above (Probe 2) the grating.</figcaption></figure>

**Figure 2:** Time evolution of the total energy and the $E_x$ field, truncated to the time period where the energy is above 1 % of the maximum peak.

decreases until reaching 1 % of the maximum at $t \approx 6 \times 10^{-14}$ s but afterwards, the energy continues to decrease slowly until reaching the 0.001 % criterion at $t \approx 9 \times 10^{-13}$ s.

Before computing the observables (mainly the reflection and transmission coefficients), we first need to make sure that the wave had the time to properly travel through the domain, interacting with the grating, then leaving the domain completely. This is what we call the temporal convergence of the solution. In order to monitor the temporal convergence of the solution, we monitor the time evolution of the electromagnetic field at certain key positions in the domain. In our case, we place a first probe point within the substrate (20 nm below the grating) and a second one in the air above (20 nm above the grating). We show in Fig. 2b the $E_x$ field measured at the two probe points. We can observe the incident pulse wave measured by Probe 1 and a similar waveform measured by Probe 2 after a $t \approx 10^{-14}$ s delay which corresponds to the wave transmitted by the grating. Later, we observe a second pulse measured by Probe 1, which corresponds to the reflection of the wave by the grating. We notice on the full $E_x$ time evolution, that the wave oscillations are slowly damped but remain visible until the end criterion is reached at $t \approx 9 \times 10^{-13}$ s.

This means that, while most of the wave energy escapes the domain at the beginning of the simulation, there remains residual energy for a long period of time, preventing the energy criterion to be fulfilled. This very long physical time translates into a very long computing time, which is a critical factor in the optimization process.

We will investigate the reason behind this very long physical time later when we will compute the diffraction efficiencies of the grating.
### 1.2.2 Reflection and transmission coefficients

Now that we investigated the temporal convergence of the numerical solution and ensured that most of the electromagnetic energy had escaped the domain at the end of the simulation, we will quantify how the grating diffracts the incident wave; in particular, which proportion of the light is transmitted towards the vertical upward direction and which proportion is reflected back in the opposite direction.

In order to quantify the power density carried by the electromagnetic wave as a function of the frequency, we use the time-averaged Poynting vector:

$$ \pi(\mathbf{x}, \omega) = \frac{1}{2} \mathcal{R} \left( \hat{\mathbf{E}}(\mathbf{x}, \omega) \times \hat{\mathbf{H}}^{*}(\mathbf{x}, \omega) \right) \quad (1) $$

where $\hat{\mathbf{E}}$ and $\hat{\mathbf{H}}$ are the Fourier transform of the electric and magnetic fields.

In the case of a periodic structure, we can define the reflection and transmission coefficients. The first is defined as the proportion of the incident field power traveling back through the incident TF/SF surface $S_i$:

$$ R(\omega) = \frac{\int_{S_i} \pi_{\text{sca}} \cdot \mathbf{n}}{\int_{S_i} \pi_{\text{inc}} \cdot \mathbf{n}} \quad (2) $$

The transmission coefficient is defined as the proportion of the incident field power being transmitted through the upper TF/SF surface $S_o$:

$$ T(\omega) = \frac{\int_{S_o} \pi_{\text{tot}} \cdot \mathbf{n}}{\int_{S_i} \pi_{\text{inc}} \cdot \mathbf{n}} \quad (3) $$

In the absence of dispersive materials, there is no absorption by the scatterer and the coefficients should verify $R(\omega) + T(\omega) = 1$. This relation will allow us to check the physical validity of the numerical solution and detect anomalies due to numerical errors.

In Fig. 3b, we plot the reflection $R$ and transmission $T$ coefficients as a function of the wavelength computed from the previous simulation (of polynomial order $Q_2$). We observe that for long wavelengths (> 700 nm), the sum of the two coefficients reaches 1, except for several specific wavelengths, where we observe drops in the curve, which coincide with drops in the transmission and peaks in the reflection. However for shorter wavelengths (400 < $\lambda$ < 700 nm), we observe a drop in the $R + T$ curve, reaching $R + T \approx 0.8$ at $\lambda = 400$ nm. We can attribute this effect to the insufficient spatial convergence of the DGTD solver for the chosen numerical parameters, causing the energy to dissipate artificially in the small wavelengths regime and thus artificially dampening the reflection and transmission coefficients.

The peaks and drops in the coefficients can be attributed to a resonance phenomenon happening for specific wavelengths and causing the wave to be reflected back instead of transmitted.

To achieve better numerical convergence, we can either refine the mesh of the numerical model, or increase the polynomial order in the numerical solver.

Without needing to remesh the model, we can first try to increase the polynomial order and see if we can obtain a more physical solution.

### 1.2.3 Influence of the polynomial order on the numerical solution

In order to compare the numerical solutions obtained with different polynomial orders, we conduct the same simulations in $Q_1$, $Q_2$, $Q_3$ and $Q_4$. Then, we compute the reflection and transmission coefficients in each case, represented in Fig. 3.

While running the simulation with a polynomial order of $Q_3$, we notice that the energy first
<figure><img src="image_7.png" /><figcaption>**Figure 3:** Wavelength dependence of the reflection and transmission coefficients for the infinite diffraction grating.</figcaption></figure>

reaches a minimum of ≈ 0.008% of the maximum peak, then diverges, never fulfilling the energy stopping criterion. We can observe the divergence of the numerical solution on the evolution of the $E_x$ field, represented in Fig. 4a. We observe the same phenomenon with a Q4 polynomial order (see Fig. 4b).

This numerical effect seems to be linked to round-off errors, taking place after iterating over a too large number of time steps (~ 71000 steps for Q3 and ~ 125000 steps for Q4). As a simple fix, we increased the energy stopping criterion to 0.01% and repeated the simulation. This allowed to stop the simulation process before reaching numerical divergence, however the CPU-time remained very large, limiting our capacity to conduct inverse design through an optimization method, which necessitate a large number of solver calls in order to converge. Using a stopping energy criterion of 0.01%, the, we obtain the total CPU-times reported in Tab. 1. The simulations were conducted on the Nef<fnref n="1" /> cluster running on 18 cores.

<footnote n="1">¹wiki.inria.fr/ClustersSophia/Clusters_Home</footnote>
<figure><img src="image_8.png" /><figcaption>Figure 4: Time evolution of the $E_z$ field for a $Q_3$ and $Q_4$ polynomial order.</figcaption></figure>

<table><caption>Table 1: Simulation CPU-times for the infinite grating with polynomial orders $Q_1$ to $Q_4$. The simulations were run on 18 cores of the Nef.</caption><thead><tr><th>Order</th><th>Q<sub>1</sub></th><th>Q<sub>2</sub></th><th>Q<sub>3</sub></th><th>Q<sub>4</sub></th></tr></thead><tbody><tr><td>Runtime (s)</td><td>25</td><td>990</td><td>13096</td><td>49544</td></tr></tbody></table>

## 1.3 Conclusion

This particular case of an infinite nano-grating in not well suited for inverse designing through optimization due to the appearance of physical resonances within the structure, increasing drastically the physical time of the simulation and thus its computational time. The choice of using hexahedral elements for the numerical model is well-suited to the geometry of the structure, but the CPU-time scaling with the polynomial order leads to very time-consuming simulations in order to obtain a physically accurate solution.

In the following section, we will study a simpler metasurface design and use a tetrahedral mesh for the numerical model.
# 2 Design of phase-gradient metasurfaces with a global optimization method

## 2.1 Motivation

We will now use the DIOGENES DGTD numerical solver combined with a Bayesian global optimization method (EGO) in order to find optimal designs for a type of metasurface called *phase-gradient metasurface*. This type of metasurface is capable of deflecting an incident beam towards an arbitrary direction, or multiple directions, due to nanopillars inducing a gradual phase-shift in the wavefront, effectively changing the direction of propagation of the wave.
In this study, we will consider a metasurface composed of rectangular *nanopillars* of spatially variable dimensions, made of gallium nitride (GaN), and deposited on a $Al_2O_3$ substrate.
We will use the implementation of the EGO method in the *diogenes-optim* Python package in order to conduct mono-objective and multi-objective optimization for the inverse design of a *metadeflector* and a *beam-splitter* respectively.

## 2.2 Numerical model

We consider a doubly periodic structure in the horizontal directions $\hat{x}$ and $\hat{y}$ composed of an infinitely repeating unit cell represented in Fig. 5. The latter is composed of four rectangular GaN nanopillars deposited on a $Al_2O_3$ substrate, as proposed in [2]. We terminate the domain vertically with PMLs and Absorbing Boundary Conditions and laterally with Periodic Boundary Conditions as for the previous case. We also define two TF/SF layers to inject an incident plane wave from below and monitor the reflection and transmission of the wave.

<figure><img src="image_9.png" /><figcaption>**Figure 5:** Numerical model for the gradient-phase metasurface.</figcaption></figure>
## 2.3 Study of numerical convergence

Prior to exploring optimal design for our gradient-phase metasurface through optimization, it is necessary to discuss the choice of numerical parameters for our simulations with the DIOGENES DGTD solver.

In particular, we study in this section the influence of the *mesh size* and the *polynomial order* on the accuracy of the numerical solution. We also compare the simulation runtimes between in each case, which is also a critical criterion to consider for the optimization of the metasurface.

<figure><img src="image_10.png" /><figcaption>**Figure 6:** Study of numerical convergence for a coarse mesh (top) and a finer mesh (bottom)</figcaption></figure>

In Fig. 6, we compare the numerical accuracy of the simulations for two types of mesh (a coarse and fine mesh) for the polynomial orders ($P_1$ to $P_4$) through the wavelength dependence of two quantities: the sum of the total reflectance and transmittance and the diffraction efficiency in transmission of the mode (0, -1), which is the quantity of interest we will optimize later.

The numerical accuracy obtained with a fine mesh (50K cells) and a polynomial order of $P_2$ or more is comparable to the one obtained with the coarse mesh (2.7K cells) in $P_4$. However the latter is far less costly in terms of computing time than the former (see Tab. 2), taking more than 6 times longer (~6 minutes against ~39 minutes). We will then favor the use of coarse mesh with a high polynomial order for the optimization process. Later, after the optimization process, we will further refine the mesh to visualize more accurately the electromagnetic field distribution within the domain.

<table><caption>**Table 2:** Simulation runtimes in seconds for polynomial orders $P_1$ to $P_4$, with a coarse mesh a fine mesh. The simulations were run on a single node of the Nef cluster on 64 cores.</caption><thead><tr><th>Polynomial order</th><th>P<sub>1</sub></th><th>P<sub>2</sub></th><th>P<sub>3</sub></th><th>P<sub>4</sub></th></tr></thead><tbody><tr><td>Coarse mesh (2.7K cells)</td><td>12s</td><td>59s</td><td>161s</td><td>382s</td></tr><tr><td>Fine mesh (50K cells)</td><td>813s</td><td>2334s</td><td>5459s</td><td>10849s</td></tr></tbody></table>
One important aspect to consider to further minimize the simulation runtimes and thus the duration of the optimization process is the CPU-time speed-up obtained by the parallel implementation of the DIOGENES DGTD solver. In Fig. 7, we show the simulation CPU-time and parallel speed-up obtained on a coarse mesh (2.7K cells) in $\mathbb{P}_2$ and $\mathbb{P}_4$.

<figure><img src="image_11.png" /><figcaption>Figure 7: CPU-time and parallel speed-up obtained on a coarse mesh with polynomial order $\mathbb{P}_2$ and $\mathbb{P}_4$. The simulation were run on a single node of the Nef cluster on up to 64 cores.</figcaption></figure>

In both cases ($\mathbb{P}_2$ and $\mathbb{P}_4$), we observe that the CPU time significantly decreases with the number of cores for a low number of CPUs, but the parallel speed-up quickly becomes less efficient and reaches a plateau. We obtain a maximum parallel speed-up of ≈ 11.3 on 48 cores for the $\mathbb{P}_2$ case, then becomes less efficient on 64 cores. For the $\mathbb{P}_4$ case, we obtain a better speed-up, reaching ≈ 14 on 64 cores. These results show that the parallel speed-up is not very efficient in our case, which can be due to the very coarse mesh we chose for our simulations. Indeed, in this case, each mesh partition contains a low number of cells and the communication cost between these partitions can become larger than the computational cost of the DGTD solver when dealing with many partitions.

For the optimization part of our study, we will choose the most efficient solution (64 cores in $\mathbb{P}_4$), but we could chose a lower number of cores without affecting significantly the optimization duration.

## 2.4 Global optimization with the EGO method

### 2.4.1 Optimization parameters and bounds

For our study, we consider an unit cell of horizontal periods $\Gamma_x = 300$ nm in the $\hat{x}$ direction and $\Gamma_y = 1500$ nm in the $\hat{y}$ direction.

We define our optimization parameters as the height $h$ of the pillars (identical for all pillars), the dimensions along $x$ and $y$ ($d_x^i, d_y^i$) of the four pillars and the three distances $L_i$ between each pillar. This gives us in total 12 optimization parameters that we can store in the 12D design vector:

$$ \xi = [h, L_1, L_2, L_3, d_x^1, d_y^1, d_x^2, d_y^2, d_x^3, d_y^3, d_x^4, d_y^4] \quad (4) $$

The positions of the optimization parameters are represented in Fig. 8.

We now need to define the parameter space, i.e. the bounds on each individual parameter. This is to ensure that no parameter configuration can lead to a non-physical design (for ex.
<figure><img src="image_12.png" /><figcaption>**Figure 8:** Top view of the unit cell and optimization parameters</figcaption></figure>

superposing pillars, a total length of the array larger than the dimension of the unit cell) during the optimization process. This would lead to either a bug in the mesh generation with Gmsh or non-physical results. The most rigorous way to proceed would be to conduct constrained optimization where the designs are constrained by the dimensions of the unit cell. Here, we choose a simpler approach, where we have to make assumptions on the bounds of the parameters to prevent non-physical designs.

For the height of the pillars, we impose $600 \text{ nm} < h < 800 \text{ nm}$. For their $x$ and $y$ dimensions, we impose a minimum feature size of 90 nm due to experimental constraints and we impose a minimum distance of 50 nm between pillars and between a pillar and the border of the unit cell ($\delta_x, \delta_y$). This prevents very small cells to be generated in the meshing, which would decrease the timestep and thus increase the CPU time. To further prevent the generation of small cells near the border of the domain, we also center the pillar array on the substrate along the two horizontal directions.

We still need to define the upper bounds of the dimensions and distances between pillars. For the $x$ dimension of pillars, accounting for the minimum distance, the upper bound is $\Gamma_x - 2\delta_x = 200 \text{ nm}$. The upper bounds on the $y$ dimension of pillars and the distances is a trade-off between the possibility of having large pillars and the possibility of having distant pillars in the designs explored by the optimization method. We fix the upper bound on the $d_y^i$ to 260 nm which constrains the maximum upper bound on the $L_i$ to 120 nm since we must verify $\sum_{i=1}^4 d_y^i + \sum_{i=1}^3 L_i + 2\delta_y \le \Gamma_y$. The lower and upper bounds of parameters are summarized in Tab. 3.

<table><caption>**Table 3:** Lower and upper bounds on the optimization parameters for the metadeflector design.</caption><thead><tr><td>Param.</td><td>Lower bound</td><td>Upper bound</td></tr></thead><tbody><tr><td>$h$</td><td>600 nm</td><td>800 nm</td></tr><tr><td>$L_i$</td><td>50 nm</td><td>120 nm</td></tr><tr><td>$d_x^i$</td><td>90 nm</td><td>200 nm</td></tr><tr><td>$d_y^i$</td><td>90 nm</td><td>260 nm</td></tr></tbody></table>

### 2.4.2 The EGO algorithm

In order to find the optimal design parameters for the gradient-phase metasurface, we search for the global minimum of a user-defined objective function (or a set of objective functions for multi-objective optimization). Since each evaluation of the objective function requires to run a CPU-intensive DGTD simulation, the number of evaluations necessary to reach convergence is a crucial parameter. Furthermore, the complex shape of the objective function in the parameter
space with possibly multiple local minimums leads us to the use of a global optimization method. We use a Bayesian optimization method named Efficient Global Optimization (EGO), which relies on building a surrogate model (most often a Gaussian Process also called Kriging model) from a set of samples evaluated from the objective function. This allows to approximate the objective function using a probabilistic representation with a mean $\mu$ and uncertainty $\sigma$ for its value at each point of the parameter space.

The model is iteratively refined by adding a new sampling point determined by the maximization of a criterion: the Expected Improvement (EI) function. The latter denotes the expected value for the improvement gained by the addition of a new design $\mathbf{x}$.

Let $Y$ be the random variable following the Gaussian process $N(\mu(\mathbf{x}), \sigma^2(\mathbf{x}))$ (surrogate model) and $\mathbf{x}_i$ the current set of sample points, the Expected Improvement is expressed as:

$$ \mathbb{E}[I(\mathbf{x})] = \mathbb{E}[\max(\min_i(f(\mathbf{x}_i) - Y(\mathbf{x})), 0)] \quad (5) $$

The next sampling point is determined by:

$$ \mathbf{x}_{i+1} = \arg\max_i(\mathbb{E}[I(\mathbf{x})]) \quad (6) $$

The method is initialized by sampling the objective functions at a number of random points, generally using a space-filling sampling called Latin Hypercube Sampling (LHS) which covers the parameter space evenly.

### 2.4.3 Mono-objective optimization

The objective of this section is to optimize the parameters of the design so that an incident TM-polarized plane wave impinging from under the substrate is deflected by the nanopillar array towards a single direction, achieving the function of a *metadeflector*. To achieve this function, we seek to maximize the diffraction efficiency $\eta^{\text{TM}}(0, -1)$ of the metasurface for the $(0, -1)$ mode, corresponding to a deflection angle of $\theta_{-1} \approx 24^\circ$ at $\lambda = 600$ nm. This quantity describes how the power transmitted by the structure divides between the different diffraction orders, normalized between 0 and 1. Thus, we can define our objective function as $f(\xi) = 1 - \eta^{\text{TM}}(0, -1)(\lambda = 600 \text{ nm})$, which we will minimize using the EGO algorithm implemented in the module *diogenes-optim* of the DIOGENES suite.

We start the optimization process by defining a Design Of Experiment (DOE) of 80 points in Latin Hypercube Sampling (LHS) and fix the number of subsequent EGO iterations to 150, totalling 230 calls to the *diogenes-dgtd* solver. After almost 20 hours of optimization (mean duration of ~5 minutes per iteration), we obtain the convergence graph in Fig. 9. We observe that the EGO method reaches quickly very optimized designs in the first iterations after the DOE phase, reaching deflection efficiencies above 85% 20 EGO steps after the DOE phase. The best deflection efficiency is found at iteration 183, reaching $\eta^{\text{TM}}(0, -1) = 90\%$. After iteration 200, we observe a sudden increase in the dispersion of the objective function values. This is due to the behavior of the EGO implementation, which aims to explore other areas of the parameter space after reaching a plateau in the convergence.

We report the values of the optimized design parameters in Tab. 4. We notice that the pillars are almost equidistant (around $L_i \approx 100$ nm) and their dimensions along $x$ and $y$ decrease spatially withing the unit cell, except for the last pillar, for which its $y$ dimension is slightly larger by 0.2 nm than the previous pillar. For the first three pillars, the design is coherent with a gradient-phase design which rely on the increase in the effective mode index along the direction transverse to the pillars.

In Fig. 10, we show the meshed model of the optimized design corresponding to the parameters
<figure><img src="image_13.png" /><figcaption>Figure 9: Convergence graph obtained with EGO for the metadeflector optimization.</figcaption></figure>

reported in Tab. 4. We also show the diffraction efficiency in transmission of the mode (0, −1) as a function of the wavelength for a TM-polarized wave for the optimized design. We observe that the design efficiency remains above 40% over the [500, 750] nm regime, with an efficiency at λ = 600 nm above 90%.

<table><caption>Table 4: Optimized parameters for the metadeflector with rectangular pillars.</caption><thead><tr><th rowspan="2">h(nm)</th><th rowspan="2">L<sub>1</sub>(nm)</th><th rowspan="2">L<sub>2</sub>(nm)</th><th rowspan="2">L<sub>3</sub>(nm)</th><th colspan="2">Pillar 1</th><th colspan="2">Pillar 2</th><th colspan="2">Pillar 3</th><th colspan="2">Pillar 4</th></tr><tr><th>d<sub>x</sub>(nm)</th><th>d<sub>y</sub>(nm)</th><th>d<sub>x</sub></th><th>d<sub>y</sub></th><th>d<sub>x</sub></th><th>d<sub>y</sub></th><th>d<sub>x</sub></th><th>d<sub>y</sub></th></tr></thead><tbody><tr><td>625.2</td><td>103.2</td><td>100.3</td><td>115.2</td><td>167.8</td><td>219.5</td><td>135.5</td><td>179.4</td><td>107.9</td><td>145.8</td><td>94.3</td><td>146.0</td></tr></tbody></table>

In Fig. 11, we show the vertical cuts of the norm of the electric and magnetic fields and the real parts of the $E_y$ and $H_x$ components. The field maps were obtained using a finer mesh that the one used during the optimization process (48K cells) and using a $\mathbb{P}_2$ polynomial order. The field maps show the deflection of the incident plane wave along what seems to be the correct expected deflection angle.

### 2.4.4 Multi-objective optimization

In this section, we will attempt to optimize a metasurface using the previous type of design, i.e. an infinitely repeating pattern of 4 rectangular pillars, in order to achieve a different optical function. We consider the optimization of a beam-splitting metasurface design, deflecting the light in the $y-z$ vertical plane along two symmetric directions with respect to $z$. To achieve this optimization, we use a multi-objective EGO method implemented in diogenes-optim in order to maximize the diffraction efficiencies in transmission of both (0, +1) and (0, −1) modes corresponding to deflection angles of $\theta = \pm 24^\circ$. The multi-objective optimization process carries the simultaneous minimization of the two following objective functions:

$$ \begin{cases} f_1(\xi) = 1 - \eta^{\text{TM}}(0, -1)(\lambda = 600 \text{ nm}) \\ f_2(\xi) = 1 - \eta^{\text{TM}}(0, +1)(\lambda = 600 \text{ nm}) \end{cases} $$
<figure><img src="image_14.png" /><figcaption>**Figure 10:** Optimized design obtained with EGO and its diffraction efficiency of mode (0, −1).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 11:** Vertical cut of the electric and magnetic fields obtained for the optimized design found with EGO.</figcaption></figure>

The result of the multi-objective optimization is the set of the optimal designs representing a trade-off between the different objectives (in our case $f_1$ and $f_2$). This set of points is called the Pareto front.

We use the same bounds on the parameter space as for the mono-objective optimization (cf. Tab. 3) and we use the same numerical parameters for the DGTD simulations with diogenes-dgtd, i.e. a coarse mesh (~3K cells) with $P_2$ polynomial order. We initialize the optimization with a DOE of 80 samples and we run the EGO method for a total of 230 iterations, totalling a runtime of 20 hours. The convergence is illustrated in Fig. 12, where we show the evolution of the Pareto front hypervolume, which is computed as the volume in the objective function space dominated by the front relative to a given reference point. This is an indicator of the quality of the optimization result, with a larger hypervolume meaning a front closer to the optimized values of the objective functions.

From the set of optimized design provided by the Pareto front, we find a design with almost equal diffraction efficiency for both modes ($\eta \approx 41\%$ for both modes), reached by EGO at iteration
<figure><img src="image_16.png" /><figcaption>**Figure 12:** Convergence graph obtained during the bi-objective optimization of the beam-splitter.</figcaption></figure>

180. We report the values of the optimized design parameters in Tab. 5. In Fig. 13, we also show the meshed model of the optimized design and the diffraction efficiency in transmission of the design for the (0, −1) and (0, +1) modes. The design appears to evidence an alternation between pillars wider along x and pillars wider along y, except for the last pillar, which has a much larger y-dimension than the others (its maximum bound of 260 nm). The wavelength dependence of the diffraction efficiencies show that the efficiency of mode (0, −1) reaches a maximum of ≈ 60% in the low wavelength regime then decreases rapidly for high wavelengths, while conversely the efficiency of the (0, +1) mode reaches a maximum of ≈ 75% in the high wavelength regime and vanishes towards λ = 500 nm. This shows that the design splits equally an incident beam into the two modes only at λ = 600 nm at the intersections of the two curves.

<table><caption>**Table 5:** Optimized parameters for the beam-splitter with rectangular pillars.</caption><thead><tr><th rowspan="2">h(nm)</th><th rowspan="2">L<sub>1</sub>(nm)</th><th rowspan="2">L<sub>2</sub>(nm)</th><th rowspan="2">L<sub>3</sub>(nm)</th><th colspan="2">Pillar 1</th><th colspan="2">Pillar 2</th><th colspan="2">Pillar 3</th><th colspan="2">Pillar 4</th></tr><tr><th>d<sub>x</sub>(nm)</th><th>d<sub>y</sub>(nm)</th><th>d<sub>x</sub></th><th>d<sub>y</sub></th><th>d<sub>x</sub></th><th>d<sub>y</sub></th><th>d<sub>x</sub></th><th>d<sub>y</sub></th></tr></thead><tbody><tr><td>665.7</td><td>83.2</td><td>89.9</td><td>50.0</td><td>183.1</td><td>139.2</td><td>90.0</td><td>184.6</td><td>197.6</td><td>120.5</td><td>191.5</td><td>260.0</td></tr></tbody></table>

In Fig. 14, we show the vertical cuts of the norm of **E** and **H** and the *E*<sub>y</sub> and *H*<sub>x</sub> components. The *E*<sub>y</sub> and *H*<sub>x</sub> components show that the transmitted wave is split into two waves phase-shifted by ≈ π, inducing a propagation of two wavefronts along opposite angles.
<figure><img src="image_17.png" /><figcaption>**Figure 13:** Optimized design obtained with EGO and diffraction efficiencies of modes (0, −1) and (0, 1)</figcaption></figure>

## 2.5 Conclusion

Through the optimization of a simple four rectangular pillars periodic structure proposed in [2], we found a design of metadeflector reaching 90% deflection efficiency at 600 nm for a TM-polarized plane wave. This results show a slight improvement from the maximum efficiency of 88% found in [2]. However, to rigorously validate this result, we would need to investigate further the accuracy of the numerical solution to make sure it is not affected by nonphysical numerical effects. One of these effects could be undesirable reflections from the PML layers, which would distort the result.

Using bi-objective optimization, we found a design for a beam-splitting metasurface with an efficiency of ≈ 82% (41% for each mode) using the same four pillars structure.

# Part II

## Numerical modeling of wave propagation
in time-varying media

### 1 Motivation

#### 1.1 Active/time-modulated metasurfaces

The first part of this report was focused on the study and design of metasurfaces through the use of numerical modeling and optimization tools. However, to be more precise, the metasurfaces described in the previous part are referred in the literature as *passive* metasurfaces, which were the first to be conceived and are still the most widespread today. Passive metasurfaces are mainly characterized by their fixed optical response, tuned during fabrication by the geometrical properties of its nano-resonators. However more recently, the concept of *active* or *time-modulated* metasurface has been introduced (see [3] for a review on active metasurfaces) which offer new opportunities for research or industrial applications due to their interesting physical properties. Active metasurfaces are characterized by having a dynamically controlled optical response, obtained by exerting an external stimulus on the optical components. This behavior is made
<figure><img src="image_18.png" /><figcaption>**Figure 14:** Vertical cut of field maps obtained for best beam splitter design found by EGO.</figcaption></figure>

possible through the use of active materials in the fabrication of the metasurface, which are able to change properties under an electric current or a temperature change.

This active tuning can be done using liquid crystals [4], phase-change materials through the thermo-optical effect [5], materials with an electro-optical response [6], or leveraging the field effect [7].

<figure><img src="image_19.png" /><figcaption>**Figure 15:** Different types of modulation regimes and physical mechanisms in active metasurfaces. Reproduced from [8].</figcaption></figure>

The timescale on which the metasurface is dynamically modulated is a crucial parameter for engineering its optical function. When the modulation speed is sufficiently slower than the frequency of the incident wave, the metasurface is said to operate in the *quasi-static* regime (Fig. 16a). When the modulation speed exceeds the spectral linewidth of the wave (around 1 MHz for an optical laser), the metasurface operates in the *time-modulated* regime (Fig. 16b). In particular, when the metasurface elements are modulated individually at high frequency, the metasurface is in the *space-time modulated* regime (Fig. 16c).
<figure><img src="image_20.png" /><figcaption>**Figure 16:** Modulation regimes of an electrically controlled metasurface. In the quasi-static regime (a), the modulation speed is too slow to change the frequency of the reflected wave. However, in the time-modulated regime (b), the metasurface elements are modulated at a high frequency $f_m$, leading to the generation of higher order harmonics $f_0 + nf_m$ in the reflected light. In the space-time modulated regime (c), a phase offset is added between the modulation signals of the individual metasurface elements, leading to a phase-gradient in the generated higher order harmonics. Reproduced from [9].</figcaption></figure>

Among the physical phenomena observed in high-frequency time-modulated metasurfaces, we can cite the frequency mixing property where higher-order frequency harmonics are generated (as illustrated in Fig. 16b and Fig. 16c). This can allow for a finer control of the light through manipulation of the individual harmonics. It has been shown that introducing temporal modulation can also break the Lorentz reciprocity constraint ([10]) leading to a non-symmetrical response under time-reversal. The nonreciprocal response of space-time metasurfaces can be exploited to achieve optical power isolation [11].

## 1.2 Existing methods for the modeling of wave propagation in time-varying media

In order to design realistic and efficient time-modulated or space-time modulated metasurfaces, it is essential to model numerically the interaction between electromagnetic waves and materials with a time-varying response. However, the literature dedicated to this subject is quite rare and only exploits a handful of numerical methods.

Most of the work dedicated to the numerical modeling of wave propagation in time-varying media is based on the Finite-Difference Time-Domain (FDTD) method [13], also named Yee's method. We can trace back to [14] the first attempt to develop a full-wave numerical scheme adapted to materials with independently time-varying properties, i.e. where the properties can be modulated in time independently of the electromagnetic field values. The numerical scheme is adapted from the FDTD method in the context of isotropic, linear and non-dispersive media, with a time-varying permittivity. The method has been applied to investigate the properties of simple space-time modulated structures [15, 16], as well as more complex structures such as space-time periodic diffraction gratings [12]. The FDTD method is a simple and efficient numerical method but it has also some drawbacks. Since the method solves the Maxwell equations on a structured grid, it can lose accuracy when dealing with complex geometries. Furthermore, particularly in the context of space-time modulated structures, simulations can be computationally costly due to the multiscale features in the spatial and temporal domains. This can be due to the important time-scale difference between the modulation frequency and the optical frequency of the wave.

We mention also the semi-analytical space-time coupled-wave analysis (RSTWA) method [17] which was proposed as an extension of the rigorous coupled-wave analysis (RCWA) method [18].
<figure><img src="image_21.png" /><figcaption>**Figure 17:** (a): Diffraction of a monochromatic wave by a classical grating (**top**) and a spacetime-modulated grating (**bottom**). In the first case, the spatial diffraction orders $-2 \le m \le 2$ share the same frequency ($\omega_0$) with the input wave while in the second case, each spatial harmonic is formed by an infinite number of temporal harmonics of frequency $\omega_0 + n\Omega$, $n \in \mathbb{Z}$ where $\Omega$ is the modulation frequency. (b): Nonreciprocal response of a spacetime-modulated reflective grating: (**top**): Forward problem. (**bottom**): Backward problem with a nonreciprocal response. Reproduced from [12].</figcaption></figure>

to space-time varying structures. The RSTWA method is based on the plane-wave expansion of Maxwell's equations and allows to compute the optical response of structures with a spatial and temporal periodicity. Its main advantage compared to the FDTD method is the very low computational cost due to its semi-analytical nature. However, it is limited to periodic two-dimensional structures with a layer-by-layer variation.

## 1.3 The Discontinuous Galerkin Time-Domain (DGTD) method

In our case, we will approach the problem using the Discontinuous Galerkin Time-Domain (DGTD) method [19], which is widely used in the field of computational electromagnetics. The discontinuous Galerkin method was introduced in 1973 by Reed and Hill [20] to solve the neutron transport equation. It has later been applied to solve the time-domain Maxwell equations. The discontinuous Galerkin (DG) method is similar to the finite element method where the unknowns are approximated on a set of basis functions. The main difference lies in the fact that the global continuity of the approximation is lifted. In this regard, the DG method can be seen as a combination of the finite element and finite volume methods, where the solution is also discontinuous but constant per cell. Due to the discontinuous nature of the solution, the DG method is well-suited to the modeling of wave propagation in heterogeneous media, where we observe natural discontinuities in the electromagnetic fields. Since the weak formulation is local to an individual element, the solving process doesn't require the inversion of a large mass matrix, leading to a
reduced computational cost compared to the finite element method.
However, the formulation given in [19] cannot handle materials with a time-varying response, as
needed to simulate wave interaction with time-modulated metasurfaces. In the following section,
we will thus extend the DGTD formulation to account for time-varying materials, in particular
by considering the permittivity to be time-dependent in the Maxwell equations.
# 2 The DGTD method for Maxwell's equations applied to time-varying media

## 2.1 Maxwell's equations in time-invariant media

First, we present the general form of Maxwell's equations in the absence of electric or magnetic losses. In this case, the time-dependent Maxwell's equations are expressed as the following system of four differential equations:

$$ \left\{ \begin{array}{ll} \nabla \times \mathbf{E} = -\frac{\partial \mathbf{B}}{\partial t} & \text{Maxwell-Faraday} \\ \nabla \times \mathbf{H} = \frac{\partial \mathbf{D}}{\partial t} + \mathbf{J} & \text{Maxwell-Ampère} \\ \nabla \cdot \mathbf{D} = \rho & \text{Maxwell-Gauss for the electric field} \\ \nabla \cdot \mathbf{B} = 0 & \text{Maxwell-Gauss for the magnetic field} \end{array} \right. \qquad (7) $$

where $\mathbf{E}$ is the electric field, $\mathbf{D}$ is the electric flux density, $\mathbf{H}$ is the magnetic field and $\mathbf{B}$ is the magnetic flux density. In general case, we also take in account the density of electric charges $\rho$ and the free electric current density $\mathbf{J}$.

In linear, isotropic and nondispersive materials, the fields ($\mathbf{E}, \mathbf{B}$) and ($\mathbf{D}, \mathbf{H}$) are related by the constitutive relations:

$$ \begin{aligned} \mathbf{D} &= \epsilon(\mathbf{x})\mathbf{E} = \epsilon_0\epsilon_r(\mathbf{x})\mathbf{E} \\ \mathbf{B} &= \mu(\mathbf{x})\mathbf{H} = \mu_0\mu_r(\mathbf{x})\mathbf{H} \end{aligned} $$

where $\epsilon_0$ and $\epsilon_r(\mathbf{x})$ are the vacuum and relative dielectric permittivity of the considered material and $\mu_0$ and $\mu_r(\mathbf{x})$ are the vacuum and relative magnetic permeability of the propagation medium. The permittivity and permeability are assumed to vary in space since we consider an heterogeneous propagation medium.

Using these constitutive relations, we can close the system formed by the two curl equations expressed using the fields $\mathbf{E}$ and $\mathbf{H}$:

$$ \begin{aligned} \nabla \times \mathbf{E} &= -\mu_0\mu_r(\mathbf{x})\frac{\partial \mathbf{H}}{\partial t} \\ \nabla \times \mathbf{H} &= \epsilon_0\epsilon_r(\mathbf{x})\frac{\partial \mathbf{E}}{\partial t} + \mathbf{J} \end{aligned} $$

We can further simplify the equations by removing the constants $\epsilon_0$ and $\mu_0$ through normalization, as detailed in [19].

Introducing the new variables $\tilde{\mathbf{H}} = \sqrt{\frac{\mu_0}{\epsilon_0}}\mathbf{H}$, $\tilde{\mathbf{E}} = \mathbf{E}$, $\tilde{t} = c_0t$, and $\tilde{\mathbf{J}} = \sqrt{\frac{\mu_0}{\epsilon_0}}\mathbf{J}$ we obtain the following normalized system:

$$ \begin{aligned} \mu_r(\mathbf{x})\frac{\partial \mathbf{H}}{\partial t} &= -\nabla \times \mathbf{E} \\ \epsilon_r(\mathbf{x})\frac{\partial \mathbf{E}}{\partial t} &= \nabla \times \mathbf{H} - \mathbf{J} \end{aligned} \qquad (8) $$

## 2.2 Maxwell's equations in spacetime-varying media

In the following, we will now assume that the permittivity of the propagation medium can vary spatially as well as temporally in the physical domain and is described by the function $\epsilon_r(\mathbf{x}, t)$.
Furthermore, we suppose that the medium is non-magnetic ($\mu_r = 1$). With this new assumption, the permittivity cannot be taken out of the time-derivative in Maxwell-Ampère's equation. We can either keep the permittivity inside the time-derivative, leading to a system involving the unknowns $\mathbf{E}$ and $\mathbf{H}$, as in the case of a time-independent medium:

$$ 
\begin{aligned}
\frac{\partial \mathbf{H}}{\partial t} &= -\nabla \times \mathbf{E} \\
\frac{\partial (\epsilon_r(\mathbf{x}, t)\mathbf{E})}{\partial t} &= \nabla \times \mathbf{H} - \mathbf{J}
\end{aligned}
\qquad (9) $$

We can also reformulate the two equations by using the definition of the electric displacement $\mathbf{D} = \epsilon_0\epsilon_r(\mathbf{x}, t)\mathbf{E}$, leading to a system involving the unknowns $\mathbf{D}$ and $\mathbf{H}$:

$$ 
\begin{aligned}
\frac{\partial \mathbf{H}}{\partial t} &= -\nabla \times (\epsilon_r^{-1}(\mathbf{x}, t)\mathbf{D}) \\
\frac{\partial \mathbf{D}}{\partial t} &= \nabla \times \mathbf{H} - \mathbf{J}
\end{aligned}
\qquad (10) $$

The first option is the simplest to implement because the permittivity variation will only be treated in the time-stepping algorithm which is more straightforward than the discontinuous Galerkin scheme. It also requires less modifications compared to an implementation of the DGTD scheme in time-invariant media.

The second option can however present some advantages. If we consider that the permittivity can vary spatially inside each cell, handling the permittivity in the spatial part of the scheme is more adequate since we can account for this spatial dependence in the computation of the DG matrices.

However, we will only consider in the following that the permittivity is spatially constant per cell. We therefore choose the first form of the system (Eq. 9) for the derivation of the numerical scheme.

## 2.3 Weak formulation

We will now write the weak formulation of the Maxwell's system for time-varying media.

Let $\Omega \subset \mathbb{R}^3$ be the domain of the problem, assumed to be bounded and convex, and delimited by its boundary $\partial\Omega$. We note $\Omega_h$ the discretization of $\Omega$ by the triangulation $\mathcal{T}_h$ composed of simplicial mesh elements $(T_i)_{i \in [1,N]}$ where $N$ is the total number of elements. The internal faces between two adjacent cells $T_i$ and $T_k$ are denoted by $a_{ik} = T_i \cap T_k$ and the unit normal vector to $a_{ik}$ oriented from $T_i$ to $T_k$ is denoted by $\mathbf{n}_{ik}$. Finally, for each cell $T_i$, we define $\mathcal{V}_i$ as the set of neighbouring cells indices which share a triangular face with $T_i$.

In order to write the local weak formulation of the problem in the cell $T_i$, we multiply the equations by a vector test function $\psi$ then integrate each term over $T_i$, yielding:

$$ 
\begin{aligned}
\int_{T_i} \frac{\partial \mathbf{H}}{\partial t} \cdot \psi + \int_{T_i} \nabla \times \mathbf{E} \cdot \psi &= 0 \\
\int_{T_i} \frac{\partial (\epsilon_r(\mathbf{x}, t)\mathbf{E})}{\partial t} \cdot \psi - \int_{T_i} \nabla \times \mathbf{H} \cdot \psi &= - \int_{T_i} \mathbf{J} \cdot \psi
\end{aligned}
$$

Using Green's formula and vector calculus, we can reformulate these equations so as to introduce an additional term over the cell boundary $\partial T_i$ for each equation:
$$ 
\begin{aligned}
& \int_{T_i} \frac{\partial \mathbf{H}}{\partial t} \cdot \psi + \int_{T_i} \mathbf{E} \cdot \nabla \times \psi = \int_{\partial T_i} (\mathbf{E} \times \mathbf{n}_i) \cdot \psi \\
& \int_{T_i} \frac{\partial (\epsilon_r(\mathbf{x}, t)\mathbf{E})}{\partial t} \cdot \psi - \int_{T_i} \mathbf{H} \cdot \nabla \times \psi = - \int_{T_i} \mathbf{J} \cdot \psi - \int_{\partial T_i} (\mathbf{H} \times \mathbf{n}_i) \cdot \psi
\end{aligned}
\quad (11) $$

where $\mathbf{n}_i$ is the unit normal vector to the cell boundary $\partial T_i$.

In order for the integrals terms over $T_i$ to be defined, we take $\mathbf{E}$ in $H_0(\text{curl}, \Omega_h)$ and $\mathbf{H}$ in $H(\text{curl}, \Omega_h)$. The choice of Sobolev space for $\mathbf{E}$ is to verify the coercivity property in order to ensure the existence and unicity of the solution according to the Lax-Milgram theorem. Furthermore, the boundary term requires the existence of the trace of $\psi$ on $\partial T_i$. We therefore take $\psi$ in $H^1(\Omega_h)$.

## 2.4 Space discretization

We first need to define the approximation space in which we will define the semi-discrete fields. We take the space of polynomials of order less than or equal to $p$ defined piecewise on each element $T_i$ of the triangulation $\mathcal{T}_h$.

$$ V_h = \{v \in (L^2(\Omega))^3, v|_{T_i} \in (\mathbb{P}_p(T_i))^3 \forall T_i \in \mathcal{T}_h\} $$

The semi-discrete fields defined in $V_h$ are denoted by $(\mathbf{E}_h, \mathbf{H}_h)$ from which we define their restriction on each element $T_i$: $(\mathbf{E}_i, \mathbf{H}_i) = (\mathbf{E}_h|_{T_i}, \mathbf{H}_h|_{T_i})$. Each component of the latter approximate vector fields will be sought for as a linear combination of scalar basis functions $(\phi_{ik})_{1 \le k \le d_i}$, where $d_i$ denotes the local number of degrees of freedom for each element $T_i$. In order to have more compact formulations, we also define the following three vector:

$$ \phi_{ik}^x = \begin{bmatrix} \phi_{ik} \\ 0 \\ 0 \end{bmatrix}, \quad \phi_{ik}^y = \begin{bmatrix} 0 \\ \phi_{ik} \\ 0 \end{bmatrix}, \quad \phi_{ik}^z = \begin{bmatrix} 0 \\ 0 \\ \phi_{ik} \end{bmatrix} $$

For the $\mathbf{E}_i$ approximate field restricted to $T_i$, the local expansion of a given component over the scalar basis functions reads:

$$ E_i^v = \sum_{j=1}^{d_i} E_{ij}^v \phi_{ij}, \quad v \in \{x, y, z\} $$

## 2.5 Numerical fluxes

To obtain the semi-discrete form of the weak formulation, we replace the exact fields $(\mathbf{E}, \mathbf{H})$ by the approximate fields $(\mathbf{E}_h, \mathbf{H}_h)$ in Eqs. 11. For the integrals over $T_i$, the latter can be replaced by their restrictions $(\mathbf{E}_i, \mathbf{H}_i)$.

The integrals over $\partial T_i$ can be decomposed as a sum of integrals over the internal faces $a_{il}$ between $T_i$ and its neighbouring cells $T_l$, $l \in \mathcal{V}_i$:

$$ \int_{\partial T_i} (\mathbf{E}_h \times \mathbf{n}_i) \cdot \psi = \sum_{l \in \mathcal{V}_i} \int_{a_{il}} (\mathbf{E}_h \times \mathbf{n}_{il}) \cdot \psi $$

$$ \int_{\partial T_i} (\mathbf{H}_h \times \mathbf{n}_i) \cdot \psi = \sum_{l \in \mathcal{V}_i} \int_{a_{il}} (\mathbf{H}_h \times \mathbf{n}_{il}) \cdot \psi $$
However, since the fields are allowed to be discontinuous at the interfaces between elements, the integrals over $a_{il}$ are not uniquely defined since $\mathbf{E}_h|_{a_{il}}$ for example, can relate to either its value in $T_i$ or $T_l$. Likewise, since the permittivity can also be discontinuous at interfaces, $\epsilon_r$ can relate to either $\epsilon_i$ or $\epsilon_l$. We introduce a numerical flux to properly define these integrals and connect the neighbouring cells in the weak formulation. We hence replace these integrals by:

$$ \int_{a_{il}} (\mathbf{E}_* \times \mathbf{n}_{il}) \cdot \psi, \int_{a_{il}} (\mathbf{H}_* \times \mathbf{n}_{il}) \cdot \psi $$

where the definition of $\mathbf{E}_*$ and $\mathbf{H}_*$ depends on the choice of numerical flux. In the following, we choose a centered flux approximation expressed as:

$$ \mathbf{E}_* = \frac{\mathbf{E}_i + \mathbf{E}_l}{2}, \quad \mathbf{H}_* = \frac{\mathbf{H}_i + \mathbf{H}_l}{2} $$

This flux approximation is non-dissipative and leads to a $L^2$ spatial convergence in $h^p$ (see [21] for proof).

In the following, we will simplify the problem by assuming that the permittivity is spatially constant per cell, while remaining discontinuous at interfaces in the general case. In the cell $T_i$, we thus have $\epsilon_r(\mathbf{x}, t) = \epsilon_i(t)$. We also replace the electric current density by its restriction $\mathbf{J}_i$. We obtain the following discrete weak DG formulation:

$$ \begin{aligned} & \int_{T_i} \frac{\partial \mathbf{H}_i}{\partial t} \cdot \psi + \int_{T_i} \mathbf{E}_i \cdot \nabla \times \psi = \sum_{l \in \mathcal{V}_i} \int_{a_{il}} (\mathbf{E}_* \times \mathbf{n}_{il}) \cdot \psi \\ & \int_{T_i} \frac{\partial(\epsilon_i(t)\mathbf{E}_i)}{\partial t} \cdot \psi - \int_{T_i} \mathbf{H}_i \cdot \nabla \times \psi = - \sum_{l \in \mathcal{V}_i} \int_{a_{il}} (\mathbf{H}_* \times \mathbf{n}_{il}) \cdot \psi - \int_{T_i} \mathbf{J}_i \cdot \psi \end{aligned} \quad (12) $$

## 2.6 DG matrix-vector formulation

Replacing $\psi$ by our choice of vector base functions $\phi_{ik}^v$, we obtain the semi-discrete DG scheme for our problem:

$$ \begin{aligned} & \int_{T_i} \frac{\partial \mathbf{H}_i}{\partial t} \cdot \phi_{ik}^v + \int_{T_i} \mathbf{E}_i \cdot \nabla \times \phi_{ik}^v = \sum_{l \in \mathcal{V}_i} \int_{a_{il}} (\mathbf{E}_* \times \mathbf{n}_{il}) \cdot \phi_{ik}^v \\ & \int_{T_i} \frac{\partial(\epsilon_i(t)\mathbf{E}_i)}{\partial t} \cdot \phi_{ik}^v - \int_{T_i} \mathbf{H}_i \cdot \nabla \times \phi_{ik}^v = - \sum_{l \in \mathcal{V}_i} \int_{a_{il}} (\mathbf{H}_* \times \mathbf{n}_{il}) \cdot \phi_{ik}^v - \int_{T_i} \mathbf{J}_i \cdot \phi_{ik}^v \end{aligned} \quad (13) $$

We will rewrite the DG semi-discrete scheme under a matrix-vector formulation.

**Mass matrix** We first reformulate the time-derivative term of the two equations to derive the mass matrix of the system. For the $x$ component of the $\mathbf{E}$ evolution equation (Maxwell-Ampère), using the expansion of the field over the set of basis functions yields:
$$
\begin{align*}
\int_{T_i} \frac{\partial(\epsilon_i(t)\mathbf{E}_i)}{\partial t} \cdot \phi_{ik}^x &= \int_{T_i} \frac{\partial(\epsilon_i(t)E_i^x)}{\partial t} \phi_{ik} \\
&= \int_{T_i} \frac{\partial}{\partial t} \left[ \epsilon_i(t) \sum_{j=1}^{d_i} E_{ij}^x \phi_{ij} \right] \cdot \phi_{ik} \\
&= \sum_{j=1}^{d_i} \frac{\partial}{\partial t} \epsilon_i(t) E_{ij}^x \int_{T_i} \phi_{ij} \phi_{ik} \\
&= \left( \mathbb{M}_i \frac{\partial(\epsilon_i(t)\bar{\mathbf{E}}_i^x)}{\partial t} \right)_k
\end{align*}
$$

where $\mathbb{M}_i$ is a square matrix of dimension $d_i \times d_i$ defined by:

$$ (\mathbb{M}_i)_{jk} = \int_{T_i} \phi_{ij} \phi_{ik} $$

and where we store the components $E_{ij}^v$ in the set of three vectors of dimension $d_i$:

$$ \bar{\mathbf{E}}_i^v = \begin{bmatrix} E_{i1}^v \\ \vdots \\ E_{id_i}^v \end{bmatrix}, \quad v \in \{x, y, z\} $$

For the **H** evolution equation (Maxwell-Faraday), we have:

$$ \int_{T_i} \frac{\partial \mathbf{H}_i}{\partial t} \cdot \phi_{ik}^x = \left( \mathbb{M}_i \frac{\partial \bar{\mathbf{H}}_i^x}{\partial t} \right)_k \quad (14) $$

**Stiffness matrix** We now derive the stiffness matrix of the DG formulation from the curl integrals of Eqs. 13. Focusing on the *x* component of the **H** evolution equation (curl term with **E**), we have:

$$
\begin{align*}
\int_{T_i} \mathbf{E}_i \cdot \nabla \times \phi_{ik}^x &= \int_{T_i} \left( E_i^y \frac{\partial \phi_{ik}}{\partial z} - E_i^z \frac{\partial \phi_{ik}}{\partial y} \right) \\
&= \int_{T_i} \sum_{j=1}^{d_i} \left( E_{ij}^y \phi_{ij} \frac{\partial \phi_{ik}}{\partial z} - E_{ij}^z \phi_{ij} \frac{\partial \phi_{ik}}{\partial y} \right) \\
&= \sum_{j=1}^{d_i} E_{ij}^y \int_{T_i} \phi_{ij} \frac{\partial \phi_{ik}}{\partial z} - \sum_{j=1}^{d_i} E_{ij}^z \int_{T_i} \phi_{ij} \frac{\partial \phi_{ik}}{\partial y} \\
&= (\mathbb{K}_i^z \bar{\mathbf{E}}_i^y - \mathbb{K}_i^y \bar{\mathbf{E}}_i^z)_k \\
&= -(\bar{\mathbb{K}}_i \times \bar{\mathbf{E}}_i)_k^x
\end{align*}
$$

where the stiffness matrices are defined by:

$$ (\mathbb{K}_i^v)_{jk} = \int_{T_i} \phi_{ij} \frac{\partial \phi_{ik}}{\partial v}, \quad v \in \{x, y, z\} $$

and the general $3d_i \times d_i$ stiffness matrix is defined by
$$ \overline{\mathbb{K}}_i = \begin{bmatrix} \mathbb{K}_i^x \\ \mathbb{K}_i^y \\ \mathbb{K}_i^z \end{bmatrix} $$

We also defined the general dimension $3d_i$ vector regrouping the three vectors $\overline{\mathbb{E}}_i^v$ as:

$$ \overline{\mathbb{E}}_i = \begin{bmatrix} \overline{\mathbb{E}}_i^x \\ \overline{\mathbb{E}}_i^y \\ \overline{\mathbb{E}}_i^z \end{bmatrix} $$

For the $\mathbb{E}$ evolutionary equation, we have similarly:

$$ \int_{T_i} \mathbf{H}_i \cdot \nabla \times \phi_{\mathbf{ik}}^x = -(\overline{\mathbb{K}}_i \times \overline{\mathbb{H}}_i)_k^x $$

**Flux matrix** Finally, we reformulate the surface integrals over the internal faces $a_{il}$ to derive the flux matrices. We use in the derivation the centered flux approximation. We also exploit the property that in a conforming mesh, we can equivalently expand a field defined on the interface $a_{il}$ using the local basis functions of $T_i$ or $T_l$.

Focusing on the $x$ component of the $\mathbf{H}$ evolutionary equation, we have:

$$ \begin{aligned} \int_{a_{il}} (\mathbf{E}_* \times \mathbf{n}_{il}) \cdot \phi_{\mathbf{ik}}^x &= \int_{a_{il}} (E_*^y n_{il}^z - E_*^z n_{il}^y) \phi_{ik} \\ &= \int_{a_{il}} \left( \frac{E_*^y + E_*^z}{2} n_{il}^z - \frac{E_*^y - E_*^z}{2} n_{il}^y \right) \phi_{ik} \\ &= \int_{a_{il}} \sum_{j=1}^{d_i} (E_{*,j}^y n_{il}^z - E_{*,j}^z n_{il}^y) \phi_{ij} \phi_{ik} \\ &= \sum_{j=1}^{d_i} (E_{*,j}^y n_{il}^z - E_{*,j}^z n_{il}^y) \int_{a_{il}} \phi_{ij} \phi_{ik} \\ &= (\mathbb{S}_{il} (\overline{\mathbb{E}}_* \times \mathbf{n}_{il}))_k^x \end{aligned} $$

where the flux matrices are defined by:

$$ (\mathbb{S}_{il})_{jk} = \int_{a_{il}} \phi_{ij} \phi_{ik} $$

For the $\mathbb{E}$ evolution equation, we have similarly:

$$ \int_{a_{il}} (\mathbf{H}_* \times \mathbf{n}_{il}) \cdot \phi_{\mathbf{ik}}^x = (\mathbb{S}_{il} (\overline{\mathbf{H}}_* \times \mathbf{n}_{il}))_k^x $$

In order to write the semi-discrete DG formulation in the most compact manner we extend the mass and flux matrices to $3d_i \times 3d_i$ matrices defined by:

$$ \overline{\mathbb{M}}_i = \begin{bmatrix} \mathbb{M}_i & \mathbf{0}_{d_i \times d_i} & \mathbf{0}_{d_i \times d_i} \\ \mathbf{0}_{d_i \times d_i} & \mathbb{M}_i & \mathbf{0}_{d_i \times d_i} \\ \mathbf{0}_{d_i \times d_i} & \mathbf{0}_{d_i \times d_i} & \mathbb{M}_i \end{bmatrix}, \quad \overline{\mathbb{S}}_{il} = \begin{bmatrix} \mathbb{S}_{il} & \mathbf{0}_{d_i \times d_i} & \mathbf{0}_{d_i \times d_i} \\ \mathbf{0}_{d_i \times d_i} & \mathbb{S}_{il} & \mathbf{0}_{d_i \times d_i} \\ \mathbf{0}_{d_i \times d_i} & \mathbf{0}_{d_i \times d_i} & \mathbb{S}_{il} \end{bmatrix} \quad (15) $$

From these matrices we express the matrix-vector formulation of the semi-discrete DG scheme for Maxwell equations in a time-varying medium:
$$ 
\begin{aligned}
\overline{\mathbf{M}}_i \frac{\partial \overline{\mathbf{H}}_i}{\partial t} &= \overline{\mathbf{K}}_i \times \overline{\mathbf{E}}_i + \sum_{l \in \mathcal{V}_i} \overline{\mathbf{S}}_{il} (\overline{\mathbf{E}}_* \times \mathbf{n}_{il}) \\
\overline{\mathbf{M}}_i \frac{\partial (\epsilon_i(t) \overline{\mathbf{E}}_i)}{\partial t} &= \overline{\mathbf{K}}_i \times \overline{\mathbf{H}}_i - \sum_{l \in \mathcal{V}_i} \overline{\mathbf{S}}_{il} (\overline{\mathbf{H}}_* \times \mathbf{n}_{il}) - \overline{\mathbf{M}}_i \overline{\mathbf{J}}_i
\end{aligned}
\qquad (16) $$

## 2.7 Time discretization

For the time-stepping scheme, we choose the order 2 Leap Frog scheme:

$$ \phi^{n+1} = \phi^n + \Delta t f(t_{n+1/2}, \phi^{n+1/2}) $$

where $\phi^n$ denotes the discrete approximation of $\phi(t_n)$ with $t_n = n\Delta t$.

Since we consider a multiple variables problem (we solve for $\mathbf{E}$ and $\mathbf{H}$), we need to sample the latter fields on staggered time grid in order for the scheme to remain explicit. Thus, we discretize $\mathbf{E}$ at $t_n = n\Delta t$ and $\mathbf{H}$ and $\mathbf{J}$ at $t_{n+1/2} = (n + \frac{1}{2})\Delta t$. Furthermore, we discretize the permittivity at $t_n$ since it is related to $\mathbf{E}$ by $\mathbf{D}(t_n) = \epsilon_r(t_n)\mathbf{E}(t_n)$.

We finally obtain the fully discretized LF2-DG scheme for Maxwell's equations in media with a spatially constant per cell and time-dependent permittivity:

$$ 
\begin{aligned}
\overline{\mathbf{M}}_i \frac{\epsilon_i^{n+1} \overline{\mathbf{E}}_i^{n+1} - \epsilon_i^n \overline{\mathbf{E}}_i^n}{\Delta t} &= \overline{\mathbf{K}}_i \times \overline{\mathbf{H}}_i^{n+1/2} - \sum_{l \in \mathcal{V}_i} \overline{\mathbf{S}}_{il} (\overline{\mathbf{H}}_*^{n+1/2} \times \mathbf{n}_{il}) - \overline{\mathbf{J}}_i (t_{n+1/2}) \\
\overline{\mathbf{M}}_i \frac{\overline{\mathbf{H}}_i^{n+3/2} - \overline{\mathbf{H}}_i^{n+1/2}}{\Delta t} &= -\overline{\mathbf{K}}_i \times \overline{\mathbf{E}}_i^{n+1} + \sum_{l \in \mathcal{V}_i} \overline{\mathbf{S}}_{il} (\overline{\mathbf{E}}_*^{n+1} \times \mathbf{n}_{il})
\end{aligned}
\qquad (17) $$

## 2.8 Reduction to two dimensions

In the context of this study, we will for now only consider problems defined in a two-dimensional domain, allowing for reduced computational costs compared to a 3D implementation of the scheme. The first step is to adapt the three-dimensional Maxwell equations to the propagation in a two-dimensional computational domain.

We consider that the propagation domain is infinite and uniform in the z-direction and the electromagnetic fields are also uniform in the z-direction (all the partial derivatives of the fields with respect to z vanish).

Following these conditions, Maxwell's curl equations (Eqs. 9) reduce to two independent sets of three equations known as the transverse-magnetic mode (TMz), which involve ($H_x, H_y, E_z$) and the transverse-electric mode (TEz), which involve ($E_x, E_y, H_z$):

$$ (\text{TM}_z): \left\{ \begin{array}{l} \frac{\partial H_x}{\partial t} = -\frac{\partial E_z}{\partial y} \\ \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} \\ \frac{\partial (\epsilon_r(\mathbf{x}, t)E_z)}{\partial t} = \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y} - J_z \end{array} \right. \qquad (\text{TE}_z): \left\{ \begin{array}{l} \frac{\partial (\epsilon_r(\mathbf{x}, t)E_x)}{\partial t} = \frac{\partial H_z}{\partial y} - J_x \\ \frac{\partial (\epsilon_r(\mathbf{x}, t)E_y)}{\partial t} = -\frac{\partial H_z}{\partial x} - J_y \\ \frac{\partial H_z}{\partial t} = \frac{\partial E_x}{\partial y} - \frac{\partial E_y}{\partial x} \end{array} \right. \qquad (19) $$
Making the assumption that the permittivity is spatially constant per cell with $\epsilon_r(\mathbf{x}, t) = \epsilon_i(t)$ in cell $T_i$, and applying the same manipulations as for the three-dimensional scheme, we obtain the fully discretized LF2-DG scheme for the two-dimensional TMz and TEz modes.

## LF2-DG scheme for the Transverse-Magnetic (TMz) mode

$$ 
\begin{aligned}
M_i \frac{\epsilon_i^{n+1} \bar{\mathbf{E}}_i^{z,n+1} - \epsilon_i^n \bar{\mathbf{E}}_i^{z,n}}{\Delta t} &= \left( \mathbb{K}_i^y \bar{\mathbf{H}}_*^{x,n+1} - \mathbb{K}_i^x \bar{\mathbf{H}}_*^{y,n+1} \right) + \sum_{l \in \mathcal{V}_i} \mathbb{S}_{il} \left( \bar{\mathbf{H}}_*^{y,n+1} n_{il}^x - \bar{\mathbf{H}}_*^{x,n+1} n_{il}^y \right) - \bar{\mathbf{J}}_i^z (t_{n+1/2}) \\
M_i \frac{\bar{\mathbf{H}}_i^{x,n+3/2} - \bar{\mathbf{H}}_i^{x,n+1/2}}{\Delta t} &= \mathbb{K}_i^y \bar{\mathbf{E}}_i^{z,n+1/2} - \sum_{l \in \mathcal{V}_i} \left( \mathbb{S}_{il} \bar{\mathbf{E}}_*^{z,n+1/2} \right) n_{il}^y \\
M_i \frac{\bar{\mathbf{H}}_i^{y,n+3/2} - \bar{\mathbf{H}}_i^{y,n+1/2}}{\Delta t} &= -\mathbb{K}_i^x \bar{\mathbf{E}}_i^{z,n+1/2} + \sum_{l \in \mathcal{V}_i} \left( \mathbb{S}_{il} \bar{\mathbf{E}}_*^{z,n+1/2} \right) n_{il}^x
\end{aligned}
\quad (20) $$

## LF2-DG scheme for the Transverse-Electric (TEz) mode

$$ 
\begin{aligned}
M_i \frac{\epsilon_i^{n+1} \bar{\mathbf{E}}_i^{z,n+1} - \epsilon_i^n \bar{\mathbf{E}}_i^{z,n}}{\Delta t} &= -\mathbb{K}_i^y \bar{\mathbf{H}}_*^{z,n+1/2} + \sum_{l \in \mathcal{V}_i} \left( \mathbb{S}_{il} \bar{\mathbf{H}}_*^{z,n+1/2} \right) n_{il}^y - \bar{\mathbf{J}}_i^x (t_{n+1/2}) \\
M_i \frac{\epsilon_i^{n+1} \bar{\mathbf{E}}_i^{y,n+1} - \epsilon_i^n \bar{\mathbf{E}}_i^{y,n}}{\Delta t} &= \mathbb{K}_i^x \bar{\mathbf{H}}_*^{z,n+1/2} - \sum_{l \in \mathcal{V}_i} \left( \mathbb{S}_{il} \bar{\mathbf{H}}_*^{z,n+1/2} \right) n_{il}^x - \bar{\mathbf{J}}_i^y (t_{n+1/2}) \\
M_i \frac{\bar{\mathbf{H}}_i^{z,n+3/2} - \bar{\mathbf{H}}_i^{z,n+1/2}}{\Delta t} &= \left( \mathbb{K}_i^x \bar{\mathbf{E}}_*^{y,n+1} - \mathbb{K}_i^y \bar{\mathbf{E}}_*^{x,n+1} \right) + \sum_{l \in \mathcal{V}_i} \mathbb{S}_{il} \left( \bar{\mathbf{E}}_*^{x,n+1} n_{il}^y - \bar{\mathbf{E}}_*^{y,n+1} n_{il}^x \right)
\end{aligned}
\quad (21) $$

In the following section, we will only implement and validate the LF2-DG scheme for the Transverse-Magnetic mode (Eqs. 20) but all the derivations and results which follow can be generalized to the Transverse-Electric mode with few modifications.
# 3 Numerical validation

In this section, we validate numerically the DGTD scheme applied to time-varying media by considering a few simple cases involving the temporal modulation of the propagation medium permittivity. Here, we only consider two-dimensional cases in order to validate the scheme simulating the propagation of the 2D Transverse Magnetic (TMz) mode (involving the field components $H_x$, $H_y$ and $E_z$).

## 3.1 Time-modulated rectangular resonant cavity

We first consider the case of a 2D rectangular resonant cavity filled with a medium of spatially homogeneous permittivity but which can vary as a function of time. Hence the medium properties are characterized by the time-dependent permittivity function $\epsilon_r(t)$. We also consider that the medium is non-magnetic ($\mu_r = 1$).

The case where the cavity is filled with a medium of time-invariant permittivity and the one where the medium permittivity is modulated in time (following a sinusoidal profile) are illustrated in Fig. 18a and 18b respectively. The first case corresponds to a well-known classical problem which solutions are the fundamental stationary modes supported by the resonant cavity. The spatial and temporal profiles of the simulated first stationary mode TMz (1, 1) for $\epsilon_r = 1$ are shown in Fig. 19a and the left of Fig. 19b, where the temporal profile is compared to the analytical temporal profile of the corresponding mode.

In the second case, however, the analytical solution for a general permittivity modulation profile $\epsilon_r(t)$ is unknown in the literature and hence we cannot validate directly the numerical solution by comparing the latter with a classical analytical solution.

For a sinusoidal permittivity modulation, the simulated spatial profile is identical to the one obtained for the first case (Fig. 19a) but the temporal profile differs greatly as shown in the right of Fig. 19b.

In the following, we employ two different methods to attempt to validate the latter case: the method of manufactured solutions and the derivation of an analytical solution for a given permittivity variation profile.

<figure><img src="image_22.png" /><figcaption>**Figure 18:** Two-dimensional rectangular resonant cavity: the case where the interior medium has a time-invariant permittivity and a null electric current ($\mathbf{J} = \mathbf{0}$) (a), a sinusoidally time-varying permittivity and a null electric current (b) and a sinusoidally time-varying permittivity and a non-null electric current ($\mathbf{J} \neq \mathbf{0}$) (c).</figcaption></figure>
### 3.1.1 Method of manufactured solutions

In a first attempt to validate numerically the DGTD scheme using the case of a time-modulated resonant cavity, we employ the so-called Method of Manufactured Solutions (MMS). This approach consists in building an artificial analytical solution to the differential equations solved by the numerical scheme by introducing additional source terms in the equations. The resulting solution is often not physically relevant, but allows us to easily validate and benchmark the implementation of the numerical scheme.

In our case, we seek to introduce a source term in the Maxwell system for time-varying media in the form of an electric current density term $J(\mathbf{x},t)$ as to impose the classical stationary modes (22) as an analytical solution. The method is illustrated in Fig. 18c.

The MMS was applied to a one-dimensional resonant cavity filled with a medium of time-varying permittivity in [22] in order to validate another numerical method for solving Maxwell equations in time-varying media. Here, we extend the MMS to a two-dimensional cavity by considering the transverse magnetic mode $(H_x, H_y, E_z)$.
We consider a rectangular resonant cavity of dimensions ($a_x, a_y$), filled with a medium of constant permittivity ($\epsilon_r = \epsilon_{r0}$) and unit permeability ($\mu_r = 1$). In the normalized unit system, the TMz stationary modes supported by the cavity take the form:

$$ 
\left\{
\begin{aligned}
H_x &= -\frac{k_y}{\omega_s} \sin(k_x x) \cos(k_y y) \sin(\omega_s t) \\
H_y &= \frac{k_x}{\omega_s} \cos(k_x x) \sin(k_y y) \sin(\omega_s t) \\
E_z &= \sin(k_x x) \sin(k_y y) \cos(\omega_s t)
\end{aligned}
\right.
\qquad (22) $$

where $k_x = \frac{n_x \pi}{a_x}$ and $k_y = \frac{n_y \pi}{a_y}$ are the components of the wavevector $k$ with ($n_x, n_y$) the integer mode numbers, $\omega_s = c_r \sqrt{k_x^2 + k_y^2}$ is the frequency of the wave with $c_r = \frac{1}{\sqrt{\epsilon_r \mu_r}} = \frac{1}{\sqrt{\epsilon_{r0} \mu_{r0}}}$ the relative speed of light in the cavity.

In the TMz mode, only the z-component of the electric current density appears in the Maxwell equations, where the Maxwell-Ampère equation for a time-varying permittivity $\epsilon_r(t)$ gives the following relation inside the time-modulated cavity:

$$ J_z(x, y, t) = \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y} - \frac{\partial(\epsilon_r(t)E_z)}{\partial t} $$

By replacing the fields components by the stationary modes (22) for a constant permittivity $\epsilon_{r0}$, we obtain an expression of $J_z$ as a function of the modulation profile $\epsilon_r(t)$ and its time derivative.

$$ J_z(x, y, t) = \left[ (\epsilon_r(t) - \epsilon_{r0})\omega_s \sin(\omega_s t) - \frac{\partial \epsilon_r(t)}{\partial t} \cos(\omega_s t) \right] \sin(k_x x) \sin(k_y y) $$

In particular, for a sinusoidal permittivity modulation around a mean value of $\epsilon_{r0}$ expressed by the function $\epsilon_r(t) = \epsilon_{r0}(1 + \delta \sin(\omega_m t))$, we obtain the expression of the term $J_z$ we need to introduce in the cavity to verify the stationary modes.

$$ J_z(x, y, t) = \epsilon_{r0} \delta \sin(k_x x) \sin(k_y y) [\omega_s \sin(\omega_m t) \sin(\omega_s t) - \omega_m \cos(\omega_m t) \cos(\omega_s t)] $$

After implementing this artificial source term in the DGTD solver, we can now validate the code by comparing the numerical solution with the analytical manufactured solution i.e. the stationary modes (22).

In Fig. 20, we represent the time evolution of the overall absolute L²-error on the electromagnetic fields ($H_x, H_y, E_z$) computed as Eq. 23 for four increasingly fine structured meshes and for polynomial orders $P_1$ to $P_4$.

$$ \mathcal{E}_{L^2}(t) = \|\mathbf{E} - \mathbf{E}_h\|_{L^2} + \|\mathbf{H} - \mathbf{H}_h\|_{L^2} \qquad (23) $$

where

$$ \|f\|_{L^p} = \left( \int_{\Omega} |f(x, y)|^p dxdy \right)^{1/p} $$

We notice in Fig. 20 that the L²-error for the sinusoidal permittivity case with the manufactured solution (in red) follows the linear trend of the L²-error for the constant permittivity
case (in blue), with some oscillations around this trend. We also observe as expected that the $L^2$-error variation rate decreases while refining the mesh and increasing the polynomial order in both cases.

In order to validate more accurately the DGTD scheme, we represent in Fig. 21 the convergence diagrams obtained for the two considered case. The error is computed as the maximum overall $L^2$-error (23) measured during the time $t_{max} = 30$ periods and is represented as a function of the minimum characteristic size $h_m$ (taken as the triangle height) of the elements for the $m^{th}$ mesh. Since we considered a structured mesh, the characteristic size is constant throughout the mesh.

We can now compute the numerical rate of convergence of the scheme for two meshes of characteristic size $h_m$ and $h_{m+1}$ as:

$$r_{\frac{m+1}{m}} = \frac{\log\left(\frac{\max_{t \in [0, t_{\max}]} ||\mathbf{E} - \mathbf{E}_{h_m}||_{L^2} + ||\mathbf{H} - \mathbf{H}_{h_m}||_{L^2}}{\max_{t \in [0, t_{\max}]} ||\mathbf{E} - \mathbf{E}_{h_{m+1}}||_{L^2} + ||\mathbf{H} - \mathbf{H}_{h_{m+1}}||_{L^2}}\right)}{\log\left(\frac{h_m}{h_{m+1}}\right)} \quad (24)$$

The results stored in Tab. 6 show that the numerical convergence rates obtained for a constant and a time-varying permittivity in the cavity are very close, except for an anomaly in $\mathbb{P}_2$. We suppose that the latter can be explained by the oscillating nature of the $L^2$-error as a function of time when the permittivity is time-varying which can lead to an extreme value in the maximum over the period $t_{max}$.

For our DGTD scheme, we used centered fluxes with an order 2 Leap-Frog time-stepping scheme. In this case, the theoretical asymptotic error is in $O(\Delta t^2 + h^p)$ for polynomial order $\mathbb{P}_p$. (see [21] for details and proof). The numerical rates are coherent with the theoretical orders although we notice that in $\mathbb{P}_1$, the rate is higher than expected (superconvergence). Since we use a second order time-stepping scheme, we cannot achieve numerical rates greater than 2 even in $\mathbb{P}_3$ and $\mathbb{P}_4$, as expected. We also notice that the convergence rates decrease when we refine the mesh too much. This can be explained by the fact that the solution converges too quickly, hence reaching the numerical precision of the machine.

<table><caption>**Table 6:** Numerical convergence rates for the rectangular resonant cavity, in the case of a constant unit permittivity $\epsilon_r = 1$ and for a sinusoidal time-dependent permittivity $\epsilon_r = f(t)$.</caption><thead><tr><th rowspan="2"></th><th colspan="2">M<sub>1</sub>(h = 0.071)</th><th colspan="2">M<sub>2</sub>(h = 0.035)</th><th colspan="2">M<sub>3</sub>(h = 0.018)</th><th colspan="2">M<sub>4</sub>(h = 0.009)</th></tr><tr><th>ε<sub>r</sub> = 1</th><th>ε<sub>r</sub> = f(t)</th><th>ε<sub>r</sub> = 1</th><th>ε<sub>r</sub> = f(t)</th><th>ε<sub>r</sub> = 1</th><th>ε<sub>r</sub> = f(t)</th><th>ε<sub>r</sub> = 1</th><th>ε<sub>r</sub> = f(t)</th></tr></thead><tbody><tr><th>P<sub>1</sub></th><td>—</td><td>—</td><td>1.88</td><td>1.85</td><td>1.88</td><td>1.88</td><td>1.71</td><td>1.72</td></tr><tr><th>P<sub>2</sub></th><td>—</td><td>—</td><td>2.07</td><td>1.91</td><td>2.00</td><td>1.59</td><td>1.93</td><td>1.28</td></tr><tr><th>P<sub>3</sub></th><td>—</td><td>—</td><td>1.99</td><td>1.99</td><td>1.94</td><td>1.94</td><td>1.79</td><td>1.79</td></tr><tr><th>P<sub>4</sub></th><td>—</td><td>—</td><td>1.94</td><td>1.94</td><td>1.79</td><td>1.78</td><td>1.24</td><td>1.24</td></tr></tbody></table>

### 3.1.2 Comparison with analytical solutions

In this section, we use another method to validate numerically the DGTD scheme for time-varying media by deriving analytical solutions for the time-modulated resonant cavity. The goal is to derive for this case a more physically meaningful solution, compared to the artificially manufactured solution described in the previous section. Indeed, we seek to derive original analytical solutions of the Maxwell equations in time-varying media without introducing an artificial source
term in the equations.

In the following, we consider two cases depending on the permittivity modulation profile $\epsilon_r(t)$, taken as a linear function and a quadratic function respectively.

**Linear permittivity modulation** Let's suppose that the interior medium of the cavity is subject to a linear permittivity variation of initial value $\epsilon_{r0}$ and variation rate $r_m$:

$$ \epsilon_r(t) = \epsilon_{r0} + r_m \cdot t \qquad (25) $$

Following the derivation detailed in the Appendices A.1 and A.2, we show that the problem has an analytical solution defined for a strictly positive or negative variation rate $r_m \le 0$:

$$ \left\{ \begin{aligned} H_x &= [c_1 J_0(u(t)) + c_2 Y_0(u(t))] \sin(k_x x) \cos(k_y y) \\ H_y &= -\frac{k_x}{k_y} [c_1 J_0(u(t)) + c_2 Y_0(u(t))] \cos(k_x x) \sin(k_y y) \\ E_z &= \frac{|\mathbf{k}|}{k_y \sqrt{\epsilon_r(t)}} [c_1 J_1(u(t)) + c_2 Y_1(u(t))] \sin(k_x x) \sin(k_y y) \end{aligned} \right. \qquad (26) $$

where

$$ u(t) = \frac{2|\mathbf{k}|}{r_m} \sqrt{\epsilon_r(t)} $$

and $J_n$ and $Y_n$ are the $n^{th}$ order Bessel's functions of the first and second kind respectively. The constants $c_1$ and $c_2$ are defined in Appendix A.2. Note that the spatial part of the solution is identical to that of the classical stationary modes (22) obtained when the permittivity is constant.

In Fig. 22, we show the temporal profile of the simulated and analytical electromagnetic fields in the cases where the permittivity of the cavity medium increases linearly and decreases linearly.

The curves show a very good agreement between the analytical and numerical solution in both cases.

**Quadratic permittivity modulation** Let's now suppose that the interior medium of the cavity is subject to a quadratic permittivity variation, such that:

$$ \epsilon_r(t) = at^2 + bt + c $$

Following the derivation in Appendix A.3, we show that the problem has an analytical solution defined for $\Delta = b^2 - 4ac \neq 0$ by:

$$ \left\{ \begin{aligned} H_x &= [c_1 P_{\lambda}(u(t)) + c_2 Q_{\lambda}(u(t))] \sin(k_x x) \cos(k_y y) \\ H_y &= -\frac{k_x}{k_y} [c_1 P_{\lambda}(u(t)) + c_2 Q_{\lambda}(u(t))] \cos(k_x x) \sin(k_y y) \\ E_z &= -\frac{2a\lambda}{k_y \sqrt{\Delta(1 - u^2(t))}} [c_1 (P_{\lambda-1}(u(t)) - uP_{\lambda}(u(t))) + c_2 (Q_{\lambda-1}(u(t)) - uQ_{\lambda}(u(t)))] \sin(k_x x) \sin(k_y y) \end{aligned} \right. $$

where
$$u(t) = \frac{1}{\sqrt{\Delta}}(2at + b)$$

and $P_λ$ and $Q_λ$ are the Legendre functions of the first and second kind respectively of non-integer degree

$$\lambda = \frac{1}{2} \left( -1 \pm \sqrt{1 - \frac{4|\mathbf{k}|^2}{a}} \right) \in \mathbb{C} \text{ for } a \neq 0$$

The constants $c_1$ and $c_2$ are defined in Appendix A.2.

In particular, we consider that the permittivity time-dependence describes a parabola which varies from $\epsilon_1$ at $t = 0$ to $\epsilon_2$ at $t = \tau$ and goes back to $\epsilon_1$ at $t = 2\tau$, which translates to:

$$\epsilon_r(t) = \frac{\epsilon_1 - \epsilon_2}{\tau^2} (t - \tau)^2 + \epsilon_2$$

In Fig. 23, we show the numerical and analytical electromagnetic fields for a parabolic permittivity variation from $\epsilon_1 = 1$ to $\epsilon_2 = 5$ at $\tau = 10$ on the left and from $\epsilon_1 = 1$ to $\epsilon_2 = 10$ at $\tau = 10$ on the right.

In both cases, the curves show very good agreement between the numerical and analytical electromagnetic fields.

# 4 Conclusion & Perspectives

In this part, we successfully derived and validated the two-dimensional implementation of a numerical scheme based on the Discontinuous Galerkin Time-Domain (DGTD) method, extended to wave propagation in time-varying media. This is an important step towards modeling the optical response of active metasurface designs by allowing us to simulate wave propagation through structures composed of time-modulated nano-resonators.

In order to design realistic active metasurfaces, there remains a few important steps in the implementation of the code. First, it is important to implement periodic boundary conditions in order to model infinite structures along one or two axes. This will allow us to model periodic structures such as gratings by considering a single unit cell. In the derivation of the scheme, we considered that the materials were non-dispersive for simplicity, but adding a dispersive model would allow for a more realistic behavior of the materials. Extending the code to 3D is also an important step, and its parallelization would reduce the consequent computational time. Lastly, coupling the DGTD code with a global optimization algorithm is crucial to finding optimal active metasurface designs.
<figure><img src="image_23.png" /><figcaption>**Figure 19:** Numerical solution for the TM<sub>z</sub> (1, 1) mode inside a two-dimensional rectangular cavity filled with medium of time-invariant and sinusoidally time-varying permittivity. (a): Spatial part of the numerical fields in both cases. (b): Temporal evolution of the numerical fields measured near the center of the cavity in the case of a time-invariant permittivity (left) and a sinusoidally varying permittivity with modulation depth δ = 0.5 and modulation frequency f<sub>m</sub> = 0.05 GHz (right). Note that we have ω<sub>m</sub> = 2πf<sub>m</sub>.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Figure 20:** Time evolution of the overall absolute $L^2$-error for the electromagnetic fields ($H_x, H_y, E_z$) between the numerical solution and the expected analytical solution as a function of the mesh size and the polynomial order. The blue curve represents the classical constant permittivity case and the red curve represents the case of a sinusoidal permittivity modulation ($\delta = 0.1, f_m = 0.05$ GHz) with the analytical manufactured solution.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 21:** Convergence diagrams obtained for the case of a resonant cavity filled with a medium of constant permittivity (a) and a medium of sinusoidally modulated permittivity where an additional electrical current has been introduced (MMS) (b).</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Figure 22:** Temporal evolution of the numerical and analytical fields measured near the center of the cavity in the cases of a linear increase (**left**) and a linear decrease (**right**) in the interior medium's permittivity.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Figure 23:** Temporal evolution of the numerical and analytical fields measured near the center of the cavity in the cases of two quadratic permittivity variations (**left**) and (**right**) in the interior medium's permittivity.</figcaption></figure>
# A Derivation of analytical solutions for a stationary wave in a resonant cavity filled with a time-varying medium

## A.1 Derivation for a general permittivity modulation profile $\epsilon_r(t)$

We consider a 2D rectangular resonant cavity filled with a medium of time-varying permittivity $\epsilon_r(t)$ and unit permeability ($\mu_r = 1$). In the case where the permittivity is time-invariant $\epsilon_r(t) = \epsilon_{r0}$, the analytical solution of the problem is well-known and corresponds to stationary modes supported by the cavity. In the Transverse Magnetic (TMz) mode and in the normalized unit system, the electric and magnetic fields take the form:

$$ 
\left\{ 
\begin{aligned}
H_x &= -\frac{k_y}{\omega_s} \sin(k_x x) \cos(k_y y) \sin(\omega_s t) \\
H_y &= \frac{k_x}{\omega_s} \cos(k_x x) \sin(k_y y) \sin(\omega_s t) \\
E_z &= \sin(k_x x) \sin(k_y y) \cos(\omega_s t)
\end{aligned}
\right. 
\qquad (27) 
$$

where the components of the wavevector $\mathbf{k}$ are $k_x = \frac{n_x \pi}{a_x}$ and $k_y = \frac{n_y \pi}{a_y}$ with $(n_x, n_y)$ the integer mode numbers and $(a_x, a_y)$ the dimensions of the cavity. For a given stationary mode $(n_x, n_y)$, the frequency of the wave is computed from the wavevector as $\omega_s = c_r |\mathbf{k}| = c_r \sqrt{k_x^2 + k_y^2} = \pi \sqrt{\left(\frac{n_x}{a_x}\right)^2 + \left(\frac{n_y}{a_y}\right)^2}$, where $c_r$ is the speed of the light in the medium of permittivity $\epsilon_{r0}$ in normalized unit.

Now, we seek to derive analytical solutions in the case where the medium inside the resonant cavity is subject to a permittivity modulation following the profile $\epsilon_r(t)$. In this case, we expect the spatial profiles of the fields to be unchanged compared to the constant permittivity case. However, the temporal profiles of the fields are no longer sinusoidal and depend on the permittivity modulation profile $\epsilon_r(t)$.

Hence, we suppose that the analytical fields are described by:

$$ 
\left\{ 
\begin{aligned}
H_x &= \sin(k_x x) \cos(k_y y) \phi_{H_x}(t) \\
H_y &= \cos(k_x x) \sin(k_y y) \phi_{H_y}(t) \\
E_z &= \sin(k_x x) \sin(k_y y) \phi_{E_z}(t)
\end{aligned}
\right. 
\qquad (28) 
$$

where we introduced the unknown temporal profiles for each field component $\phi_{H_x}(t)$, $\phi_{H_y}(t)$, $\phi_{E_z}(t)$.

Injecting the new field expressions in Maxwell's equations for the 2D TMz mode in time-varying media, we obtain for the first Maxwell equation:

$$ 
\begin{aligned}
\frac{\partial H_x}{\partial t} = -\frac{\partial E_z}{\partial y} &\iff \phi'_{H_x}(t) \sin(k_x x) \cos(k_y y) = -\phi_{E_z}(t) k_y \sin(k_x x) \cos(k_y y) \\
&\iff \boxed{\phi'_{H_x}(t) = -k_y \phi_{E_z}(t)}
\end{aligned} 
$$
The second Maxwell equation gives:

$$ \begin{aligned} \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x} &\iff \phi'_{H_y}(t) \cos(k_x x) \sin(k_y y) = \phi_{E_z}(t) k_x \cos(k_x x) \sin(k_y y) \\ &\iff \boxed{\phi'_{H_y}(t) = k_x \phi_{E_z}(t)} \end{aligned} $$

Finally, the third Maxwell equation gives:

$$ \begin{aligned} \frac{\partial(\epsilon_r(t)E_z)}{\partial t} = \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y} &\iff \epsilon'_r(t)\phi_{E_z}(t) \sin(k_x x) \sin(k_y y) + \epsilon_r(t)\phi'_{E_z}(t) \sin(k_x x) \sin(k_y y) \\ &= -\phi_{H_y}(t)k_x \sin(k_x x) \sin(k_y y) + \phi_{H_x}(t)k_y \sin(k_x x) \sin(k_y y) \\ &\iff \boxed{\epsilon_r(t)\phi'_{E_z}(t) + \epsilon'_r(t)\phi_{E_z}(t) = k_y\phi_{H_x}(t) - k_x\phi_{H_y}(t)} \end{aligned} $$

However, from the two first relations, we obtain:

$$ \begin{aligned} \phi'_{H_y}(t) = -\frac{k_x}{k_y}\phi'_{H_x}(t) &\iff \phi_{H_y}(t) = -\frac{k_x}{k_y}\phi_{H_x}(t) + C \\ &\iff \phi_{H_y}(t) = -\frac{k_x}{k_y}\phi_{H_x}(t) \end{aligned} $$

since we assume that $\phi_{H_x}(t=0) = \phi_{H_y}(t=0) = 0$ as in the analytical solution for the constant permittivity case. Replacing $\phi_{E_z}$ and $\phi_{H_y}$ by their expressions as a function of $\phi_{H_x}$, we obtain an ordinary differential equation of second order where $\phi_{H_x}(t)$ is the unknown:

$$ \boxed{\epsilon_r(t)\phi''_{H_x}(t) + \epsilon'_r(t)\phi'_{H_x}(t) + |\mathbf{k}|^2\phi_{H_x}(t) = 0} \qquad (29) $$

This equation can also be expressed equivalently as a function of the unknown $\phi_{H_y}$:

$$ \epsilon_r(t)\phi''_{H_y}(t) + \epsilon'_r(t)\phi'_{H_y}(t) + |\mathbf{k}|^2\phi_{H_y}(t) = 0 $$

The temporal profiles for the remaining fields can be computed from the relations:

$$ \phi_{H_y}(t) = -\frac{k_x}{k_y}\phi_{H_x}(t) $$

$$ \phi_{E_z}(t) = -\frac{1}{k_y}\phi'_{H_x}(t) $$

## A.2 Linear permittivity modulation

In this section, we seek to derive an analytical solution for the problem of a rectangular resonant cavity which interior medium is subjected to a linear permittivity variation, such that:

$$ \epsilon_r(t) = \epsilon_{r0} + r_m \cdot t $$

where $\epsilon_{r0}$ is the initial value of the permittivity and $r_m$ is the rate of permittivity modulation. Replacing the expression of $\epsilon_r(t)$ in Eq. 29 yields the following order 2 ODE:

$$ (\epsilon_{r0} + r_m t)\phi''(t) + r_m \phi'(t) + |\mathbf{k}|^2 \phi(t) = 0 \qquad (30) $$
where the unknown $\phi(t)$ denotes the temporal profile of the $H_x$ field component.

In order to reformulate the equation into a known form, we introduce the following change of variable:

$$t \to u(t) = \frac{2|\mathbf{k}|}{r_m} \sqrt{\epsilon_r(t)} = \frac{2|\mathbf{k}|}{r_m} \sqrt{\epsilon_{r0} + r_m t}$$

Using the chain rule, we can express the derivatives of the unknown $\phi(t)$ with respect to $t$ as a function of the derivatives of $\phi(u(t))$ with respect to the new variable $u(t)$. We get for the first derivative:

$$\begin{aligned} \frac{d\phi}{dt} &= \frac{d\phi}{du} \cdot \frac{du}{dt} \\ &= \frac{d\phi}{du} \frac{|\mathbf{k}|}{\sqrt{\epsilon_r(t)}} \end{aligned}$$

For the second derivative, we have:

$$\begin{aligned} \frac{d^2\phi}{dt^2} &= \frac{d^2\phi}{du^2} \left(\frac{du}{dt}\right)^2 + \frac{d\phi}{du} \cdot \frac{d^2u}{dt^2} \\ &= \frac{d^2\phi}{du^2} \frac{|\mathbf{k}|^2}{\epsilon_r(t)} - \frac{d\phi}{du} \frac{r_m|\mathbf{k}|}{2[\epsilon_r(t)]^{3/2}} \end{aligned}$$

Replacing these expressions in Eq. 30 yields:

$$\begin{align*} & \epsilon_r(t) \left( \frac{d^2\phi}{du^2} \frac{|\mathbf{k}|^2}{\epsilon_r(t)} - \frac{d\phi}{du} \frac{r_m|\mathbf{k}|}{2[\epsilon_r(t)]^{3/2}} \right) + r_m \left( \frac{d\phi}{du} \frac{|\mathbf{k}|}{\sqrt{\epsilon_r(t)}} \right) + |\mathbf{k}|^2 \phi(u) = 0 \\ \iff & |\mathbf{k}|^2 \frac{d^2\phi}{du^2} + \frac{r_m|\mathbf{k}|}{2\sqrt{\epsilon_r(t)}} \frac{d\phi}{du} + |\mathbf{k}|^2 \phi(u) = 0 \\ \iff & \frac{d^2\phi}{du^2} + \frac{1}{u} \frac{d\phi}{du} + \phi(u) = 0 \\ \iff & u^2 \frac{d^2\phi}{du^2} + u \frac{d\phi}{du} + u^2 \phi(u) = 0 \end{align*}$$

The last expression takes the form of Bessel's differential equation of order 0. Hence, the general solution of this equation is a linear combination of the 0th order Bessel's functions of the first $J_0(u)$ and second kind $Y_0(u)$.

Therefore the analytical expression for the temporal profile of the $H_x$ field component is:

$$\boxed{\phi_{H_x}(t) = c_1 J_0(u(t)) + c_2 Y_0(u(t))} \qquad (31)$$

The constants $c_1$ and $c_2$ are to be determined from the initial conditions $\phi_{H_x}(t=0)$ and $\phi'_{H_x}(t=0)$. Assuming a null magnetic field and a unit electric field in their temporal profile at $t=0$, the initial conditions on $\phi_{H_x}$ are:

$$\begin{aligned} \phi_{H_x}(t=0) &= 0 \\ \phi'_{H_x}(t=0) &= -k_y \phi_{E_z}(t=0) = -k_y \end{aligned}$$

The first condition gives:

$$c_2 = -c_1 \frac{J_0(u_0)}{Y_0(u_0)}$$
where we pose $u_0 = u(t = 0)$.

From the second condition, we get:

$$c_1 J_1(u_0) + c_2 Y_1(u_0) = \frac{k_y}{|\mathbf{k}|} \sqrt{\epsilon_{r0}}$$

Combining the two relations yields the explicit expressions for the two coefficients:

$$c_1 = \frac{k_y \sqrt{\epsilon_{r0}}}{|\mathbf{k}|} \frac{Y_0(u_0)}{J_1(u_0)Y_0(u_0) - J_0(u_0)Y_1(u_0)} \\ c_2 = - \frac{k_y \sqrt{\epsilon_{r0}}}{|\mathbf{k}|} \frac{J_0(u_0)}{J_1(u_0)Y_0(u_0) - J_0(u_0)Y_1(u_0)}$$

The analytical expressions for the temporal profiles of the remaining field components can be easily obtained from $\phi_{H_x}(t)$ as:

$$\phi_{H_y}(t) = -\frac{k_x}{k_y} [c_1 J_0(u(t)) + c_2 Y_0(u(t))] \quad (32)$$

$$\phi_{E_z}(t) = \frac{|\mathbf{k}|}{k_y \sqrt{\epsilon_r(t)}} [c_1 J_1(u(t)) + c_2 Y_1(u(t))] \quad (33)$$

In the case of a **positive permittivity variation** ($r_m > 0$), the argument of the Bessel functions is also positive ($u(t) > 0$). However, for real values of the parameter $\nu$ and for a positive argument $z$, the Bessel functions $J_\nu(z)$ and $Y_\nu(z)$ are real-valued (ref. bib. ?), hence the analytical expressions for the 3 field components are also real-valued.

However, if we now consider a **negative permittivity variation** ($r_m < 0$), leading to a negative argument ($u(t) < 0$), not all the Bessel functions are real-valued. Indeed, we have the following relations:

$$\text{For } z \in \mathbb{R}, \quad \begin{cases} J_0(-z) = J_0(z) \in \mathbb{R} \\ J_1(-z) = -J_1(z) \in \mathbb{R} \end{cases} \quad \text{For } z > 0, \quad \begin{cases} Y_0(-z) = Y_0(z) + 2iJ_0(z) \in \mathbb{C} \\ Y_1(-z) = -Y_1(z) - 2iJ_1(z) \in \mathbb{C} \end{cases}$$

Defining the constant $c_0 = J_1(u_0)Y_0(u_0) - J_0(u_0)Y_1(u_0)$ and replacing the Bessel functions by the above expressions for a negative argument $u_0 = -|u_0| < 0$, we obtain:

$$\begin{aligned} c_0 &= J_1(-|u_0|)Y_0(-|u_0|) - J_0(-|u_0|)Y_1(-|u_0|) \\ &= -J_1(|u_0|) [Y_0(|u_0|) + 2iJ_0(|u_0|)] - J_0(|u_0|) [-Y_1(|u_0|) - 2iJ_1(|u_0|)] \\ &= -[J_1(|u_0|)Y_0(|u_0|) - J_0(|u_0|)Y_1(|u_0|)] \in \mathbb{R} \end{aligned}$$

For the temporal profile of $H_x$, we have:

$$\begin{aligned} \phi_{H_x}(t) &= \frac{k_y/|\mathbf{k}|}{c_0} \sqrt{\epsilon_{r0}} [Y_0(-|u_0|)J_0(-|u(t)|) - J_0(-|u_0|)Y_0(-|u(t)|)] \\ &= \frac{k_y/|\mathbf{k}|}{c_0} \sqrt{\epsilon_{r0}} [(Y_0(|u_0|) + 2iJ_0(|u_0|))J_0(|u(t)|) - J_0(|u_0|)(Y_0(|u(t)|) + 2iJ_0(|u(t)|))] \\ &= \frac{k_y/|\mathbf{k}|}{c_0} \sqrt{\epsilon_{r0}} [Y_0(|u_0|)J_0(|u(t)|) - J_0(|u_0|)Y_0(|u(t)|)] \in \mathbb{R} \end{aligned}$$

Inria
For the temporal profile of $E_z$, we have:

$$ 
\begin{aligned}
\phi_{E_z}(t) &= \frac{\sqrt{\epsilon_{r0}/\epsilon_r(t)}}{c_0} [Y_0(-|u_0|)J_1(-|u(t)|) - J_0(-|u_0|)Y_1(-|u(t)|)] \\
&= \frac{\sqrt{\epsilon_{r0}/\epsilon_r(t)}}{c_0} [(Y_0(|u_0|) + 2iJ_0(|u_0|))(-J_1(|u(t)|)) - J_0(|u_0|)(-Y_1(|u(t)|) - 2iJ_1(|u(t)|))] \\
&= -\frac{\sqrt{\epsilon_{r0}/\epsilon_r(t)}}{c_0} [Y_0(|u_0|)J_1(|u(t)|) - J_0(|u_0|)Y_1(|u(t)|)]
\end{aligned}
$$

In the general case, for a strictly positive or negative variation rate ($r_m \le 0$), we pose the new variable:

$$ u(t) = \left| \frac{2|\mathbf{k}|}{r_m} \sqrt{\epsilon_r(t)} \right| = \left| \frac{2|\mathbf{k}|}{|r_m|} \sqrt{\epsilon_r(t)} \right| $$

The temporal profiles of the field components take the more general expressions:

$$ \phi_{H_x}(t) = \operatorname{sgn}(r_m) \times [c_1 J_0(u(t)) + c_2 Y_0(u(t))] \quad (34) $$

$$ \phi_{H_y}(t) = -\operatorname{sgn}(r_m) \frac{k_x}{k_y} \times [c_1 J_0(u(t)) + c_2 Y_0(u(t))] \quad (35) $$

$$ \phi_{E_z}(t) = \frac{|\mathbf{k}|}{k_y \sqrt{\epsilon_r(t)}} [c_1 J_1(u(t)) + c_2 Y_1(u(t))] \quad (36) $$

### A.3 Quadratic permittivity modulation

In this section, we seek to derive an analytical solution in the case where the interior medium of the cavity is subjected to a quadratic permittivity variation, such that:

$$ \epsilon_r(t) = at^2 + bt + c, \quad a \neq 0 $$

Replacing the expression of $\epsilon_r(t)$ in Eq. 29 yields:

$$ (at^2 + bt + c)\phi''(t) + (2at + b)\phi'(t) + |\mathbf{k}|^2\phi(t) = 0 \quad (37) $$

where the unknown $\phi(t)$ denotes the temporal profile of the $H_x$ field component.

Let's pose the new variable:

$$ t \rightarrow u(t) = \frac{1}{\sqrt{b^2 - 4ac}} (2at + b) = \frac{1}{\sqrt{\Delta}} (2at + b) $$

Using the chain rule, we can express the derivatives of the unknown $\phi(t)$ with respect to $t$ as a function of the derivatives of $\phi(u(t))$ with respect to the new variable $u(t)$. We get for the first derivative:

$$ 
\begin{aligned}
\frac{d\phi}{dt} &= \frac{d\phi}{du} \cdot \frac{du}{dt} \\
&= \frac{d\phi}{du} \frac{2a}{\sqrt{\Delta}}
\end{aligned}
$$
For the second derivative, we have:

$$ \begin{aligned} \frac{d^2\phi}{dt^2} &= \frac{d^2\phi}{du^2} \left(\frac{du}{dt}\right)^2 + \frac{d\phi}{du} \cdot \frac{d^2u}{dt^2} \\ &= \frac{d^2\phi}{du^2} \frac{4a^2}{\Delta} \end{aligned} $$

Replacing these expressions in Eq. 37 yields:

$$ \begin{aligned} & \underbrace{(at^2 + bt + c)}_{a \left(t + \frac{b}{2a}\right)^2 - \frac{\Delta}{4a}} \frac{4a^2 d^2\phi}{\Delta du^2} + \underbrace{(2at + b) \frac{2a}{\sqrt{\Delta}} \frac{d\phi}{du}}_{2au} + |\mathbf{k}|^2 \phi(u) = 0 \\ \iff & \underbrace{\left(\frac{4a^2}{\Delta} \left(t + \frac{b}{2a}\right)^2 - 1\right)}_{u^2} \frac{d^2\phi}{du^2} + 2u \frac{d\phi}{du} + \frac{|\mathbf{k}|^2}{a} \phi(u) = 0 \\ \iff & (1 - u^2) \frac{d^2\phi}{du^2} - 2u \frac{d\phi}{du} \underbrace{- \frac{|\mathbf{k}|^2}{a}}_{\lambda(\lambda + 1)} \phi(u) = 0 \end{aligned} $$

This last expression corresponds to the expression of the Legendre differential equation of non-integer degree $\lambda$ such that:

$$ \lambda(\lambda + 1) = -\frac{|\mathbf{k}|^2}{a} \\ \iff \lambda = \frac{1}{2} \left( -1 \pm \sqrt{1 - \frac{4|\mathbf{k}|^2}{a}} \right) \in \mathbb{C} \text{ for } a \neq 0 $$

We can choose either value for the degree $\lambda$ of the Legendre equation which is a complex number for a general value of $a \neq 0$.

The Legendre equation has two linearly independent solution which are the Legendre functions of the first and second kind of non-integer degree $\lambda$: $P_{\lambda}(u)$ and $Q_{\lambda}(u)$ respectively. Hence the general solution for Eq. 37 is:

$$ \phi_{H_x} = c_1 P_{\lambda}(u(t)) + c_2 Q_{\lambda}(u(t)) \quad (38) $$

where the constants $c_1$ and $c_2$ can be computed from the initial conditions $\phi_{H_x}(t = 0)$ and $\phi'_{H_x}(t = 0)$.

We still assume that a null magnetic field and a unit electric field at $t = 0$, hence:

$$ \begin{aligned} \phi_{H_x}(t = 0) &= 0 \\ \phi'_{H_x}(t = 0) &= -k_y \phi_{E_z}(t = 0) = -k_y \end{aligned} $$

From the first relation, we deduce:

$$ c_2 = -c_1 \frac{P_{\lambda}(u_0)}{Q_{\lambda}(u_0)} $$
where $u_0 = u(t = 0) = \frac{b}{\sqrt{\Delta}}$.

Differentiating Eq. 38 with respect to time yields:

$$ \frac{d\phi_{H_x}}{dt} = \frac{du}{dt} \left[ c_1 \frac{d}{du} P_{\lambda}(u) + c_2 \underbrace{\frac{d}{du} Q_{\lambda}(u)}_{\frac{2a}{\Delta}} \right] $$

where the first derivative of the Legendre functions can be expressed from the recurrence formulae (see [23] for the derivation):

$$ \frac{d}{du} P_{\lambda}(u) = \frac{\lambda(P_{\lambda-1}(u) - uP_{\lambda}(u))}{1 - u^2} $$

$$ \frac{d}{du} Q_{\lambda}(u) = \frac{\lambda(Q_{\lambda-1}(u) - uQ_{\lambda}(u))}{1 - u^2} $$

Replacing these formulae in the previous expression and gives at $t = 0$:

$$ \phi'_{H_x}(t) = \frac{2a}{\sqrt{\Delta}} \frac{\lambda}{1 - u_0^2} [c_1(P_{\lambda-1}(u_0) - u_0P_{\lambda}(u_0)) + c_2(Q_{\lambda-1}(u_0) - u_0Q_{\lambda}(u_0))] = -k_y $$

Finally, by replacing $c_2$ by the first relation, we obtain the expressions for the constants:

$$ c_1 = \frac{k_y \sqrt{\Delta} (1 - u_0^2)}{2a\lambda} \frac{Q_{\lambda}(u_0)}{P_{\lambda}(u_0)Q_{\lambda-1}(u_0) - Q_{\lambda}(u_0)P_{\lambda-1}(u_0)} $$

$$ c_2 = - \frac{k_y \sqrt{\Delta} (1 - u_0^2)}{2a\lambda} \frac{P_{\lambda}(u_0)}{P_{\lambda}(u_0)Q_{\lambda-1}(u_0) - Q_{\lambda}(u_0)P_{\lambda-1}(u_0)} $$

The expressions for the temporal profiles of the remaining field components can be easily obtained from $\phi_{H_x}(t)$ as:

$$ \phi_{H_y} = -\frac{k_x}{k_y} [c_1 P_{\lambda}(u(t)) + c_2 Q_{\lambda}(u(t))] \quad (39) $$

$$ \phi_{E_z} = -\frac{2a\lambda}{k_y \sqrt{\Delta}(1 - u^2(t))} \left[ c_1 \left( P_{\lambda-1}(u(t)) - uP_{\lambda}(u(t)) \right) + c_2 \left( Q_{\lambda-1}(u(t)) - uQ_{\lambda}(u(t)) \right) \right] \quad (40) $$
References

[1] Jean-Pierre Bérenger. "Perfectly Matched Layer (PML) for Computational Electromagnetics". en. In: Synthesis Lectures on Computational Electromagnetics. Cham: Springer International Publishing, 2007. ISBN: 978-3-031-00568-8 978-3-031-01696-7. DOI: 10.1007/978-3-031-01696-7.

[2] Mahmoud M. R. Elsawy et al. "Global optimization of metasurface designs using statistical learning methods". en. In: *Scientific Reports* 9.1 (Nov. 2019), p. 17918. ISSN: 2045-2322. DOI: 10.1038/s41598-019-53878-9.

[3] Elena Mikheeva et al. "Space and Time Modulations of Light with Metasurfaces: Recent Progress and Future Prospects". In: *ACS Photonics* 9.5 (2022), pp. 1458-1482. DOI: 10.1021/acsphotonics.1c01833.

[4] Andrei Komar et al. "Dynamic Beam Switching by Liquid Crystal Tunable Dielectric Metasurfaces". en. In: *ACS Photonics* 5.5 (May 2018), pp. 1742-1748. ISSN: 2330-4022, 2330-4022. DOI: 10.1021/acsphotonics.7b01343.

[5] Tomer Lewi et al. "Ultrawide Thermo-optic Tuning of PbTe Meta-Atoms". In: *Nano Letters* 17.6 (June 2017), pp. 3940-3945. ISSN: 1530-6984. DOI: 10.1021/acs.nanolett.7b01529.

[6] Ileana-Cristina Benea-Chelmus et al. "Electro-optic spatial light modulator from an engineered organic layer". en. In: *Nature Communications* 12.1 (Oct. 2021), p. 5928. ISSN: 2041-1723. DOI: 10.1038/s41467-021-26035-y.

[7] Yao-Wei Huang et al. "Gate-Tunable Conducting Oxide Metasurfaces". In: *Nano Letters* 16.9 (Sept. 2016), pp. 5319-5325. ISSN: 1530-6984. DOI: 10.1021/acs.nanolett.6b00555.

[8] Arseniy I. Kuznetsov et al. "Roadmap for Optical Metasurfaces". en. In: *ACS Photonics* (Feb. 2024), acsphotonics.3c00457. ISSN: 2330-4022, 2330-4022. DOI: 10.1021/acsphotonics.3c00457.

[9] Prachi Thureja et al. "Toward a universal metasurface for optical imaging, communication, and computation". en. In: *Nanophotonics* 11.17 (Sept. 2022), pp. 3745-3768. ISSN: 2192-8614. DOI: 10.1515/nanoph-2022-0155.

[10] Amr Shaltout, Alexander Kildishev, and Vladimir Shalaev. "Time-varying metasurfaces and Lorentz non-reciprocity". en. In: *Optical Materials Express* 5.11 (Nov. 2015), pp. 2459-2467. ISSN: 2159-3930. DOI: 10.1364/OME.5.002459.

[11] Raana Sabri and Hossein Mosallaei. "High-quality-factor space-time metasurface for free-space power isolation at near-infrared regime". en. In: *Advanced Photonics Nexus* 2.06 (Dec. 2023). ISSN: 2791-1519. DOI: 10.1117/1.APN.2.6.066008.

[12] Sajjad Taravati and George V. Eleftheriades. "Generalized Space-Time Periodic Diffraction Gratings: Theory and Applications". en. In: *Physical Review Applied* 12.2 (Aug. 2019). arXiv:1902.09885 [physics], p. 024026. ISSN: 2331-7019. DOI: 10.1103/PhysRevApplied.12.024026.

[13] Kane Yee. "Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media". en. In: *IEEE Transactions on Antennas and Propagation* 14.3 (May 1966), pp. 302-307. ISSN: 0018-926X, 1558-2221. DOI: 10.1109/TAP.1966.1138693.

[14] Xiaowen Liu. "The use of the FDTD method for electromagnetic analysis in the presence of time-varying media". en. PhD thesis. University of Ottawa (Canada), 2004.
[15] Mahdi Chegnizadeh, Mohammad Memarian, and Khashayar Mehrany. "Wave Propagation in Metallic Slab Waveguides Undergoing Arbitrary Temporal Variations of Permittivity". en. In: *2019 27th Iranian Conference on Electrical Engineering (ICEE)*. Yazd, Iran: IEEE, Apr. 2019, pp. 1529-1533. ISBN: 978-1-72811-508-5. DOI: 10.1109/IranianCEE.2019.8786763.

[16] Anand Kumar and Debdeep Sarkar. "FDTD Analysis of Space-time Metamaterials using Modulated TVTLs for Frequency Translation, Mixing and Non-reciprocity". en. In: *2023 53rd European Microwave Conference (EuMC)*. Berlin, Germany: IEEE, Sept. 2023, pp. 311-314. ISBN: 978-2-87487-072-9. DOI: 10.23919/EuMC58039.2023.10290288.

[17] Sandeep Inampudi et al. "Rigorous space-time coupled-wave analysis for patterned surfaces with temporal permittivity modulation [Invited]". en. In: *Optical Materials Express* 9.1 (Jan. 2019), p. 162. ISSN: 2159-3930. DOI: 10.1364/OME.9.000162.

[18] M. G. Moharam et al. "Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings". en. In: *JOSA A* 12.5 (May 1995), pp. 1068-1076. ISSN: 1520-8532. DOI: 10.1364/JOSAA.12.001068.

[19] Jonathan Viquerat. "Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method". PhD thesis.

[20] W. H. Reed and T. R. Hill. *Triangular mesh methods for the neutron transport equation*. English. LA-UR-73-479; CONF-730414-2. Oct. 1973.

[21] Loula Fezoui et al. "Convergence and stability of a discontinuous Galerkin time-domain method for the 3D heterogeneous Maxwell equations on unstructured meshes". en. In: *ESAIM: Mathematical Modelling and Numerical Analysis* 39.6 (Nov. 2005), pp. 1149-1176. ISSN: 0764-583X, 1290-3841. DOI: 10.1051/m2an:2005049.

[22] István Faragó, Ágnes Havasi, and Robert Horváth. "Numerical solution of the Maxwell equations in time-varying media using Magnus expansion". en. In: *Central European Journal of Mathematics* 10.1 (Feb. 2012), pp. 137-149. ISSN: 1895-1074, 1644-3616. DOI: 10.2478/s11533-011-0074-3.

[23] E. T. Whittaker and G. N. Watson. "Legendre Functions". In: *A Course of Modern Analysis*. Cambridge Mathematical Library. Cambridge University Press, 1996, pp. 302-336.
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