Application of Physics-Informed Neural Networks to
Maxwell equations

Alexandre Pugin, Stéphane Lanteri, Stéphane Descombes, Mahmoud Elsawy

► To cite this version:

Alexandre Pugin, Stéphane Lanteri, Stéphane Descombes, Mahmoud Elsawy. Application of Physics-Informed Neural Networks to Maxwell equations. RR-9588, Centre Inria d'Université Côte d'Azur. 2025. hal-05073106

HAL Id: hal-05073106

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

Submitted on 20 May 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

Application of
Physics-Informed Neural
Networks to Maxwell
equations

Alexandre Pugin, Stéphane Lanteri, Stéphane Descombes, Mahmoud Elsawy

RESEARCH
REPORT

Nº 9588

May 2025

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

# Application of Physics-Informed Neural Networks to Maxwell equations

Alexandre Pugin, Stéphane Lanteri, Stéphane Descombes,
Mahmoud Elsawy

Project-Teams Atlantis

Research Report n° 9588 — May 2025 — 59 pages

**Abstract:** In the last decade, the study of electromagnetic waves has been a primary concern in many applications such as designing radar or antenna technology as well as cloaking devices and metalenses which fill the role of a lens through the use of metamaterials. The electromagnetic field follows the Maxwell equations that in most applications have no analytic solutions. Traditional ways to solve this system, which are based on numerical methods such as the finite difference or the finite element methods, are computationally expensive. This becomes a problem when a lot of simulations need to be done such as when doing numerical optimisation. Physics-informed neural networks (PINNs) have been introduced only recently in the world of Scientific Machine Learning (SciML) as a novel way to solve systems of partial differential equations. Therefore, in this study we would like to study the concept of PINNs for the Maxwell equations to generate surrogate models, trading accuracy for faster results. We pursue three objectives: (1) establish an appropriate PINN framework for the Maxwell equations ; (2) assess this PINN framework on a collection of simple cases to eventually build towards complex and realistic applications ; (3) complete the basic PINN implementation with additional techniques to improve results and test them on the previous cases. In the end, we found that our framework was effective to solve the simpler cases considered but that more difficult cases benefited greatly from the addition of techniques.

**Key-words:** physics-informed neural networks, partial differential equations, Maxwell equations, deep learning, scientific machine learning

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

2004 route des Lucioles - BP 93
06902 Sophia Antipolis Cedex
# Application des réseaux de neurones informés par la physique aux équations de Maxwell

**Résumé** : Durant la dernière décennie, l'étude des ondes électromagnétiques a été une préoccupation majeure dans de nombreux domaines. Notamment pour la conception de radars ou d'antennes ainsi que pour les technologies furtives et les métalentilles qui remplissent le rôle de lentille grâce à l'utilisation de métamatériaux. Le champ électromagnétique suit les équations de Maxwell qui, dans la plupart des applications, n'ont pas de solution analytique. La façon traditionnelle de résoudre ce système, avec la méthode des éléments finis, est coûteuse en termes de temps de calcul. Cela devient un problème lorsque de nombreuses simulations doivent être effectuées, par exemple lors d'une optimisation. Les réseaux de neurones informés par la physique (PINN) n'ont été introduits que récemment dans le monde de l'apprentissage automatique scientifique comme nouvelle façon de résoudre des systèmes d'équations aux dérivées partielles. Par conséquent, dans cette étude, nous aimerions utiliser les PINN avec les équations de Maxwell pour générer des modèles de substitution, quitte à perdre un peu de précision pour avoir des résultats plus rapidement. Nous avons trois objectifs : (1) établir un modèle PINN pour les équations de Maxwell ; (2) évaluer ce modèle sur un ensemble de cas simples pour éventuellement aller vers des applications complexes et réalistes ; (3) compléter le PINN avec des techniques supplémentaires pour améliorer les résultats et les tester sur les cas précédents. En fin de compte, nous avons constaté que les PINNs pouvaient résoudre les cas les plus simples, mais que les plus difficiles nécessitaient l'ajout de techniques appropriées.

**Mots-clés** : réseaux neuronaux informés par la physique, équations aux dérivées partielles, apprentissage profond, apprentissage automatique scientifique
# Contents

<table><tbody><tr><td>1</td><td>Mathematical Framework</td><td>6</td></tr><tr><td>1.1</td><td>List of symbols</td><td>6</td></tr><tr><td>1.2</td><td>Maxwell equations</td><td>6</td></tr><tr><td>1.3</td><td>Boundary and interface conditions</td><td>7</td></tr><tr><td>1.3.1</td><td>Boundary conditions for a perfect electric conductor</td><td>7</td></tr><tr><td>1.3.2</td><td>Absorbing boundary condition</td><td>7</td></tr><tr><td>1.3.3</td><td>Interface conditions</td><td>7</td></tr><tr><td>1.4</td><td>Transverse magnetic and transverse electric waves</td><td>8</td></tr><tr><td>1.4.1</td><td>System of equations for the propagation</td><td>8</td></tr><tr><td>1.4.2</td><td>Boundary and interface conditions</td><td>9</td></tr><tr><td>1.5</td><td>Renormalization</td><td>9</td></tr><tr><td>1.5.1</td><td>Usual adimensionalization for numerical methods</td><td>9</td></tr><tr><td>1.5.2</td><td>Further scaling for the inputs to a neural network</td><td>10</td></tr><tr><td>1.5.3</td><td>Systems for transverse waves</td><td>11</td></tr><tr><td>1.5.4</td><td>Boundary and interface conditions</td><td>11</td></tr><tr><td>1.6</td><td>Physics-Informed Neural Networks</td><td>12</td></tr><tr><td>2</td><td>Test cases</td><td>13</td></tr><tr><td>2.1</td><td>C1: Stationary wave in a metallic cavity</td><td>13</td></tr><tr><td>2.1.1</td><td>Resonant mode (1,1) in a square cavity</td><td>13</td></tr><tr><td>2.1.2</td><td>Higher-order resonant mode in a square cavity</td><td>14</td></tr><tr><td>2.2</td><td>C2: Propagation of a plane wave in a homogeneous domain</td><td>15</td></tr><tr><td>2.2.1</td><td>Incident monochromatic plane wave</td><td>17</td></tr><tr><td>2.2.2</td><td>Incident gaussian plane wave</td><td>17</td></tr><tr><td>2.3</td><td>C3: Flat interface between two dielectric materials</td><td>18</td></tr><tr><td>2.4</td><td>C4: Dielectric slab</td><td>18</td></tr><tr><td>2.5</td><td>C5: Diffraction of an incident wave by a dielectric cylinder</td><td>20</td></tr><tr><td>3</td><td>Sampling of the collocation points</td><td>22</td></tr><tr><td>3.1</td><td>Sampling methods</td><td>22</td></tr><tr><td>3.2</td><td>Sampling strategies</td><td>23</td></tr><tr><td>4</td><td>Neural networks architecture</td><td>24</td></tr><tr><td>4.1</td><td>Vanilla MLP</td><td>24</td></tr><tr><td>4.2</td><td>SIREN</td><td>24</td></tr><tr><td>4.3</td><td>Fourier embedding</td><td>24</td></tr><tr><td>4.4</td><td>Imposing periodicity</td><td>24</td></tr><tr><td>4.5</td><td>Domain decomposition for heterogeneous media</td><td>25</td></tr><tr><td>4.6</td><td>Scattering problem</td><td>25</td></tr><tr><td>5</td><td>Training dynamics</td><td>27</td></tr><tr><td>5.1</td><td>Convergence and conditioning</td><td>27</td></tr><tr><td>5.2</td><td>Energy Natural Gradients</td><td>28</td></tr></tbody></table>
<table><tr><td><b>6</b></td><td><b>Results</b></td><td><b>30</b></td></tr><tr><td>6.1</td><td>Architecture</td><td>30</td></tr><tr><td>6.2</td><td>Resampling</td><td>32</td></tr><tr><td>6.3</td><td>Optimizer</td><td>35</td></tr><tr><td>6.4</td><td>Results for each case</td><td>37</td></tr><tr><td>6.4.1</td><td>C1: Stationary wave in a metallic cavity</td><td>37</td></tr><tr><td>6.4.2</td><td>C2: Propagation of a plane wave in a homogeneous domain</td><td>40</td></tr><tr><td>6.4.3</td><td>C3: Flat interface between two dielectric materials</td><td>42</td></tr><tr><td>6.4.4</td><td>C4: Dielectric slab</td><td>46</td></tr><tr><td>6.4.5</td><td>C5: Diffraction of an incident wave by a dielectric cylinder</td><td>47</td></tr><tr><td>6.5</td><td>Discussion</td><td>52</td></tr><tr><td><b>7</b></td><td><b>Conclusions and future developments</b></td><td><b>53</b></td></tr><tr><td><b>A</b></td><td><b>Analytical solutions</b></td><td><b>55</b></td></tr><tr><td>A.1</td><td>Flat interface</td><td>56</td></tr><tr><td>A.2</td><td>Dielectric slab</td><td>57</td></tr></table>
# Introduction

Maxwell's equations govern the propagation of electromagnetic waves. Solving them is necessary to design metamaterials in radar, antenna, cloaking and metalenses technologies for example. Currently, simulations are predominantly conducted using numerical methods such as the Finite Element Method (FEM) and the Finite Difference Method (FDM). Simulations with these methods can become computationally intensive in the case of three-dimensional domains with complex geometries. This represents a strong limitation if the ultimate goal is to perform numerical optimization for example, which requires many simulations for a large set of design parameters. Hence, we would like to use neural networks as surrogate models, possibly with lower accuracy for shorter simulation times.

Physics-Informed Neural Networks (PINNs) have seen a surge in interest over the last five years after the original article from Raissi et al.[12] was published. This paradigm allows one to model the solution of systems of partial differential equations with little to no data provided. PINNs have been applied to many physical systems, such as Navier-Stokes[1] or the heat transfer equations[13]. There is also extensive research on ways to improve PINNs in all settings, no matter what partial differential equations are used. Wang et al.[21] describe many techniques and general approaches, which yield success for a variety of systems of partial differential equations they tested, though Maxwell's equations were not one of them. Sukumar et al.[16] present a way to impose the boundary conditions through hard-constraints which allows the removal of one of the criterion to optimize. Also, Dong et al.[4] show how to impose periodicity in a neural network which is very useful if the solution to be found is known to be periodic. Finally, there are also surveys attempting to list PINNs advances and to summarize the disparate and complex literature for PINNs, Toscano et al.[18] compiles a comprehensive overview of the latest advancements in PINNs and Physics-Informed Kolmogorov-Arnold Networks (PIKANs).

However, Maxwell's equations remain poorly explored with PINNs despite the incentives. There are only a few works dealing with PINNs to solve Maxwell equations. The goal of Chen et al.[3] resembles ours, as they want to use PINNs to design metamaterials. However, in contrast to us, they solve inverse problems, i.e. they have data dictating what the electromagnetic field should look like and from there they use PINNs to retrieve the effective permittivity of metamaterials. This is different from our goal, since we instead want to determine the electromagnetic field given the permittivity, focusing on the design of parametric PINNs. In this respect, Baldan et al.[2] is closer to us as they solve a parametric problem, i.e. for a set of given parameters, they predict the electromagnetic field. Here, the parameters are the positions of the wires. But unlike us, they do so in the magnetostatic setting which means the magnetic field is supposed constant in time. We want to solve Maxwell's equations in a more general setting where the magnetic field is time-dependent.

The goal of this report is to lay the groundwork to eventually apply PINNs to solve parametric electromagnetic wave propagation problems in three-dimensional domains. To do so, we will start by solving forward problems in two-dimensional spaces. First, we will delve into the mathematical framework for the equations we are trying to solve as well as the concept of PINNs with their original formulation. Second, we will describe some modifications that can be brought to PINNs to improve their performance. Then, we will present some results of simulations with PINNs.
# 1 Mathematical Framework

## 1.1 List of symbols

<table><thead><tr><th>Variable</th><th>Description</th><th>SI unit</th></tr></thead><tbody><tr><td><b>E</b></td><td>Electric field</td><td>V m<sup>-1</sup></td></tr><tr><td><b>D</b></td><td>Electric induction</td><td>C m<sup>-2</sup></td></tr><tr><td><b>H</b></td><td>Magnetic field</td><td>A m<sup>-1</sup></td></tr><tr><td><b>B</b></td><td>Magnetic induction</td><td>T</td></tr><tr><td><b>x</b></td><td>Position</td><td>m</td></tr><tr><td><b>t</b></td><td>Time</td><td>s</td></tr><tr><td><b>ρ</b></td><td>Density of electric charges</td><td>C m<sup>-3</sup></td></tr><tr><td><b>j</b></td><td>Density of electric currents</td><td>A m<sup>-2</sup></td></tr><tr><td><b>ε</b></td><td>Permittivity</td><td>F m<sup>-1</sup></td></tr><tr><td><b>μ</b></td><td>Permeability</td><td>H m<sup>-1</sup></td></tr><tr><td><b>Z</b></td><td>Impedance</td><td>Ω</td></tr><tr><td><b>c</b></td><td>Celerity</td><td>m s<sup>-1</sup></td></tr></tbody></table>

## 1.2 Maxwell equations

To be able to use PINNs with the Maxwell equations, we first have to change them so that they only use scalar quantities. To do so, we need to start from the most general formulation:

$$ \left\{ \begin{aligned} \frac{\partial \mathbf{D}}{\partial t} - \nabla \times \mathbf{H} &= \mathbf{j} \\ \frac{\partial \mathbf{B}}{\partial t} + \nabla \times \mathbf{E} &= \mathbf{0} \\ \nabla \cdot \mathbf{D} &= \rho \\ \nabla \cdot \mathbf{B} &= 0 \end{aligned} \right. \qquad (1) $$

where **E** and **D** are the electric field and its induction, **H** and **B** are the magnetic field and its induction, ρ is the density of electric charges and **j** is the electric current density. For linear and isotropic materials, the constitutive relations allow us to link the fields and their respective inductions:

$$ \left\{ \begin{aligned} \mathbf{D}(\mathbf{x}, t) &= \varepsilon(\mathbf{x})\mathbf{E}(\mathbf{x}, t) \\ \mathbf{B}(\mathbf{x}, t) &= \mu(\mathbf{x})\mathbf{H}(\mathbf{x}, t) \end{aligned} \right. $$
where $\varepsilon(\mathbf{x})$ and $\mu(\mathbf{x})$ are respectively the permittivity and the permeability of the medium. With the constitutive relations, system (1) becomes:

$$ 
\left\{
\begin{aligned}
\varepsilon(\mathbf{x}) \frac{\partial \mathbf{E}}{\partial t} - \mathbf{rot}(\mathbf{H}) &= \mathbf{j} \\
\mu(\mathbf{x}) \frac{\partial \mathbf{H}}{\partial t} + \mathbf{rot}(\mathbf{E}) &= \mathbf{0} \\
\mathrm{div}(\varepsilon(\mathbf{x})\mathbf{E}) &= \rho \\
\mathrm{div}(\mu(\mathbf{x})\mathbf{H}) &= \mathbf{0}
\end{aligned}
\right.
\qquad (2) $$

Moreover, the last two equations are "auxiliary" in the sense that if they are verified for the initial condition then they will also be verified for any future state. One only needs to consider the first two equations to simulate the propagation of electromagnetic waves.

## 1.3 Boundary and interface conditions

In the context of metasurfaces, we often have discontinuous media since we use multiple materials to scatter the wave. For this next part, $S$ is an oriented surface and $\mathbf{n}$ is the unit normal to $S$ oriented outwards the domain $\Omega$.

### 1.3.1 Boundary conditions for a perfect electric conductor

For a perfect conductor with infinite conductivity, we get this system of equations:

$$ 
\left\{
\begin{aligned}
\mathbf{n} \times \mathbf{H}_S &= \mathbf{j}_S \\
\mathbf{n} \times \mathbf{E}_S &= \mathbf{0} \\
\mathbf{n} \cdot \mathbf{E}_S &= \rho_S / \varepsilon \\
\mathbf{n} \cdot \mathbf{H}_S &= 0
\end{aligned}
\right.
$$

where $\mathbf{E}_S$ and $\mathbf{H}_S$ are the electric and magnetic fields near the surface $S$. The surface current $\mathbf{j}_S$ and the surface charge $\rho_S$ are unknown.

### 1.3.2 Absorbing boundary condition

In our simulations, we will often consider the medium of propagation to be infinite. However, the domain we simulate must be finite. The absorbing condition of Silver-Müller acts as a virtual frontier that lets electromagnetics waves go through with few unwanted reflections:

$$ \mathbf{n} \times \mathbf{E} = -c\mu\mathbf{n} \times (\mathbf{n} \times \mathbf{H}) $$

where $c$ is the speed of light in the medium and $\mu$ is the medium's permeability.

### 1.3.3 Interface conditions

At the interface between two materials, 1 and 2, we have the system:

$$ 
\left\{
\begin{aligned}
\mathbf{n}_{1 \to 2} \cdot (\varepsilon_1 \mathbf{E}_1 - \varepsilon_2 \mathbf{E}_2) &= 0 \\
\mathbf{n}_{1 \to 2} \times (\mathbf{E}_1 - \mathbf{E}_2) &= \mathbf{0} \\
\mathbf{n}_{1 \to 2} \cdot (\mu_1 \mathbf{H}_1 - \mu_2 \mathbf{H}_2) &= 0 \\
\mathbf{n}_{1 \to 2} \times (\mathbf{H}_1 - \mathbf{H}_2) &= \mathbf{0}
\end{aligned}
\right.
$$
<figure><img src="image_3.png" /><figcaption>Figure 1: A schema for the interface between two media.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2: A schema showing the components and the wave vector in TM and TE modes.</figcaption></figure>

We illustrate such an interface with the normal in figure 1.

## 1.4 Transverse magnetic and transverse electric waves

By taking into account only the first two equations of (2), we have two equations for three components (x, y and z) giving us a total of six scalar equations.

If we suppose that the fields are constant following (Oz) (let us say we have an invariance following (Oz)), we get two decoupled systems of three equations each corresponding to two modes of propagation: transverse magnetic (TM) and transverse electric (TE).

The magnetic field of a $TM_z$ wave has a null component on z. Thus, its non-null components are $(E_z, H_x, H_y)$.

Equivalently, the electric field of a $TE_z$ wave has a null component on z. Thus, its non-null components are $(H_z, E_x, E_y)$.

A sketch of this configuration is shown in Figure 2.

### 1.4.1 System of equations for the propagation

$TM_z$ case:

$$ egin{cases} \varepsilon \frac{\partial E_z}{\partial t} - \frac{\partial H_y}{\partial x} + \frac{\partial H_x}{\partial y} = -j_z \\ \mu \frac{\partial H_x}{\partial t} + \frac{\partial E_z}{\partial y} = 0 \\ \mu \frac{\partial H_y}{\partial t} - \frac{\partial E_z}{\partial x} = 0 \end{cases} $$
TEz case:

$$
\left\{
\begin{aligned}
& \mu \frac{\partial H_z}{\partial t} + \frac{\partial E_y}{\partial x} - \frac{\partial E_x}{\partial y} = 0 \\
& \varepsilon \frac{\partial E_x}{\partial t} - \frac{\partial H_z}{\partial y} = -j_x \\
& \varepsilon \frac{\partial E_y}{\partial t} + \frac{\partial H_z}{\partial x} = -j_y
\end{aligned}
\right.
$$

1.4.2 Boundary and interface conditions

• Perfect conductor

- TMz case:

$E_z = 0$

- TEz case:

$n_y E_x - n_x E_y = 0$

with $\mathbf{n} = (n_x, n_y)^t$

• Absorbing condition of Silver-Müller

- TMz case:

$E_z = c\mu(n_y H_x - n_x H_y)$

- TEz case:

$c\mu H_z = -n_y E_x + n_x E_y$

• Interface between two materials

- TMz case:

$$
\begin{cases}
E_{z,1} - E_{z,2} = 0 \\
n_x(\mu_1 H_{x,1} - \mu_2 H_{x,2}) + n_y(\mu_1 H_{y,1} - \mu_2 H_{y,2}) = 0 \\
n_x(H_{y,1} - H_{y,2}) - n_y(H_{x,1} - H_{x,2}) = 0
\end{cases}
$$

- TEz case:

$$
\begin{cases}
H_{z,1} - H_{z,2} = 0 \\
n_x(\varepsilon_1 E_{x,1} - \varepsilon_2 E_{x,2}) + n_y(\varepsilon_1 E_{y,1} - \varepsilon_2 E_{y,2}) = 0 \\
n_x(E_{y,1} - E_{y,2}) - n_y(E_{x,1} - E_{x,2}) = 0
\end{cases}
$$

1.5 Renormalization

1.5.1 Usual adimensionalization for numerical methods

For a physical variable X, we define $\hat{X}$, the normalized variable. We describe the method used
in Viquerat's thesis[19] for adimensionning:

$$
\begin{cases}
\hat{\mathbf{H}} = Z_0 \mathbf{H} \\
\hat{\mathbf{E}} = \mathbf{E} \\
\hat{t} = c_0 t \\
\hat{\mathbf{j}} = Z_0 \mathbf{j}
\end{cases}
$$
where $Z_0 = \sqrt{\mu_0/\varepsilon_0} = 376.730 \, \Omega$ is the vacuum impedance with $\varepsilon_0 = 8.8541878188 \times 10^{-12} \, \text{F m}^{-1}$, the vacuum permittivity and $\mu_0 = 1.25663706127 \times 10^{-6} \, \text{H m}^{-1}$, the vacuum permeability. Finally, $c_0 = 1/\sqrt{\varepsilon_0\mu_0} = 2.99792458 \times 10^8 \, \text{m s}^{-1}$ is the speed of light in vacuum. The units of the physical variables are given in Table 1.

<table><caption>Table 1: Original and normalized units of the variables for Maxwell equations</caption><thead><tr><th></th><th>Normalized unit</th><th>Original units</th></tr></thead><tbody><tr><th><b>H</b></th><td>V m<sup>-1</sup></td><td>A m<sup>-1</sup></td></tr><tr><th><b>E</b></th><td>V m<sup>-1</sup></td><td>V m<sup>-1</sup></td></tr><tr><th><b>t</b></th><td>m</td><td>s</td></tr><tr><th><b>j</b></th><td>V m<sup>-2</sup></td><td>A m<sup>-2</sup></td></tr></tbody></table>

The system (2) becomes:

$$ \left\{ \begin{aligned} \varepsilon_r \frac{\partial \hat{\mathbf{E}}}{\partial \hat{t}} - \mathbf{rot}(\hat{\mathbf{H}}) &= \hat{\mathbf{j}} \\ \mu_r \frac{\partial \hat{\mathbf{H}}}{\partial \hat{t}} + \mathbf{rot}(\hat{\mathbf{E}}) &= \mathbf{0} \end{aligned} \right. \qquad (3) $$

where $\varepsilon_r$ and $\mu_r$ are respectively the relative permittivity and the relative permeability of the medium defined as $\varepsilon = \varepsilon_r\varepsilon_0$ and $\mu = \mu_r\mu_0$.

### 1.5.2 Further scaling for the inputs to a neural network

Let us recall that we intend to use a neural network to solve our problem, so ideally, we want the inputs, outputs and the loss functions to be properly scaled. The previous adimensionalization allows to scale the outputs and the partial differential equations which provide the loss. However, the inputs are far from being well-scaled.

When initializing a neural network, the goal is to avoid using weights that would cause vanishing or exploding gradients to appear[7]. To ensure this, the initialization scheme generates weights that keep the scale of signals constant throughout the neural network. Without further scaling, we would have an output between -1 and 1 (the electromagnetic field) while the input is between $-10^{-6}$ and $10^{-6}$ since the domain we would usually consider in a numerical simulation may be around a few micrometers large. These input and output scales are not compatible.

Hence, we need a scaling that does not change the partial differential equations but still scales the inputs of the neural network, we propose to use the following system:

$$ \left\{ \begin{aligned} \tilde{\mathbf{H}} &= \frac{Z_0}{E_0} \mathbf{H}, \quad \tilde{\mathbf{E}} = \frac{1}{E_0} \mathbf{E} \\ \tilde{t} &= \frac{c_0}{L_0} t \\ \tilde{\mathbf{j}} &= \frac{Z_0 L_0}{E_0} \mathbf{j} \\ \tilde{x} &= \frac{1}{L_0} x, \quad \tilde{y} = \frac{1}{L_0} y \end{aligned} \right. $$

Inria
where $E_0$ is the characteristic amplitude of the electrical field component chosen to make $\tilde{E}$ and $\tilde{H}$ well scaled while $L_0$ is the characteristic length chosen to make $\tilde{x}$, $\tilde{y}$ and $\tilde{t}$ well scaled.

This scaling still follows the system of equations 3 and allows us to scale our inputs. It also leaves our variables with no dimensions which is generally preferable.

### 1.5.3 Systems for transverse waves

We get these systems for the transverse cases: $TM_z$ case:

$$ egin{cases} 
\varepsilon_r \frac{\partial \tilde{E}_z}{\partial \tilde{t}} - \frac{\partial \tilde{H}_y}{\partial \tilde{x}} + \frac{\partial \tilde{H}_x}{\partial \tilde{y}} = -\tilde{j}_z \\ 
\mu_r \frac{\partial \tilde{H}_x}{\partial \tilde{t}} + \frac{\partial \tilde{E}_z}{\partial \tilde{y}} = 0 \\ 
\mu_r \frac{\partial \tilde{H}_y}{\partial \tilde{t}} - \frac{\partial \tilde{E}_z}{\partial \tilde{x}} = 0 
\end{cases} 
\quad (4) $$

$TE_z$ case:

$$ egin{cases} 
\mu_r \frac{\partial \tilde{H}_z}{\partial \tilde{t}} + \frac{\partial \tilde{E}_y}{\partial \tilde{x}} - \frac{\partial \tilde{E}_x}{\partial \tilde{y}} = 0 \\ 
\varepsilon_r \frac{\partial \tilde{E}_x}{\partial \tilde{t}} - \frac{\partial \tilde{H}_z}{\partial \tilde{y}} = -\tilde{j}_x \\ 
\varepsilon_r \frac{\partial \tilde{E}_y}{\partial \tilde{t}} + \frac{\partial \tilde{H}_z}{\partial \tilde{x}} = -\tilde{j}_y 
\end{cases} 
\quad (5) $$

### 1.5.4 Boundary and interface conditions

* Perfect conductor

- $TM_z$ case:

$$ \tilde{E}_z = 0 $$

- $TE_z$ case:

$$ n_y \tilde{E}_x - n_x \tilde{E}_y = 0 $$

* Absorbing condition of Silver-Müller

- $TM_z$ case:

$$ \tilde{E}_z = Z_r(n_y \tilde{H}_x - n_x \tilde{H}_y) $$

Where $Z_r = \frac{Z}{Z_0}$ is the relative impedance of the medium.

- $TE_z$ case:

$$ Z_r \tilde{H}_z = -n_y \tilde{E}_x + n_x \tilde{E}_y $$

* Interface between two materials

- $TM_z$ case:

$$ egin{cases} 
\tilde{E}_{z,1} - \tilde{E}_{z,2} = 0 \\ 
n_x(\mu_{r,1} \tilde{H}_{x,1} - \mu_{r,2} \tilde{H}_{x,2}) + n_y(\mu_{r,1} \tilde{H}_{y,1} - \mu_{r,2} \tilde{H}_{y,2}) = 0 \\ 
n_x(\tilde{H}_{y,1} - \tilde{H}_{y,2}) - n_y(\tilde{H}_{x,1} - \tilde{H}_{x,2}) = 0 
\end{cases} $$
– TEz case:

$$
\begin{cases}
\tilde{H}_{z,1} - \tilde{H}_{z,2} = 0 \\
n_x(\varepsilon_{r,1}\tilde{E}_{x,1} - \varepsilon_{r,2}\tilde{E}_{x,2}) + n_y(\varepsilon_{r,1}\tilde{E}_{y,1} - \varepsilon_{r,2}\tilde{E}_{y,2}) = 0 \\
n_x(\tilde{E}_{y,1} - \tilde{E}_{y,2}) - n_y(\tilde{E}_{x,1} - \tilde{E}_{x,2}) = 0
\end{cases}
$$

1.6 Physics-Informed Neural Networks

The concept of PINNs was popularized by Raissi et al.[12]. With PINNs, we are trying to solve a
system of partial differential equations with the associated initial and boundary conditions. Let
us first define the spatial domain as Ω ⊂ R n (here, n = 2) and the temporal domain as [0, T]
where T ∈ R ∗ $_{+}$. We want to find the solution u ∈ L $^{2}$(Ω × [0, T ] → R $^{m}$) (here, m = 3) such that:

$$
\left\{
\begin{array}{l}
\mathcal{N}(\mathbf{u})(\mathbf{x}, t) = \mathbf{f}(\mathbf{x}, t) \quad \text{for } \mathbf{x} \in \Omega \text{ and } t \in [0, T] \\
\mathcal{B}(\mathbf{u})(\mathbf{x}, t) = \mathbf{b}(\mathbf{x}, t) \quad \text{for } \mathbf{x} \in \partial\Omega \text{ and } t \in [0, T] \\
\mathbf{u}(\mathbf{x}, t=0) = \mathbf{u}_0(\mathbf{x}) \quad \text{for } \mathbf{x} \in \Omega
\end{array}
\right.
$$

where $\mathcal{N}$ is a linear differential operator and $\mathcal{B}$ is a linear operator in the present case.
We can represent the solution $\mathbf{u}$ as the outputs of a neural network which takes $\mathbf{x}$ and $t$ as inputs.
We call this approximation $\mathbf{u}_{\boldsymbol{\theta}}$ where $\boldsymbol{\theta}$ are the parameters of the neural network. From there,
we can use automatic differentiation to compute the partial derivatives of $\mathbf{u}_{\boldsymbol{\theta}}$ and use them to
compute residuals.

Moreover, for $g \in L^2(\Omega \to \mathbb{R})$:

$$
g = 0 \quad \text{almost everywhere} \iff \int_{x \in \Omega} g(x)^2 dx = 0
$$

We can then use Monte Carlo integration to approximate this integral:

$$
\int_{x \in \Omega} g(x)^2 dx \approx \frac{1}{N} \sum_{i=1}^{N} g(x_i)^2 \text{ for } \{x_i\}_{i=1}^{N} \text{ chosen at random in } \Omega
$$

This allows us to establish loss functions that we have to minimize in order to solve the problem:

$$
\mathcal{L}(\boldsymbol{\theta}) = \mathcal{L}_r(\boldsymbol{\theta}) + \mathcal{L}_b(\boldsymbol{\theta}) + \mathcal{L}_0(\boldsymbol{\theta})
$$

with the loss functions defined as:

$$
\mathcal{L}_r(\boldsymbol{\theta}) = \frac{1}{N_r} \sum_{i=1}^{N_r} |\mathcal{N}(\mathbf{u}_{\boldsymbol{\theta}})(\mathbf{x}_r^i, t_r^i) - \mathbf{f}(\mathbf{x}_r^i, t_r^i)|^2
$$

$$
\mathcal{L}_b(\boldsymbol{\theta}) = \frac{1}{N_b} \sum_{i=1}^{N_b} |\mathcal{B}(\mathbf{u}_{\boldsymbol{\theta}})(\mathbf{x}_b^i, t_b^i) - \mathbf{b}(\mathbf{x}_b^i, t_b^i)|^2
$$

$$
\mathcal{L}_0(\boldsymbol{\theta}) = \frac{1}{N_0} \sum_{i=1}^{N_0} |\mathbf{u}_{\boldsymbol{\theta}}(\mathbf{x}_0^i, t=0) - \mathbf{u}_0(\mathbf{x}_0^i)|^2
$$

where {$\mathbf{x}_r^i, t_r^i$}$_{i=1}^{N_r} \in (\Omega \times [0, T])^{N_r}$, {$\mathbf{x}_b^i, t_b^i$}$_{i=1}^{N_b} \in (\partial\Omega \times [0, T])^{N_b}$ and {$\mathbf{x}_0^i$}$_{i=1}^{N_0} \in \Omega^{N_0}$ correspond respectively to databases for the partial differential equations, for the boundary conditions and for the initial condition.

From there, the PINN is trained the same way as any other neural networks. We use gradient
back-propagation to update the parameters and try to minimize the loss.
<figure><img src="image_5.png" /><figcaption>Figure 3: The neural network has to model the solution so it takes as inputs spatio-temporal points and outputs the components of the solution. We use collocation points inside the domain, at the boundary and at the initial to compute the residuals which we have to minimize.</figcaption></figure>

## 2 Test cases

As explained in section 1.6, the NN is trained without using any annotated data. However, for validation purposes, we allow ourselves to compare the PINN prediction with the true solution. We only treat cases where the analytical solution is known but if no analytical solution is available, the classical approach is to use the result of a simulation with a numerical method for comparison.

Our metric of success is the L² relative error:

$$L_2(\mathbf{y}_{\text{true}}, \mathbf{y}_{\text{pred}}) = \frac{\|\mathbf{y}_{\text{true}} - \mathbf{y}_{\text{pred}}\|_2}{\|\mathbf{y}_{\text{true}}\|_2}$$

where ||.||₂ is the L² norm and y is a one dimensional vector made of the values of a scalar for multiple positions x and multiple times t either chosen randomly or on a grid. For example, for the TM mode, we check the L² relative error on $H_x$, $H_y$ and $E_z$. Meanwhile for the TE mode, we look at the L² relative error on $E_x$, $E_y$ and $H_z$.

### 2.1 C1: Stationary wave in a metallic cavity

We first consider the case of a square cavity of length L inside of a perfect conductor. A sketch of this problem is shown in Figure 4.

#### 2.1.1 Resonant mode (1,1) in a square cavity

In this case, the initial conditions are set as one of the resonant modes which are known for this geometry. We first model the case of the mode (1,1) in the TMz case:

$$\begin{cases} H_x = 0 \\ H_y = 0 \\ E_z = E_0 \cdot \sin\left(\frac{\pi x}{L}\right) \sin\left(\frac{\pi y}{L}\right) \end{cases}$$
<figure><img src="image_6.png" /><figcaption>Figure 4: A sketch of the resonant cavity case.</figcaption></figure>

Then, the solution is:

$$ egin{cases} H_x = -\frac{\sqrt{2} E_0}{2 Z} \sin\left(\frac{\pi x}{L}\right) \cos\left(\frac{\pi y}{L}\right) \sin(\omega t) \\ H_y = \frac{\sqrt{2} E_0}{2 Z} \cos\left(\frac{\pi x}{L}\right) \sin\left(\frac{\pi y}{L}\right) \sin(\omega t) \\ E_z = E_0 \sin\left(\frac{\pi x}{L}\right) \sin\left(\frac{\pi y}{L}\right) \cos(\omega t) \end{cases} $$

where $\omega = c\pi\sqrt{2}/L$.

When normalized, the initial condition becomes:

$$ egin{cases} H_x = 0 \\ H_y = 0 \\ E_z = \sin(\pi \tilde{x}) \sin(\pi \tilde{y}) \end{cases} $$

And the solution becomes:

$$ egin{cases} \tilde{H}_x = -\frac{\sqrt{2}}{2} \frac{1}{Z_r} \sin(\pi \tilde{x}) \cos(\pi \tilde{y}) \sin(\pi \sqrt{2} c_r \tilde{t}) \\ \tilde{H}_y = \frac{\sqrt{2}}{2} \frac{1}{Z_r} \cos(\pi \tilde{x}) \sin(\pi \tilde{y}) \sin(\pi \sqrt{2} c_r \tilde{t}) \\ \tilde{E}_z = \sin(\pi \tilde{x}) \sin(\pi \tilde{y}) \cos(\pi \sqrt{2} c_r \tilde{t}) \end{cases} $$

Where $c_r = c/c_0$ is the relative speed of light in the medium.

### 2.1.2 Higher-order resonant mode in a square cavity

Let $m$ be the number of antinodes along $x$ and $n$ be the number of antinodes along $y$. Inside of a rectangular cavity, the resonant modes are known for a transverse magnetic wave. Let us
consider the mode $(m, n)$ as the initial condition:

$$ 
\begin{cases} 
H_x = 0 \\ 
H_y = 0 \\ 
E_z = E_0 \cdot \sin\left(\frac{m\pi x}{L}\right) \sin\left(\frac{n\pi y}{L}\right) 
\end{cases} 
$$

Then, the solution is:

$$ 
\begin{cases} 
H_x = -\frac{n}{\sqrt{\gamma}} \frac{E_0}{Z} \sin\left(\frac{m\pi x}{L}\right) \cos\left(\frac{n\pi y}{L}\right) \sin(\omega t) \\ 
H_y = \frac{m}{\sqrt{\gamma}} \frac{E_0}{Z} \cos\left(\frac{m\pi x}{L}\right) \sin\left(\frac{n\pi y}{L}\right) \sin(\omega t) \\ 
E_z = E_0 \cdot \sin\left(\frac{m\pi x}{L}\right) \sin\left(\frac{n\pi y}{L}\right) \cos(\omega t) 
\end{cases} 
$$

where $\omega = c\pi\sqrt{\gamma}/L$ and $\gamma = m^2 + n^2$ as explained in *Classical electrodynamics* [8] in subsection 8.4 "Modes in a Rectangular Wave Guide".

When normalized, the initial condition becomes:

$$ 
\begin{cases} 
H_x = 0 \\ 
H_y = 0 \\ 
E_z = \sin(m\pi\tilde{x}) \sin(n\pi\tilde{y}) 
\end{cases} 
$$

And the solution becomes:

$$ 
\begin{cases} 
\tilde{H}_x = -\frac{n}{\sqrt{\gamma}} \frac{1}{Z_r} \sin(m\pi\tilde{x}) \cos(n\pi\tilde{y}) \sin(\pi\sqrt{\gamma}c_r\tilde{t}) \\ 
\tilde{H}_y = \frac{m}{\sqrt{\gamma}} \frac{1}{Z_r} \cos(m\pi\tilde{x}) \sin(n\pi\tilde{y}) \sin(\pi\sqrt{\gamma}c_r\tilde{t}) \\ 
\tilde{E}_z = \sin(m\pi\tilde{x}) \sin(n\pi\tilde{y}) \cos(\pi\sqrt{\gamma}c_r\tilde{t}) 
\end{cases} 
$$

## 2.2 C2: Propagation of a plane wave in a homogeneous domain

We consider once again a square domain of length L but this time, the propagation of an incident wave is supposed to take place in an open domain. A sketch of this configuration is shown in Figure 5.

In practice, we need to impose a boundary condition to artificially truncate the propagation domain. We use the Silver-Müller boundary condition and we take into account the incident wave:

$$ \mathbf{n} \times (\mathbf{E} - \mathbf{E}_{\text{inc}}) = -c_0 \mu_0 \mathbf{n} \times (\mathbf{n} \times (\mathbf{H} - \mathbf{H}_{\text{inc}})) $$

where $\mathbf{E}_{\text{inc}}$ and $\mathbf{H}_{\text{inc}}$ are the electric and magnetic fields of the incident wave near the boundary of the domain. Note that when using the Silver-Müller boundary condition, the problem is well posed.

For a $TM_z$ wave, the normalized constraint is:

$$ \tilde{E}_z - n_y \tilde{H}_x + n_x \tilde{H}_y = \tilde{E}_{z,inc} - n_y \tilde{H}_{x,inc} + n_x \tilde{H}_{y,inc} $$
<figure><img src="image_7.png" /><figcaption>Figure 5: A sketch for the case of the propagation through an empty domain.</figcaption></figure>

While for a TEz wave, we have:

$$ \tilde{H}_z + n_y \tilde{E}_x - n_x \tilde{E}_y = \tilde{H}_{z,inc} + n_y \tilde{E}_{x,inc} - n_x \tilde{E}_{y,inc} $$

As for the incident wave, a plane wave is defined by a single parameter differentiable function f as well as its direction of propagation $k_u = (k_{ux}, k_{uy}, 0)^T$ which is a unit vector. Such plane waves can be written as:

$$ u(x, y, t) = f(k_{ux}x + k_{uy}y - ct) $$

In the TMz case, let $E_z = E_0.f(k_{ux}x + k_{uy}y - ct)$, then, the incident wave follows:

$$ \left\{ \begin{aligned} H_x &= k_{uy} \frac{E_0}{Z} f(k_{ux}x + k_{uy}y - ct) \\ H_y &= -k_{ux} \frac{E_0}{Z} f(k_{ux}x + k_{uy}y - ct) \\ E_z &= E_0.f(k_{ux}x + k_{uy}y - ct) \end{aligned} \right. $$

In the normalized system, this becomes:

$$ \left\{ \begin{aligned} \tilde{H}_x &= k_{uy} \frac{1}{Z_r} f(L \times [k_{ux}\tilde{x} + k_{uy}\tilde{y} - c_r \tilde{t}]) \\ \tilde{H}_y &= -k_{ux} \frac{1}{Z_r} f(L \times [k_{ux}\tilde{x} + k_{uy}\tilde{y} - c_r \tilde{t}]) \\ \tilde{E}_z &= f(L \times [k_{ux}\tilde{x} + k_{uy}\tilde{y} - c_r \tilde{t}]) \end{aligned} \right. $$

For the TEz case, let us say that $H_z = E_0/Z \times f(k_{ux}x + k_{uy}y - ct)$, then, the incident wave
follows:

$$
\left\{
\begin{aligned}
E_x &= -k_{uy} E_0 \cdot f(k_{ux} x + k_{uy} y - ct) \\
E_y &= k_{ux} E_0 \cdot f(k_{ux} x + k_{uy} y - ct) \\
H_z &= \frac{E_0}{Z} \cdot f(k_{ux} x + k_{uy} y - ct)
\end{aligned}
\right.
$$

In the normalized system, this becomes:

$$
\left\{
\begin{aligned}
\tilde{E}_x &= -k_{uy} \cdot f(L \times [k_{ux}\tilde{x} + k_{uy}\tilde{y} - c_r \tilde{t}]) \\
\tilde{E}_y &= k_{ux} \cdot f(L \times [k_{ux}\tilde{x} + k_{uy}\tilde{y} - c_r \tilde{t}]) \\
\tilde{H}_z &= \frac{1}{Z_r} f(L \times [k_{ux}\tilde{x} + k_{uy}\tilde{y} - c_r \tilde{t}])
\end{aligned}
\right.
$$

2.2.1 Incident monochromatic plane wave

For a monochromatic wave we have f(.) = cos(k.) where k is the wave number. Hence, the
following system:

$$
\left\{
\begin{aligned}
H_x &= k_{uy} \frac{E_0}{Z} \cos(k \times [k_{ux}x + k_{uy}y - ct]) \\
H_y &= -k_{ux} \frac{E_0}{Z} \cos(k \times [k_{ux}x + k_{uy}y - ct]) \\
E_z &= E_0 \cdot \cos(k \times [k_{ux}x + k_{uy}y - ct])
\end{aligned}
\right.
$$

We normalize:

$$
\left\{
\begin{aligned}
\tilde{H}_x &= k_{uy} \frac{1}{Z_r} \cos(\tilde{k} \times [k_{ux}\tilde{x} + k_{uy}\tilde{y} - c_r \tilde{t}]) \\
\tilde{H}_y &= -k_{ux} \frac{1}{Z_r} \cos(\tilde{k} \times [k_{ux}\tilde{x} + k_{uy}\tilde{y} - c_r \tilde{t}]) \\
\tilde{E}_z &= \cos(\tilde{k} \times [k_{ux}\tilde{x} + k_{uy}\tilde{y} - c_r \tilde{t}])
\end{aligned}
\right.
$$

where $\tilde{k} = kL$ is the normalized wave vector.

2.2.2 Incident gaussian plane wave

Doing a simulation in the time-domain with monochromatic waves is a bit of a waste since it is
more appropriate to a frequency-domain setting. The point of being in the time-domain is being
able to simulate signals of continuous spectrum. A very common choice is the gaussian plane
wave which is the product of a cosine with a gaussian function. Thus we have:

$$
f(\cdot) = \exp\left(\frac{-(\cdot - d_0)^2}{2\sigma^2}\right) \cos(k\cdot)
$$

where $d_0$ is the displacement of the center of the envelop at $t = 0$ and $\sigma$ is the standard deviation of the gaussian envelop.
Hence the following system:

$$
\left\{
\begin{aligned}
H_x &= \frac{k_{u,y}}{Z} E \cdot \exp \left( \frac{-(\mathbf{k}_u \cdot \mathbf{x} - x_0 - ct)^2}{2\sigma^2} \right) \cos(k(\mathbf{k}_u \cdot \mathbf{x} - ct)) \\
H_y &= -\frac{k_{u,x}}{Z} E \cdot \exp \left( \frac{-(\mathbf{k}_u \cdot \mathbf{x} - x_0 - ct)^2}{2\sigma^2} \right) \cos(k(\mathbf{k}_u \cdot \mathbf{x} - ct)) \\
E_z &= E \cdot \exp \left( \frac{-(\mathbf{k}_u \cdot \mathbf{x} - x_0 - ct)^2}{2\sigma^2} \right) \cos(k(\mathbf{k}_u \cdot \mathbf{x} - ct))
\end{aligned}
\right.
$$

We normalize this system and we get:

$$
\left\{
\begin{aligned}
\tilde{H}_x &= \frac{k_{u,y}}{Z_r} \exp \left( \frac{-(\mathbf{k}_u \cdot \tilde{\mathbf{x}} - \tilde{x}_0 - c_r \tilde{t})^2}{2\tilde{\sigma}^2} \right) \cos(\tilde{k}(\mathbf{k}_u \cdot \tilde{\mathbf{x}} - c_r \tilde{t})) \\
\tilde{H}_y &= -\frac{k_{u,x}}{Z_r} \exp \left( \frac{-(\mathbf{k}_u \cdot \tilde{\mathbf{x}} - \tilde{x}_0 - c_r \tilde{t})^2}{2\tilde{\sigma}^2} \right) \cos(\tilde{k}(\mathbf{k}_u \cdot \tilde{\mathbf{x}} - c_r \tilde{t})) \\
\tilde{E}_z &= \exp \left( \frac{-(\mathbf{k}_u \cdot \tilde{\mathbf{x}} - \tilde{x}_0 - c_r \tilde{t})^2}{2\tilde{\sigma}^2} \right) \cos(\tilde{k}(\mathbf{k}_u \cdot \tilde{\mathbf{x}} - c_r \tilde{t}))
\end{aligned}
\right.
$$

where $\tilde{\sigma} = \sigma/L$ is the normalized standard deviation and $\tilde{x}_0 = x_0/L$ is the normalized starting position.

## 2.3 C3: Flat interface between two dielectric materials

In most applications, materials with different permittivity are used to control the propagation of an electromagnetic wave. Here, we want to study if the PINN can accurately model the interface between two media. There are absorbing boundary conditions on the left and right of the domain and periodic boundary conditions on the lower and upper limits. The incident wave comes from the left with normal incidence to the interface. A sketch of this configuration is shown in Figure 6.

The wave has normal incidence and the interface is parallel to the axis (Oy). This means that the propagates according to x only. The analytical solution is the following (proof in appendix A.1):

$$
\left\{
\begin{aligned}
f_t(c_2 t - x) &= \frac{2n_1}{n_1 + n_2} f_i \left( c_1 t - \frac{c_1}{c_2} x + \frac{c_1 - c_2}{c_2} x_{int} \right) \\
f_r(c_1 t + x) &= \frac{n_1 - n_2}{n_1 + n_2} f_i(c_1 t + x - 2x_{int})
\end{aligned}
\right.
$$

In general, the refractive indices are defined as $n_i = c_0/c_i = \sqrt{\varepsilon_i \mu_i}/\sqrt{\varepsilon_0 \mu_0}$.

## 2.4 C4: Dielectric slab

The boundary conditions are the same as before and the incident wave still has normal incidence. A sketch of the configuration is shown in Figure 7.

Once again, the wave propagates only alongside x since it has normal incidence. The analytical solution can be established for a monochromatic wave. We will work in the Fourier space with the complex notations. In total, there are 5 waves of interest:
<figure><img src="image_8.png" /><figcaption>Figure 6: A sketch of the flat interface case. The incident wave is plane and comes with normal incidence.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: A sketch of the dielectric slab case. The incident wave is plane monochromatic and comes with normal incidence.</figcaption></figure>
• The incident wave in medium 1: $\bar{E}_{i,1} = \bar{A}_{i,1} \exp(i\omega \left[t - \frac{x}{c_1}\right])$

• The reflected wave in medium 1: $\bar{E}_{r,1} = \bar{A}_{r,1} \exp(i\omega \left[t + \frac{x}{c_1}\right])$

• The transmitted wave in medium 2: $\bar{E}_{t,2} = \bar{A}_{t,2} \exp(i\omega \left[t - \frac{x}{c_2}\right])$

• The reflected wave in medium 2: $\bar{E}_{r,2} = \bar{A}_{r,2} \exp(i\omega \left[t + \frac{x}{c_2}\right])$

• The transmitted wave in medium 3: $\bar{E}_{t,3} = \bar{A}_{t,3} \exp(i\omega \left[t - \frac{x}{c_3}\right])$

The $\bar{A}$'s are constant complex numbers. Since we work in the Fourier space, each $\bar{A}$ represents the amplitude and the phase shift of each wave.

The analytical solution follows this system (proof in appendix A.2):

$$ 
\left\{ 
\begin{aligned}
\bar{A}_{t,2} &= \frac{T_{1 \to 2} \exp\left(i\omega\left[\frac{x_1}{c_2} - \frac{x_1}{c_1}\right]\right)}{1 + R_{1 \to 2} R_{2 \to 3} \exp\left(i\omega\left[\frac{2x_1}{c_2} - \frac{2x_2}{c_2}\right]\right)} \bar{A}_{t,1} \\
\bar{A}_{r,1} &= \frac{R_{1 \to 2} + R_{2 \to 3} \exp\left(i\omega\left[\frac{2x_1}{c_2} - \frac{2x_2}{c_2}\right]\right)}{T_{1 \to 2} \exp\left(i\omega\left[\frac{x_1}{c_2} + \frac{x_1}{c_1}\right]\right)} \bar{A}_{t,2} \\
\bar{A}_{t,3} &= T_{2 \to 3} \exp\left(i\omega\left[\frac{x_2}{c_3} - \frac{x_2}{c_2}\right]\right) \bar{A}_{t,2} \\
\bar{A}_{r,2} &= R_{2 \to 3} \exp\left(-i\omega\frac{2x_2}{c_2}\right) \bar{A}_{t,2}
\end{aligned}
\right. 
$$

With this system, we have the amplitude and phase shift of each wave. But once again, we only have the solution in the frequency-domain but we work in the time-domain. If the incident wave is monochromatic, we can simply convert the analytical solution to the temporal domain by taking the real part of the complex solution.

For incident waves with continuous spectrum like the gaussian plane wave, it is impossible to do something similar. So, instead, one can compare the solution and the prediction for multiple frequencies by using the Fourier transform on the prediction and comparing those spectra with our solution in the Fourier space.

## 2.5 C5: Diffraction of an incident wave by a dielectric cylinder

In the previous cases, we have considered situations for which the wave was invariant alongside the axis (Oy) which made the problem much simpler. On our way toward more complicated applications, we choose to model a cylinder since we have an analytical solution to this problem in the frequency-domain. A sketch of the configuration is shown in Figure 8.

Once again, we cannot compare the exact and the predicted solutions since they are not both in the same physical setting. To fix this, we have to convert either the solution to the temporal domain or the prediction to the frequency domain. We apply the same methods as in subsection
<figure><img src="image_10.png" /><figcaption>Figure 8: A sketch of the diffraction by a dielectric cylinder case.</figcaption></figure>

2.4.

We will consider two cases, one where the cylinder is dielectric with a different permittivity and one where the cylinder is a perfect electric conductor. The analytical solutions are taken from the appendix A of [5].
<table><caption>Table 2: Examples of sampling for each method as well as their convergence rates.</caption><thead><tr><th>Method</th><th>Grid</th><th>Pseudo-random</th><th>Halton</th><th>Hammersley</th></tr></thead><tbody><tr><td>100 points</td><td></td><td></td><td></td><td></td></tr><tr><td>400 points</td><td></td><td></td><td></td><td></td></tr><tr><td>900 points</td><td></td><td></td><td></td><td></td></tr><tr><td>Error convergence</td><td>O(1/N<sup>1/d</sup>)</td><td>O(1/√N)</td><td>O(ln(N)<sup>d</sup>/N)</td><td>O(ln(N)<sup>d</sup>/N)</td></tr></tbody></table>

# 3 Sampling of the collocation points

## 3.1 Sampling methods

Training the PINN requires to estimate the loss functions through Monte Carlo integration. Hence, we need a method to sample the points where we evaluate the residuals. Many sampling methods are available, but we can classify them in three families:

* **Grid.** We simply take a regular grid over the domain. This method is very inefficient for dimensions higher than 2 since the error of the Monte Carlo estimate decreases as $O(N^{-1/d})$ where $d$ is the number of dimensions. This method also has the problem of being deterministic so it is impossible to sample multiple datasets.

* **Pseudo-random.** Points are sampled according to a uniform distribution over the domain. The error converges as $O(N^{-1/2})$ regardless of the number of dimensions. It is easy to use and can generate large datasets on complex geometries quickly. Its main weakness is the lack of consistent discrepancy which results in some parts of the domain getting overcrowded with samples while others are left barren.

* **Low-discrepancy sequences.** Points are sampled from a deterministic low-discrepancy sequence to which we add some noise to get randomized samples. The error decreases as $O(\ln(N)^d N^{-1})$ which is almost linear for low numbers of dimensions. As one can see, the density of samples seem much more uniform over the domain.

The above convergence rates are taken from Morokoff et al. [10]. We illustrate the methods in Table 2. There are multiple low-discrepancy sequences available and we chose to use the Halton and Hammersley sequences for our purposes.
## 3.2 Sampling strategies

Now, there is the matter of the sampling strategy. For example, do we sample the collocation points a single time or do we update the dataset during training? If so, how often? Etc. Here are some of those strategies as presented by Wu et al. [22]:

* Fixed residual points: we sample once at the start of the training.

* Uniform resampling: every N epochs, we sample a new dataset to replace the old one.

* Non-uniform resampling: once again, we sample a new dataset every N epochs but instead of using a uniform distribution, we use the residuals to sample our new dataset:

- RAR-G: Residual adaptive refinement with greed, we sample M points from the domain and we select the m points with the largest residuals and add them to the database.

- RAR-D: Residual adaptive refinement with distribution, we sample m points according to the probability density function defined afterwards in equation 3.2. Then, we add them to the dataset.

- RAD: Residual adaptive distribution, we sample m points according to the probability density function defined in equation 3.2. Then we replace the previous dataset with these new points.

Let us define the probability density function used for the residual adaptive distribution on a domain $\Omega^*$:

$$p(x) \propto \frac{|r(x)|^k}{\mathbb{E}_{\Omega^*}[|r(x)|^k]} + c \quad \text{for } \mathbf{x} \in \Omega^*$$

where $r(x)$ is the residual at $x$. Additionally, we choose $k > 0$ and $c > 0$.

In general, using more complex sampling strategies can allow to get better results with smaller databases which means quicker training. However, these strategies do require to tinker some hyperparameters to have optimal results.

In practice, we will prefer using strategies that do not change the size of the dataset. It allows to control the training time more effectively, prevents the dataset from growing too much and avoids the "inertia" caused by adding new samples without removing any. Finally, it also makes one less hyperparameter to worry about. Hence, we will focus on three strategies: fixed sampling, uniform resampling and residual adaptive distribution.
# 4 Neural networks architecture

Let us discuss possible choices for the architecture of the neural network ranging from vanilla to specifically engineered architectures for our problems.

## 4.1 Vanilla MLP

The Multilayer Perceptron (MLP) is the most commonly used architecture in Deep Learning. It consists of multiple layers of so called neurons that take as inputs a linear combination of the outputs of the previous layer and apply an activation function:

$$a_i^{(l)} = \sigma \left( \sum_j \omega_{i,j}^{(l)} a_j^{(l-1)} + b_{i,j}^{(l)} \right)$$

where $a_i^{(l)}$ is the i-th output of the l-th layer, $\omega_{i,j}^{(l)}$ are the weights for the linear combination and $b_{i,j}^{(l)}$ are the biases. Finally, $\sigma$ is the activation function. In our case, we use hyperbolic tangent. These kinds of neural networks can approximate very complicated functions if given enough layers but they struggle when used as PINNs. This is due to many reasons, one of these being spectral bias as explained by Wang et al. [20]. Spectral bias is the tendency of neural networks to fit lower frequencies much faster than higher frequencies which causes problems when the expected solution exhibits said high frequencies.

## 4.2 SIREN

The so-called Sinusoidal Representation Networks (SIREN) were introduced by Sitzmann et al. [15] as a way to create implicit neural representations of shapes and images. The main difference between this architecture and a regular MLP is that of the activation function, which is the sine function. It also multiplies the input and output of each layer by coefficients to encourage learning more diverse frequencies. This allows the network to mitigate spectral bias.

## 4.3 Fourier embedding

As stated before, spectral bias causes dense neural networks to learn lower frequencies first and struggle to learn higher frequencies. To help with this issue, Wang et al. [21] suggests using Random Fourier features for PINNs as presented by Tancik et al. [17] in the context of 3D imagery. The inputs are embedded using the following encoding:

$$\gamma(\mathbf{X}) = \begin{bmatrix} \cos(B\mathbf{X}) \\ \sin(B\mathbf{X}) \end{bmatrix}$$

where $\mathbf{X}$ are the inputs of the neural network, $B$ is a matrix whose coefficients are sampled from a Gaussian distribution $N(0, \sigma^2)$ with $\sigma$ being a hyper-parameter. This encourages the neural network to learn various frequencies, which helps to negate spectral bias.

## 4.4 Imposing periodicity

Dong et al. [4] show how to impose periodic boundary conditions by modifying the inputs of the neural network. Let us say that we want to impose periodicity on time to the neural network.
<figure><img src="image_11.png" /><figcaption>Figure 9: Domain decomposition for heterogeneous media: we divide the medium in homogeneous subdomains and use a different neural network for each.</figcaption></figure>

It needs to satisfy:

$$ \forall \mathbf{x} \in \Omega, \forall t, \mathbf{u}(\mathbf{x}, t + P) = \mathbf{u}(\mathbf{x}, t) $$

where P is the period we want to impose.

All we need to do is embed the time feature using:

$$ \mathbf{p}(t) = (\cos(\omega t), \sin(\omega t)) $$

where $\omega = \frac{2\pi}{P}$. The input of the neural network becomes $(\mathbf{x}, \mathbf{p}(t)) = (x, y, \cos(\omega t), \sin(\omega t))$.

It is also possible to enforce periodicity along multiple dimensions at once.

## 4.5 Domain decomposition for heterogeneous media

In realistic applications, we are often confronted to heterogeneous media. This creates additional problems for the PINN as the solution follows different systems of partial differential equations in the domain. This results in solutions that are not $C^\infty$ while the output of a neural network with $C^\infty$ activation functions is always $C^\infty$.

To fix this, we propose using domain decomposition so that each subdomain is homogeneous. We then add the interface conditions between each material to the loss function. That way, the solution computed by the neural networks is $C^\infty$ inside each material while it can be irregular at the interfaces. We illustrate this approach in figure 9.

## 4.6 Scattering problem

A common class of problem for us will be the scattering of an incoming wave by a heterogeneous structure. An advantageous reformulation of the problem is to have the PINN solve the scattering problem instead. This involves decomposing the total field as two components, the incident field and the scattered field: $\vec{F}_{total} = \vec{F}_{inc} + \vec{F}_{scatter}$. A classic resolution consists in solving for $\vec{F}_{total}$ while in the scattered problem we are trying to model $\vec{F}_{scatter}$. In classical numerical methods,
there are Total Field/Scattered Field (TFSF) surfaces where this decomposition is applied and where observables are computed. For PINNs, the goal is instead to solve a simpler task by modeling only the scattered field.
# 5 Training dynamics

In this part, we study the way PINNs train through gradient descent to explain why it is hard to train these models. We then try to alleviate those problems with new algorithms accounting for those difficulties.

## 5.1 Convergence and conditioning

Let us follow the reasoning presented by Ryck et al. [14]. To train neural networks, we use gradient descent which gives us:

$$ \theta_{k+1} = \theta_k - \eta \nabla_{\theta} L(\theta_k) \quad (6) $$

where $\eta$ is the learning rate and $\theta_k$ are the parameters of the model at step $k$ of the training. First, the second-order Taylor expansion of $u_{\theta_k}(x)$ gives:

$$ u_{\theta_k}(x) = u_{\theta_0}(x) + \nabla_{\theta} u_{\theta_0}(x)^T (\theta_k - \theta_0) + \frac{1}{2} (\theta_k - \theta_0)^T H_k (\theta_k - \theta_0) $$

where $\theta_0$ is a fixed set of parameters (like the initial parameters) and $H_k = Hess_{\theta}(u_{\alpha\theta_k+(1-\alpha)\theta_0}(x))$ for some $\alpha$ between 0 and 1.

We take the hypothesis that $H_k$ is close to null (this is very similar to being in the NTK regime). We now replace $u_{\theta_k}(x)$ by its Taylor expansion in the losses:

$$ \left\{ \begin{aligned} R(\theta_k) &= \frac{1}{2} \int_{x \in \Omega} (\mathcal{L}(u_{\theta_0})(x) + \mathcal{L}(\nabla_{\theta} u_{\theta_0}))(x)^T (\theta_k - \theta_0) - f(x))^2 dx \\ B(\theta_k) &= \frac{1}{2} \int_{x \in \partial \Omega} (u_{\theta_0}(x) + \nabla_{\theta} u_{\theta_0})(x)^T (\theta_k - \theta_0) - g(x))^2 dx \end{aligned} \right. $$

Now, we compute their gradients. First for $\nabla_{\theta} R(\theta_k)$:

$$ \begin{aligned} \nabla_{\theta} R(\theta_k) &= \int_{x \in \Omega} \nabla_{\theta} (\mathcal{L}(u_{\theta_0})(x) + \mathcal{L}(\nabla_{\theta} u_{\theta_0}))(x)^T (\theta_k - \theta_0) - f(x)) \\ & \quad \cdot (\mathcal{L}(u_{\theta_0})(x) + \mathcal{L}(\nabla_{\theta} u_{\theta_0}))(x)^T (\theta_k - \theta_0) - f(x)) dx \\ &= \int_{x \in \Omega} \mathcal{L}(\nabla_{\theta} u_{\theta_0})(x) \cdot (\mathcal{L}(u_{\theta_0})(x) + \mathcal{L}(\nabla_{\theta} u_{\theta_0}))(x)^T (\theta_k - \theta_0) - f(x)) dx \\ &= \int_{x \in \Omega} [\mathcal{L}(\nabla_{\theta} u_{\theta_0})(x) \mathcal{L}(\nabla_{\theta} u_{\theta_0})(x)^T (\theta_k - \theta_0) + \mathcal{L}(\nabla_{\theta} u_{\theta_0})(x) (\mathcal{L}(u_{\theta_0})(x) - f(x))] dx \\ &= \left[ \int_{x \in \Omega} \mathcal{L}(\nabla_{\theta} u_{\theta_0})(x) \mathcal{L}(\nabla_{\theta} u_{\theta_0})(x)^T dx \right] (\theta_k - \theta_0) + \left[ \int_{x \in \Omega} \mathcal{L}(\nabla_{\theta} u_{\theta_0})(x) (\mathcal{L}(u_{\theta_0})(x) - f(x)) \right] \end{aligned} $$

In the same way:

$$ \begin{aligned} \nabla_{\theta} B(\theta_k) &= \int_{x \in \partial \Omega} \nabla_{\theta} (u_{\theta_0}(x) + \nabla_{\theta} u_{\theta_0}(x)^T (\theta_k - \theta_0) - g(x)) \\ & \quad \cdot (u_{\theta_0}(x) + \nabla_{\theta} u_{\theta_0}(x)^T (\theta_k - \theta_0) - g(x)) dx \\ &= \left[ \int_{x \in \partial \Omega} \nabla_{\theta} u_{\theta_0}(x) \nabla_{\theta} u_{\theta_0}(x)^T dx \right] (\theta_k - \theta_0) + \left[ \int_{x \in \partial \Omega} \nabla_{\theta} u_{\theta_0}(x) (u_{\theta_0}(x) - g(x)) dx \right] \end{aligned} $$
Hence, we have:

$$ 
\begin{aligned}
\nabla_{\boldsymbol{\theta}} L(\boldsymbol{\theta}_k) &= \left[ \int_{x \in \Omega} \mathcal{L}(\nabla_{\boldsymbol{\theta}} u_{\boldsymbol{\theta}_0})(x) \mathcal{L}(\nabla_{\boldsymbol{\theta}} u_{\boldsymbol{\theta}_0})(x)^T dx + \int_{x \in \partial \Omega} \nabla_{\boldsymbol{\theta}} u_{\boldsymbol{\theta}_0}(x) \nabla_{\boldsymbol{\theta}} u_{\boldsymbol{\theta}_0}(x)^T dx \right] (\boldsymbol{\theta}_k - \boldsymbol{\theta}_0) \\
&\quad + \left[ \int_{x \in \Omega} \mathcal{L}(\nabla_{\boldsymbol{\theta}} u_{\boldsymbol{\theta}_0})(x) (\mathcal{L}(u_{\boldsymbol{\theta}_0})(x) - f(x)) + \int_{x \in \partial \Omega} \nabla_{\boldsymbol{\theta}} u_{\boldsymbol{\theta}_0}(x) (u_{\boldsymbol{\theta}_0}(x) - g(x)) dx \right] \\
&= \mathbf{A}_{\boldsymbol{\theta}_0}(\boldsymbol{\theta}_k - \boldsymbol{\theta}_0) + \mathbf{B}_{\boldsymbol{\theta}_0}
\end{aligned}
$$

where $\mathbf{A}_{\boldsymbol{\theta}_0}$ and $\mathbf{B}_{\boldsymbol{\theta}_0}$ only depend on $\boldsymbol{\theta}_0$ and $\mathbf{B}_{\boldsymbol{\theta}_0} = \nabla_{\boldsymbol{\theta}} L(\boldsymbol{\theta}_0)$. The equation (6) becomes:

$$ \boldsymbol{\theta}_{k+1} = \boldsymbol{\theta}_k - \eta \mathbf{A}_{\boldsymbol{\theta}_0}(\boldsymbol{\theta}_k - \boldsymbol{\theta}_0) - \eta \mathbf{B}_{\boldsymbol{\theta}_0} $$

Which we can rewrite:

$$ \boldsymbol{\theta}_{k+1} = (\mathbf{I} - \eta \mathbf{A}_{\boldsymbol{\theta}_0}) \boldsymbol{\theta}_k + \eta \mathbf{A}_{\boldsymbol{\theta}_0} \boldsymbol{\theta}_0 - \eta \mathbf{B}_{\boldsymbol{\theta}_0} $$

Let us suppose $\mathbf{A}_{\boldsymbol{\theta}_0}$ is invertible, we have:

$$ \mathbf{v}_{k+1} = (\mathbf{I} - \eta \mathbf{A}_{\boldsymbol{\theta}_0}) \mathbf{v}_k $$

where $\mathbf{v}_k = \boldsymbol{\theta}_k - \boldsymbol{\theta}_0 + \mathbf{A}_{\boldsymbol{\theta}_0}^{-1} \mathbf{B}_{\boldsymbol{\theta}_0}$.

It then follows that:

$$ \mathbf{v}_k = (\mathbf{I} - \eta \mathbf{A}_{\boldsymbol{\theta}_0})^k \mathbf{v}_0 $$

As shown in the appendix A.4 of [14], we have:

$$ 
\begin{aligned}
\|\mathbf{v}_k\|_2 &= \|(\mathbf{I} - \eta \mathbf{A}_{\boldsymbol{\theta}_0})^k \mathbf{v}_0\|_2 \\
&\le \|(\mathbf{I} - \eta \mathbf{A}_{\boldsymbol{\theta}_0})\|_2^k \cdot \|\mathbf{v}_0\|_2 \\
&\le \max_j |1 - \eta \lambda_j(\mathbf{A}_{\boldsymbol{\theta}_0})|^k \cdot \|\mathbf{v}_0\|_2
\end{aligned}
$$

where the set $\{\lambda_j(\mathbf{A}_{\boldsymbol{\theta}_0})\}$ refers to the eigenvalues of $\mathbf{A}_{\boldsymbol{\theta}_0}$.

Let $0 < c < 1$. We set $\eta = c/\lambda_{\max}(\mathbf{A}_{\boldsymbol{\theta}_0})$, we have:

$$ 
\begin{aligned}
\|\mathbf{v}_k\|_2 &\le \left(1 - \frac{c \lambda_{\min}(\mathbf{A}_{\boldsymbol{\theta}_0})}{\lambda_{\max}(\mathbf{A}_{\boldsymbol{\theta}_0})}\right)^k \|\mathbf{v}_0\|_2 \\
&= \left(1 - \frac{c}{\kappa(\mathbf{A}_{\boldsymbol{\theta}_0})}\right)^k \|\mathbf{v}_0\|_2
\end{aligned}
$$

where $\kappa(\mathbf{A}_{\boldsymbol{\theta}_0}) = \lambda_{\max}(\mathbf{A}_{\boldsymbol{\theta}_0})/\lambda_{\min}(\mathbf{A}_{\boldsymbol{\theta}_0})$ is the condition number of $\mathbf{A}_{\boldsymbol{\theta}_0}$.

In conclusion, $\|\mathbf{v}_k\|_2$ tends to 0 and the condition number limits the speed at which we converge.

We also conclude that $\boldsymbol{\theta}_k$ converges to $\boldsymbol{\theta}_0 - \mathbf{A}_{\boldsymbol{\theta}_0}^{-1} \mathbf{B}_{\boldsymbol{\theta}_0}$.

## 5.2 Energy Natural Gradients

Müller et al. [11] introduces Energy Natural Gradients (ENG) in the context of PINNs which aims to find more suitable gradients according to the energy function. The new gradient descent algorithm is the following:

$$ \boldsymbol{\theta}_{k+1} = \boldsymbol{\theta}_k - \alpha \mathbf{A}_{\boldsymbol{\theta}_k}^{-1} \nabla_{\boldsymbol{\theta}} L(\boldsymbol{\theta}_k) $$

where $\alpha$ is the learning rate.

This solves our concerns about convergence as we directly take the limit of the sequence studied
previously. The condition number still plays a big role as we have to invert the matrix. More precisely, we view $A_{\theta_k}^{-1} \nabla_{\theta} L(\theta_k)$ as the solution $\mathbf{x}$ to the system $A_{\theta_k} \mathbf{x} = \nabla_{\theta} L(\theta_k)$, we use the least square solution for this system.

To be as efficient as possible and to avoid computing the least square solution too many times, we do a line search on the learning rate (either we sweep across possible learning rates and select the best one or we use the Wolfe conditions as stopping criteria [6]).

Another way we can tweak our algorithm is by introducing a damping parameter to the matrix we invert. Instead of using $A_{\theta_k}$, we consider $A_{\theta_k}^{\lambda} = A_{\theta_k} + \lambda \mathbf{I}$.

Algorithm 1 sums up each step, from the computation of the matrices to the application of the gradients.

<caption>**Algorithm 1** Revised gradient descent with line search</caption>

$\theta \leftarrow \theta_0$
**for** epoch = 1, ..., N **do**
    $\mathbf{A} \leftarrow \mathbf{A}_{\theta} + \lambda \mathbf{I}$
    $\mathbf{B} \leftarrow \nabla_{\theta} L(\theta)$
    $\mathbf{g} \leftarrow \text{LeastSquare}(\mathbf{A}, \mathbf{B})$
    $\alpha \leftarrow \text{LineSearch}(L, \theta, \mathbf{g})$
    $\theta \leftarrow \theta - \alpha \mathbf{g}$
**end for**
# 6 Results

Let us now apply the techniques presented in the earlier sections on our selected test cases. First, we will perform ablation studies on some of the techniques presented in Section 4 to assert their effectiveness. Next, we will test different sampling strategies from Section 3.2 with an increasing number of collocation points. Finally, we will review the Energy Natural Gradient Descent optimizer from Section 5.2 on two different cases. In these experiments, we will use the L² relative error with the analytical solution as the metric for success (see Eq. 2 for the definition).

## 6.1 Architecture

We have presented several techniques in Section 4. We can test their efficacy by doing ablation studies. In such an experimental setting, we compare the evolution of the L² relative error during training for multiple cases. First, we plot the error in the case where no techniques are used, a basic PINN. Second, we look at the error in the case where all techniques are used. Then, we verify what happens to the error when we remove one technique at a time.

We do an ablation study on the flat interface case of Section 2.3 and we test the following techniques: domain decomposition, Fourier embedding (embedding dimension:64, standard deviation: 4.0) and SIREN architecture (instead of FCNN). We give details on the model and training in Table 3 on the physical problem in Table 4 and on the number of collocation points in Table 5.

<table><caption>Table 3: Configuration for the model and training for the flat interface ablation study.</caption><thead><tr><th>Neurons in hidden layers</th><th>Epochs with Adam</th><th>Learning rate</th><th>Epochs with L-BFGS</th></tr></thead><tbody><tr><td>[32, 32, 3]</td><td>1000</td><td>10<sup>-4</sup></td><td>3000</td></tr></tbody></table>

<table><caption>Table 4: Configuration for the physical problem and the incoming gaussian wave for the flat interface ablation study.</caption><thead><tr><th>ε<sub>r,1</sub></th><th>ε<sub>r,2</sub></th><th>Size</th><th>Time</th><th>Frequency (wave)</th><th>Spatial width (wave)</th></tr></thead><tbody><tr><td>1.0</td><td>8.0</td><td>1 μm</td><td>15 fs</td><td>550 THz</td><td>0.5 μm</td></tr></tbody></table>

<table><caption>Table 5: Number of collocation points for each residual for the flat interface ablation study.</caption><thead><tr><th>PDEs</th><th>Left boundary</th><th>Right boundary</th><th>Interface</th><th>Initial</th><th>Periodic boundary</th></tr></thead><tbody><tr><td>8000</td><td>800</td><td>800</td><td>800</td><td>800</td><td>800</td></tr></tbody></table>

The results are shown on Figure 10. We can see that using Fourier embedding degrades the results in that case, possibly because we did not set the right parameters for it. This shows that using this technique is not always a net positive. We do a second ablation study without Fourier embedding to be more thorough on the other techniques in Figure 11. In this second ablation
<figure><img src="image_12.png" /><figcaption>Figure 10: Ablation study performed on the flat interface case: relative error during training for multiple cases (DD: domain decomposition, FE: Fourier embedding).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11: Revised ablation study performed on the flat interface case: relative error during training for multiple cases (DD: domain decomposition).</figcaption></figure>
study, we can conclude that both the SIREN architecture and the domain decomposition are invaluable to converge to a good solution. In theory, SIREN helps learning the high frequencies in the solution while domain decomposition helps to handle the discontinuity at the interface. Then, we test the other techniques in a second ablation study. This time, it is performed on the dielectric cylinder case of Section 2.5 and we test the following techniques: Fourier embedding, enforcing periodicity on time and solving the scattering problem. We give details on the model and training in Table 6, the physical problem in Table 7 and on the number of collocation points in Table 8.

<table><caption>Table 6: Configuration for the model and training for the dielectric cylinder ablation study.</caption><thead><tr><th>Neurons in hidden layers</th><th>Model</th><th>Domain decomposition</th><th>Epochs with Adam</th><th>Learning rate</th><th>Epochs with L-BFGS</th></tr></thead><tbody><tr><td>[64, 64, 3]</td><td>SIREN</td><td>Yes</td><td>2000</td><td>10<sup>-4</sup></td><td>5000</td></tr></tbody></table>

<table><caption>Table 7: Configuration for the physical problem and the incoming monochromatic wave for the dielectric cylinder ablation study/sampling comparison.</caption><thead><tr><th>ε<sub>r,1</sub></th><th>ε<sub>r,2</sub></th><th>Size</th><th>Radius</th><th>Time</th><th>Frequency (wave)</th></tr></thead><tbody><tr><td>1.0</td><td>2.0</td><td>1 μm</td><td>0.1 μm</td><td>15 fs</td><td>550 THz</td></tr></tbody></table>

<table><caption>Table 8: Number of collocation points for each residual for the dielectric cylinder ablation study.</caption><thead><tr><th>PDEs</th><th>Boundary</th><th>Interface</th></tr></thead><tbody><tr><td>16000</td><td>1000</td><td>1000</td></tr></tbody></table>

The results are shown on Figure 12. We see that while solving the scattering problem is neither very helpful or harmful, the other two techniques are very important to our success. Hence, Fourier embedding can be a boon if correctly tuned and enforcing periodicity on time is a very good idea.

## 6.2 Resampling

We will review the effects of the resampling strategies on the square cavity case of Section 2.1. To do so, we will keep the number of collocation points for the boundary and initial conditions fixed at 400 each while we change the number of collocation points for the Maxwell equations residuals as well as the sampling strategy. We have three strategies: no resampling, uniform resampling and RAD (Residual Adaptive Distribution). We compare the final relative error after the training is done.

We give details on the model and training in Table 9 and on the physical problem in Table 10. The resampling period is also kept at 300. Thus, every 300 epochs, we resample the collocation points. We give the numbers of collocation points in Table 11.
<figure><img src="image_14.png" /><figcaption>Figure 12: Ablation study performed on the dielectric cylinder case: relative error during training for multiple cases (FE: Fourier embedding, PT: enforcing periodicity in time, SP: solving the scattering problem).</figcaption></figure>

<table><caption>Table 9: Configuration for the model and training for the resampling comparison on the cavity.</caption><thead><tr><td>Neurons in hidden layers</td><td>Model</td><td>Epochs with Adam</td><td>Learning rate</td><td>Epochs with L-BFGS</td></tr></thead><tbody><tr><td>[32, 32, 3]</td><td>SIREN</td><td>3000</td><td>10<sup>-4</sup></td><td>3000</td></tr></tbody></table>

<table><caption>Table 10: Configuration for the physical problem for the resampling/optimizer comparison on the cavity.</caption><thead><tr><td>ε<sub>r</sub></td><td>Size</td><td>Time</td><td>Mode</td></tr></thead><tbody><tr><td>1.0</td><td>1 µm</td><td>5 fs</td><td>(1, 1)</td></tr></tbody></table>

<table><caption>Table 11: Number of collocation points for each residual for the resampling comparison on the cavity.</caption><thead><tr><td>PDEs</td><td>Boundary</td><td>Initial</td></tr></thead><tbody><tr><td>[5, 10, 25, 50, 100, 200, 400, 800, 1600]</td><td>400</td><td>400</td></tr></tbody></table>

The results are shown on Figure 13. In this particular case, resampling does not provide massive benefits for a too high or a too low amount of collocation points. However, the final relative error decreases faster with the number of collocation points when using resampling strategies. This could allow us to use fewer collocation points which speeds up the training process. So, on the one hand, resampling is likely to not cause any hindrance. But on the other hand, it does not seem very helpful if we use a lot of collocation points.
<figure><img src="image_15.png" /><figcaption>Figure 13: Resampling comparison on the cavity: final relative error after training depending on the number of collocation points.</figcaption></figure>

Let us apply those same sampling strategies to the dielectric cylinder case to test them. We give details on the model and training in Table 12 and on the physical problem in Table 7. The resampling period is still kept at 300. We give the numbers of collocation points in Table 13.

<table><caption>Table 12: Configuration for the model and training for resampling comparison on the dielectric cylinder.</caption><thead><tr><th>Neurons in hidden layers</th><th>Model</th><th>Domain decomposition</th><th>Fourier embedding dimension</th><th>Enforcing time period</th><th>Epochs with Adam</th><th>Learning rate</th><th>Epochs with L-BFGS</th></tr></thead><tbody><tr><td>[64, 64, 3]</td><td>SIREN</td><td>Yes</td><td>32</td><td>Yes</td><td>2000</td><td>10<sup>-4</sup></td><td>5000</td></tr></tbody></table>

<table><caption>Table 13: Number of collocation points for each residual for the resampling comparison on the dielectric cylinder.</caption><thead><tr><th>PDEs</th><th>Boundary</th><th>Interface</th></tr></thead><tbody><tr><td>[5, 10, 25, 50, 100, 200, 400, 800, 1600, 3200, 6400, 12800, 24000]</td><td>1000</td><td>1000</td></tr></tbody></table>

The results are shown on Figure 14. Once again, we see the same trends as in the cavity case but this time, the resampling strategies have even less of a lead, and perform very similarly with the no resampling case.
<figure><img src="image_16.png" /><figcaption>Figure 14: Resampling comparison on the dielectric cylinder: final relative error after training depending on the number of collocation points.</figcaption></figure>

## 6.3 Optimizer

We will test three different optimizers on the square cavity case of Section 2.1. Those optimizers are the classic Adam and L-BFGS methods as well as Energy Natural Gradient Descent (ENGD) which is more novel. First, we will work with a single-layered neural network which works to the strength of ENGD. Second, we will try to solve the same problem with a two-layered neural network.

We give details on the single-layered model and training in Table 14 and on the physical problem in Table 10. We give the numbers of collocation points in Table 15.

<table><caption>Table 14: Configuration for the single-layered model and training for the optimizer comparison on the cavity.</caption><thead><tr><th>Neurons in hidden layers</th><th>Model</th><th>Epochs with Adam</th><th>Learning rate</th><th>Epochs with L-BFGS</th><th>Epochs with ENGD</th></tr></thead><tbody><tr><td>[32, 3]</td><td>SIREN</td><td>6000</td><td>10<sup>-4</sup></td><td>3000</td><td>500</td></tr></tbody></table>

<table><caption>Table 15: Number of collocation points for each residual for the optimizer comparison on the cavity.</caption><thead><tr><th>PDEs</th><th>Boundary</th><th>Initial</th></tr></thead><tbody><tr><td>4000</td><td>400</td><td>400</td></tr></tbody></table>
<figure><img src="image_17.png" /><figcaption>Figure 15: Optimizer comparison on the cavity with a single-layered model: relative error during training for three different optimizers.</figcaption></figure>

<table><caption>Table 16: Execution times for each optimizer for the comparison on the cavity with a single-layered model (6000 epochs for Adam, 3000 epochs for L-BFGS and 500 epochs for ENGD).</caption><thead><tr><th>Adam</th><th>L-BFGS</th><th>ENGD</th></tr></thead><tbody><tr><td>43 s</td><td>45 s</td><td>51 s</td></tr></tbody></table>

The results are presented in Figure 15 and the execution times are in Table 16. We can see that the total times are pretty similar but the results vary greatly. Adam descends slowly to around $10^{-2}$ relative error. In just 1000 epochs, L-BFGS has converged to almost $10^{-4}$ relative error. Finally, in less than 20 epochs, ENGD almost got the relative error down to machine precision at $10^{-12}$. We can see here the huge advantage that second-order methods have over first-order methods when they are applied to PINNs.

Let us now consider a network with two layers instead of one. We give details on the two-layered model and training in Table 17.

<table><caption>Table 17: Configuration for the single-layered model and training for the optimizer comparison on the cavity.</caption><thead><tr><th>Neurons in hidden layers</th><th>Model</th><th>Epochs with Adam</th><th>Learning rate</th><th>Epochs with L-BFGS</th><th>Epochs with ENGD</th></tr></thead><tbody><tr><td>[32, 32, 3]</td><td>SIREN</td><td>6000</td><td>$10^{-4}$</td><td>3000</td><td>500</td></tr></tbody></table>

<table><caption>Table 18: Execution times for each optimizer for the comparison on the cavity with a two-layered model (6000 epochs for Adam, 3000 epochs for L-BFGS and 500 epochs for ENGD).</caption><thead><tr><th>Adam</th><th>L-BFGS</th><th>ENGD</th></tr></thead><tbody><tr><td>45 s</td><td>55 s</td><td>363 s</td></tr></tbody></table>
<figure><img src="image_18.png" /><figcaption>Figure 16: Optimizer comparison on the cavity with a two-layered model: relative error during training for three different optimizers.</figcaption></figure>

The results are presented in Figure 16 and the execution times are in Table 18. Here, we see that the optimizers struggled with the increase in the number of parameters resulting in slower convergence and worse results. ENGD is hit the hardest by this. It performs the best in the end but its execution time is severely degraded since the number of parameters went from $3 \times (32+1) + 32 \times (3+1) = 227$ to $3 \times (32+1) + 32 \times (3+1) + 32 \times (3+1) = 1283$. The part of ENGD that slows us down the most is solving the linear system to find the new direction for gradient descent, the size of the system is equal to the number of parameters. Hence, while the computational cost of Adam and L-BFGS are $O(n)$ where n is the number of parameters, the computational cost of ENGD is $O(n^3)$. Hence, using ENGD for a model with more than a single hidden layer is a difficult task.

## 6.4 Results for each case

In this subsection, we present some of the best models we got for each case. Every model has a SIREN architecture and uses domain decomposition if the medium is heterogeneous.

### 6.4.1 C1: Stationary wave in a metallic cavity

The model's configuration is given in table 19, the physical problem parameters are in table 20 and the number of points sampled for each residual is in table 21. We show the evolution of the loss functions during training in figure 17 as well as the relative error during training in figure 18. The training times with each optimizer are given in table 22. Finally, we show a snapshot of the electromagnetic fields in table 23.

<table><caption>Table 19: Configuration for the model and training for the cavity case (C1).</caption><thead><tr><th>Neurons in hidden layers</th><th>Model</th><th>Epochs with Adam</th><th>Learning rate</th><th>Epochs with L-BFGS</th><th>Epochs with ENGD</th></tr></thead><tbody><tr><td>[64, 3]</td><td>SIREN</td><td>1000</td><td>10<sup>-3</sup></td><td>1000</td><td>1000</td></tr></tbody></table>
<figure><img src="image_19.png" /><figcaption>Table 23: Colormaps of the electromagnetic fields at t=2.5 fs for the cavity case (C1).</figcaption></figure>

<table><caption>Table 20: Configuration for the physical problem for the cavity case (C1).</caption><thead><tr><td>ε<sub>r</sub></td><td>Size</td><td>Time</td><td>Mode</td></tr></thead><tbody><tr><td>1.0</td><td>1 µm</td><td>5 fs</td><td>(1, 1)</td></tr></tbody></table>

<table><caption>Table 21: Number of collocation points for each residual for the cavity case (C1).</caption><thead><tr><td>PDEs</td><td>Boundary</td><td>Initial</td></tr></thead><tbody><tr><td>4000</td><td>400</td><td>400</td></tr></tbody></table>

<table><caption>Table 22: Training times for each optimizer for the cavity case (C1).</caption><thead><tr><td>Adam</td><td>L-BFGS</td><td>ENGD</td><td>Total</td></tr></thead><tbody><tr><td>10.8 s</td><td>20.8 s</td><td>114.0 s</td><td>145.6 s</td></tr></tbody></table>

The Adam and L-BFGS optimizers quickly converge to a very good solution. Since a shallow model with a single hidden layer is enough to solve this problem, we can use the ENGD optimizer to achieve extremely high precision. This case presents no difficulty to be solved.
<figure><img src="image_20.png" /><figcaption>Figure 17: The training loss (subfigure 17(a)) and test loss (subfigure 17(b)) during training for the cavity case (C1). The legend is given in subfigure 17(c).</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 18: The L² relative error during training for the cavity case (C1).</figcaption></figure>

### 6.4.2 C2: Propagation of a plane wave in a homogeneous domain

The model's configuration is given in table 24, the physical problem parameters are in table 25 and the number of points sampled for each residual is in table 26. We show the evolution of the loss functions during training in figure 19 as well as the relative error during training in figure 20. The training times with each optimizer are given in table 27. Finally, we show a snapshot of the electromagnetic fields in table 28.

<table><caption>Table 24: Configuration for the model and training for the propagation case (C2).</caption><thead><tr><th>Neurons in hidden layers</th><th>Model</th><th>Epochs with Adam</th><th>Learning rate</th><th>Epochs with L-BFGS</th><th>Epochs with ENGD</th></tr></thead><tbody><tr><td>[64, 3]</td><td>SIREN</td><td>1000</td><td>10<sup>-3</sup></td><td>1000</td><td>1000</td></tr></tbody></table>

<table><caption>Table 25: Configuration for the physical problem and the incoming gaussian wave for the propagation case (C2).</caption><thead><tr><th>ε<sub>r,1</sub></th><th>Size</th><th>Time</th><th>Frequency (wave)</th><th>Spatial width (wave)</th></tr></thead><tbody><tr><td>1.0</td><td>1 μm</td><td>15 fs</td><td>550 THz</td><td>0.5 μm</td></tr></tbody></table>

<table><caption>Table 26: Number of collocation points for each residual for the propagation case (C2).</caption><thead><tr><th>PDEs</th><th>Boundary</th><th>Initial</th></tr></thead><tbody><tr><td>4000</td><td>400</td><td>400</td></tr></tbody></table>
<figure><img src="image_22.png" /><figcaption>Figure 19: The training loss (subfigure 19(a)) and test loss (subfigure 19(b)) during training for the propagation case (C2). The legend is given in subfigure 19(c).</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 20: The L² relative error during training for the propagation case (C2).</figcaption></figure>

<table><caption>Table 27: Training times for each optimizer for the propagation case (C2).</caption><thead><tr><td>Adam</td><td>L-BFGS</td><td>ENGD</td><td>Total</td></tr></thead><tbody><tr><td>10.0 s</td><td>21.3 s</td><td>131.6 s</td><td>162.9 s</td></tr></tbody></table>

Adam and L-BFGS struggle finding a solution at first but ENGD eventually converges to a good solution. Anecdotally, Adam and L-BFGS struggle more on this case when using a single hidden layer but perform much better when using two hidden layer. Once again, the case is fairly simple.

### 6.4.3 C3: Flat interface between two dielectric materials

The model's configuration is given in table 29, the physical problem parameters are in table 30 and the number of points sampled for each residual is in table 31. We show the evolution of the loss functions during training in figure 21 as well as the relative error during training in figure 22. The training times with each optimizer are given in table 32. Finally, we show a snapshot of the electromagnetic fields in table 33.

<table><caption>Table 29: Configuration for the model and training for the flat interface case (C3).</caption><thead><tr><td>Neurons in hidden layers</td><td>Epochs with Adam</td><td>Learning rate</td><td>Epochs with L-BFGS</td></tr></thead><tbody><tr><td>[32, 32, 3]</td><td>2000</td><td>10<sup>-3</sup></td><td>4000</td></tr></tbody></table>
<figure><img src="image_24.png" /><figcaption>Table 28: Colormaps of the electromagnetic fields at t=7.5 fs for the propagation case (C2).</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Table 33: Colormaps of the electromagnetic fields at t=7.5 fs for the flat interface case (C3).</figcaption></figure>

<table><caption>Table 30: Configuration for the physical problem and the incoming gaussian wave for the flat interface case (C3).</caption><thead><tr><th>ε<sub>r,1</sub></th><th>ε<sub>r,2</sub></th><th>Size</th><th>Time</th><th>Frequency (wave)</th><th>Spatial width (wave)</th></tr></thead><tbody><tr><td>1.0</td><td>8.0</td><td>1 μm</td><td>15 fs</td><td>550 THz</td><td>0.5 μm</td></tr></tbody></table>

<table><caption>Table 31: Number of collocation points for each residual for the flat interface case (C3).</caption><thead><tr><th>PDEs</th><th>Left boundary</th><th>Right boundary</th><th>Interface</th><th>Initial</th><th>Periodic boundary</th></tr></thead><tbody><tr><td>8000</td><td>1000</td><td>1000</td><td>1000</td><td>1000</td><td>1000</td></tr></tbody></table>

<table><caption>Table 32: Training times for each optimizer for the flat interface case (C3).</caption><thead><tr><th>Adam</th><th>L-BFGS</th><th>Total</th></tr></thead><tbody><tr><td>79.4 s</td><td>280.8 s</td><td>360.2 s</td></tr></tbody></table>

This is the first case that requires the use of particular techniques to be solved, domain
<figure><img src="image_26.png" /><figcaption>Figure 21: The training loss (subfigure 21(a)) and test loss (subfigure 21(b)) during training for the flat interface case (C3). The legend is given in subfigure 21(c).</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Figure 22: The L² relative error during training for the flat interface case (C3).</figcaption></figure>

decomposition is a huge help, especially when the permittivity changes drastically across the medium. We also use a model with two hidden layers and hence we avoid using ENGD which would take too long. Even with just Adam and L-BFGS, we achieve a fairly low relative error (below one percent). This case is harder to solve than the preceding ones but the PINN still performs well.

### 6.4.4 C4: Dielectric slab

The model's configuration is given in table 34, the physical problem parameters are in table 35 and the number of points sampled for each residual is in table 36. We show the evolution of the loss functions during training in figure 23 as well as the relative error during training in figure 24. The training times with each optimizer are given in table 37. Finally, we show a snapshot of the electromagnetic fields in table 38.

<table><caption>Table 34: Configuration for the model and training for the dielectric slab case (C4).</caption><thead><tr><th>Neurons in hidden layers</th><th>Enforcing time period</th><th>Epochs with Adam</th><th>Learning rate</th><th>Epochs with L-BFGS</th></tr></thead><tbody><tr><td>[32, 32, 3]</td><td>Yes</td><td>2000</td><td>10<sup>-3</sup></td><td>4000</td></tr></tbody></table>

<table><caption>Table 35: Configuration for the physical problem and the incoming monochromatic wave for the dielectric slab case (C4).</caption><thead><tr><th>List of ε<sub>r</sub></th><th>Interface placements</th><th>Size</th><th>Time</th><th>Frequency (wave)</th></tr></thead><tbody><tr><td>[1.0, 8.0, 1.0]</td><td>[0.33 μm, 0.66 μm]</td><td>1 μm</td><td>15 fs</td><td>550 THz</td></tr></tbody></table>
<figure><img src="image_28.png" /><figcaption>Table 38: Colormaps of the electromagnetic fields at t=7.5 fs for the dielectric slab case (C4).</figcaption></figure>

<table><caption>Table 36: Number of collocation points for each residual for the dielectric slab case (C4).</caption><thead><tr><th>PDEs</th><th>Left boundary</th><th>Right boundary</th><th>Interface 1</th><th>Interface 2</th><th>Periodic boundary</th></tr></thead><tbody><tr><td>8000</td><td>1000</td><td>1000</td><td>1000</td><td>1000</td><td>1000</td></tr></tbody></table>

<table><caption>Table 37: Training times for each optimizer for the dielectric slab case (C4).</caption><thead><tr><th>Adam</th><th>L-BFGS</th><th>Total</th></tr></thead><tbody><tr><td>108.3 s</td><td>391.5 s</td><td>499.8 s</td></tr></tbody></table>

This case is very similar to the flat interface case and the results are about the same. Hence, our comments are the same as for the flat interface case (C3), domain decomposition is essential and allows us to solve even with a relatively high disparity in permittivity.

### 6.4.5 C5: Diffraction of an incident wave by a dielectric cylinder

The model's configuration is given in table 39, the physical problem parameters are in table 40 and the number of points sampled for each residual is in table 41. We show the evolution of the loss functions during training in figure 25 as well as the relative error during training in figure 26. The training times with each optimizer are given in table 42. Finally, we show a snapshot of the electromagnetic fields in table 43.
<figure><img src="image_29.png" /><figcaption>Figure 23: The training loss (subfigure 23(a)) and test loss (subfigure 23(b)) during training for the dielectric slab case (C4). The legend is given in subfigure 23(c).</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Figure 24: The L² relative error during training for the dielectric slab case (C4).</figcaption></figure>

<table><caption>Table 39: Configuration for the model and training for the dielectric cylinder case (C5).</caption><thead><tr><th>Neurons in hidden layers</th><th>Enforcing time period</th><th>Fourier embedding dimension</th><th>Epochs with Adam</th><th>Learning rate</th><th>Epochs with L-BFGS</th></tr></thead><tbody><tr><td>[64, 64, 3]</td><td>Yes</td><td>32</td><td>2000</td><td>10<sup>-4</sup></td><td>5000</td></tr></tbody></table>

<table><caption>Table 40: Configuration for the physical problem and the incoming monochromatic wave for the dielectric cylinder case (C5).</caption><thead><tr><th>ε<sub>r,1</sub></th><th>ε<sub>r,2</sub></th><th>Size</th><th>Radius</th><th>Time</th><th>Frequency (wave)</th></tr></thead><tbody><tr><td>1.0</td><td>2.0</td><td>1 μm</td><td>0.1 μm</td><td>15 fs</td><td>550 THz</td></tr></tbody></table>

<table><caption>Table 41: Number of collocation points for each residual for the dielectric cylinder case (C5).</caption><thead><tr><th>PDEs</th><th>Boundary</th><th>Interface</th></tr></thead><tbody><tr><td>16000</td><td>1000</td><td>1000</td></tr></tbody></table>

<table><caption>Table 42: Training times for each optimizer for the flat interface case (C3).</caption><thead><tr><th>Adam</th><th>L-BFGS</th><th>Total</th></tr></thead><tbody><tr><td>192.6 s</td><td>927.6 s</td><td>1120.2 s</td></tr></tbody></table>

This is by far and away the most troublesome case to handle. It requires careful tuning of the Fourier embedding, as well as a lot of sampled points to approach a good solution. We often
<figure><img src="image_31.png" /><figcaption>Figure 25: The training loss (subfigure 25(a)) and test loss (subfigure 25(b)) during training for the dielectric slab cylinder (C5). The legend is given in subfigure 25(c).</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Figure 26: The L² relative error during training for the dielectric cylinder case (C5).</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Table 43: Colormaps of the electromagnetic fields at t=7.5 fs for the dielectric cylinder case (C5).</figcaption></figure>
fail to converge to a good solution and even when we do, that solution is not as good as what we had in other cases. It is difficult to know the exact reason for this failure. Do the imperfect absorbing boundary conditions generate additional waves? Or is the training at fault? Maybe the neural network itself is not well-suited and needs more layers? This case is still not solved well despite our progress.

## 6.5 Discussion

From the previous parts, we can conclude that the PINN is very sensitive to choices in architecture. Many configurations can be tested and the efficacy of each technique varies case by case. However, some ideas are never truly bad. First, using domain decomposition to handle heterogeneities seems almost mandatory in the cases where that applies. Second, using either SIREN or Fourier embedding to prevent spectral bias is also very helpful and produces much faster convergence to better models.

However, for sampling strategies, it seems we can simply ignore them if we use enough collocation points. Resampling produces better results only for some kind of middle ground between too few and too many collocation points. Otherwise, it has only been noteworthy for potentially lowering training times by letting us use fewer points.

Finally, the optimizer presents a difficult choice. Second-order gradient descent methods are enticing because they offer a way out of the bad preconditionning problems that come from using partial differential equations as loss functions, which is why PINNs are often trained using Adam for a bit before switching to L-BFGS in the literature. This works somewhat well though it eventually also has its limits. Energy Natural Gradient Descent has a lot of theoretical backing way beyond what is explained in Section 5. However, it is too slow to work outside of very small scale neural networks because the computational cost of solving a linear system scales terribly with the number of parameters which leaves us at an impasse for more difficult physical problems that require much more complex models than just a single-layered neural network. There may be potential but additional work will be needed.
# 7 Conclusions and future developments

In this work, we have established a framework to solve Maxwell equations using PINNs in the two dimensional case. In this framework, we have included a scaling which allows us to work with physical values making us closer to real applications.

This was enough to obtain good results on simple cases such as the resonant cavity for low resonant modes or the transmission of a wave through flat interfaces with low enough disparity in permittivity. However, when we complicated those cases or used a more complex geometry, we fell short and the PINN could not converge to a good solution.

Thanks to the additional techniques we introduced, we could correct this problem. We managed to simulate more difficult situations where the base formulation was insufficient. However, the case of the diffraction by the cylinder still remained difficult despite showing progress with the techniques.

In the end, PINNs show promising results in modeling the Maxwell equations though some kinks still need to be ironed out. Now, we must introduce parametrization to the PINNs so that they can gain an advantage over classical methods. Some articles show promising leads. Arthurs et al.[1] use the results of simulations from numerical methods to train a parametrized PINN. The partial differential equations residuals are then used to estimate for which value of the parameters the model performs the worst. Then, we do a new simulation for this value of the parameters which we can use as training data for the PINN. Mezzadri et al.[9] proposes to use parametric non-uniform rational B-Spline (NURBS) to encode complex shapes in a mapping from the unit square to the complex shape.
References

[1] C. J. Arthurs and A. P. King. Active training of physics-informed neural networks to aggregate and interpolate parametric solutions to the navier-stokes equations. *Journal of Computational Physics*, 438:110364, Aug. 2021. ISSN 0021-9991. URL http://dx.doi.org/10.1016/j.jcp.2021.110364.

[2] M. Baldan, P. Di Barba, and D. A. Lowther. Physics-informed neural networks for inverse electromagnetic problems. *IEEE Transactions on Magnetics*, 59(5):1-5, 2023. doi: 10.1109/TMAG.2023.3247023.

[3] Y. Chen, L. Lu, G. E. Karniadakis, and L. Dal Negro. Physics-informed neural networks for inverse problems in nano-optics and metamaterials. *Optics Express*, 28(8):11618, Apr. 2020. ISSN 1094-4087. URL http://dx.doi.org/10.1364/OE.384875.

[4] S. Dong and N. Ni. A method for representing periodic functions and enforcing exactly periodic boundary conditions with deep neural networks. *Journal of Computational Physics*, 435:110242, June 2021. ISSN 0021-9991. URL http://dx.doi.org/10.1016/j.jcp.2021.110242.

[5] M. Duruflé. *Intégration numérique et éléments finis d'ordre élevé appliqués aux équations de Maxwell en régime harmonique*. Thèse, ENSTA ParisTech, Feb. 2006. URL https://pastel.hal.science/tel-00068590.

[6] M. Frank and P. Wolfe. An algorithm for quadratic programming. *Naval Research Logistics Quarterly*, 3(1-2):95-110, 1956. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nav.3800030109.

[7] X. Glorot and Y. Bengio. Understanding the difficulty of training deep feedforward neural networks. In Y. W. Teh and M. Titterington, editors, *Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics*, volume 9 of *Proceedings of Machine Learning Research*, pages 249-256, Chia Laguna Resort, Sardinia, Italy, 13-15 May 2010. PMLR. URL https://proceedings.mlr.press/v9/glorot10a.html.

[8] J. D. Jackson. *Classical electrodynamics*. Wiley, New York, NY, 3rd edition, 1999. ISBN 9780471309321. URL http://cdsweb.cern.ch/record/490457.

[9] F. Mezzadri, J. Gasick, and X. Qian. A framework for physics-informed deep learning over freeform domains. *Computer-Aided Design*, 160:103520, 2023. ISSN 0010-4485. URL https://www.sciencedirect.com/science/article/pii/S0010448523000520.

[10] W. J. Morokoff and R. E. Caflisch. Quasi-monte carlo integration. *Journal of Computational Physics*, 122(2):218-230, 1995. URL https://www.sciencedirect.com/science/article/pii/S0021999185712090.

[11] J. Müller and M. Zeinhofer. Achieving high accuracy with pinns via energy natural gradients, 2023. URL https://arxiv.org/abs/2302.13163.

[12] M. Raissi, P. Perdikaris, and G. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. *Journal of Computational Physics*, 378:686-707, 2019. ISSN 0021-9991. URL https://www.sciencedirect.com/science/article/pii/S0021999118307125.
[13] S. Rezaei, A. Harandi, A. Moeineddin, B.-X. Xu, and S. Reese. A mixed formula-
tion for physics-informed neural networks as a potential solver for engineering problems
in heterogeneous domains: Comparison with finite element method. *Computer Methods
in Applied Mechanics and Engineering*, 401:115616, Nov. 2022. ISSN 0045-7825. URL
http://dx.doi.org/10.1016/j.cma.2022.115616.

[14] T. D. Ryck, F. Bonnet, S. Mishra, and E. de Bézenac. An operator preconditioning perspective on training in physics-informed machine learning. *ArXiv*, abs/2310.05801, 2023. URL https://api.semanticscholar.org/CorpusID:263831268.

[15] V. Sitzmann, J. Martel, A. Bergman, D. Lindell, and G. Wetzstein. Implicit neural representations with periodic activation functions. In H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, editors, *Advances in Neural Information Processing Systems*, volume 33, pages 7462–7473. Curran Associates, Inc., 2020. URL https://proceedings.neurips.cc/paper_files/paper/2020/file/53c04118df112c13a8c34b38343b9c10-Paper.pdf.

[16] N. Sukumar and A. Srivastava. Exact imposition of boundary conditions with distance func-
tions in physics-informed deep neural networks. *Computer Methods in Applied Mechanics
and Engineering*, 389:114333, Feb. 2022. ISSN 0045-7825. URL http://dx.doi.org/10.
1016/j.cma.2021.114333.

[17] M. Tancik, P. Srinivasan, B. Mildenhall, S. Fridovich-Keil, N. Raghavan, U. Singhal, R. Rammamoorthi, J. Barron, and R. Ng. Fourier features let networks learn high frequency functions in low dimensional domains. In H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, editors, *Advances in Neural Information Processing Systems*, volume 33, pages 7537–7547. Curran Associates, Inc., 2020. URL https://proceedings.neurips.cc/paper_files/paper/2020/file/55053683268957697aa39fba6f231c68-Paper.pdf.

[18] J. D. Toscano, V. Oommen, A. J. Varghese, Z. Zou, N. A. Daryakenari, C. Wu, and G. E. Karniadakis. From pinns to pikans: Recent advances in physics-informed machine learning, 2024. URL https://arxiv.org/abs/2410.13228.

[19] J. Viquerat. *Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method*. Thèse, Université Nice Sophia Antipolis, Dec. 2015. URL https://hal.science/tel-01272010.

[20] S. Wang, X. Yu, and P. Perdikaris. When and why pinns fail to train: A neural tangent kernel perspective, 2020. URL https://arxiv.org/abs/2007.14527.

[21] S. Wang, S. Sankaran, H. Wang, and P. Perdikaris. An expert's guide to training physics-informed neural networks, 2023. URL https://arxiv.org/abs/2308.08468.

[22] C. Wu, M. Zhu, Q. Tan, Y. Kartha, and L. Lu. A comprehensive study of non-adaptive and residual-based adaptive sampling for physics-informed neural networks. *Computer Methods in Applied Mechanics and Engineering*, 403:115671, Jan. 2023. ISSN 0045-7825. URL http://dx.doi.org/10.1016/j.cma.2022.115671.

A Analytical solutions

Here, we demonstrate the analytical solutions for some of the cases considered in the report.
## A.1 Flat interface

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

Let us say that the incoming wave is plane transverse magnetic with normal incidence relative to the interface. Then, so are the reflected and transmitted waves. The incident wave follows the system:

$$ egin{cases} H_{x,i} = 0 \\ H_{y,i} = -\frac{1}{Z_1} f_i(c_1 t - x) \\ E_{z,i} = f_i(c_1 t - x) \end{cases} $$

In the same way, we have the transmitted wave which propagates in medium 2:

$$ egin{cases} H_{x,t} = 0 \\ H_{y,t} = -\frac{1}{Z_2} f_t(c_2 t - x) \\ E_{z,t} = f_t(c_2 t - x) \end{cases} $$

And the reflected wave which propagates in the opposite direction:

$$ egin{cases} H_{x,r} = 0 \\ H_{y,r} = \frac{1}{Z_1} f_r(c_1 t + x) \\ E_{z,r} = f_r(c_1 t + x) \end{cases} $$

where the functions $f_i$, $f_t$ and $f_r$ are differentiable.

Now, we want to find what are the reflected and transmitted wave. To do so, we use the interface conditions:

$$ egin{cases} H_{y,i} + H_{y,r} = H_{y,t} \\ E_{z,i} + E_{z,r} = E_{z,t} \end{cases} $$

Let us say the interface is at $x_{int}$, then, the system becomes:

$$ egin{cases} -\frac{1}{Z_1} f_i(c_1 t - x_{int}) + \frac{1}{Z_1} f_r(c_1 t + x_{int}) = -\frac{1}{Z_2} f_t(c_2 t - x_{int}) \\ f_i(c_1 t - x_{int}) + f_r(c_1 t + x_{int}) = f_t(c_2 t - x_{int}) \end{cases} $$
We define the refractive index of a medium i as $n_i = \sqrt{\epsilon_i}$. We now have the following system:

$$ egin{cases} n_1 f_i(c_1 t - x_{int}) - n_1 f_r(c_1 t + x_{int}) = n_2 f_t(c_2 t - x_{int}) \\ f_i(c_1 t - x_{int}) + f_r(c_1 t + x_{int}) = f_t(c_2 t - x_{int}) \end{cases} $$

This leads to:

$$ egin{cases} f_t(c_2 t - x_{int}) = \frac{2n_1}{n_1 + n_2} f_i(c_1 t - x_{int}) \\ f_r(c_1 t + x_{int}) = \frac{n_1 - n_2}{n_1 + n_2} f_i(c_1 t - x_{int}) \end{cases} $$

The system is true for any time $t$. We can use two changes of variables: $u = c_2 t - x_{int}$ and $v = c_1 t + x_{int}$:

$$ egin{cases} f_t(u) = \frac{2n_1}{n_1 + n_2} f_i \left( \frac{c_1}{c_2} u + \frac{c_1 - c_2}{c_2} x_{int} \right) \\ f_r(v) = \frac{n_1 - n_2}{n_1 + n_2} f_i(v - 2x_{int}) \end{cases} $$

This system is true for any $u$ and $v$, we do another change of variable so that $u = c_2 t - x$ and $v = c_1 t + x$. We finally get the analytical solution:

$$ egin{cases} f_t(c_2 t - x) = \frac{2n_1}{n_1 + n_2} f_i \left( c_1 t - \frac{c_1}{c_2} x + \frac{c_1 - c_2}{c_2} x_{int} \right) \\ f_r(c_1 t + x) = \frac{n_1 - n_2}{n_1 + n_2} f_i(c_1 t + x - 2x_{int}) \end{cases} $$

## A.2 Dielectric slab

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

We establish the solution for monochromatic waves with normal incidence. We will work with the complex notations. In total, there are 5 waves of interest:

* The incident wave in medium 1: $\bar{E}_{i,1} = \bar{A}_{i,1} \exp(i\omega \left[ t - \frac{x}{c_1} \right])$

* The reflected wave in medium 1: $\bar{E}_{r,1} = \bar{A}_{r,1} \exp(i\omega \left[ t + \frac{x}{c_1} \right])$
• The transmitted wave in medium 2: $\bar{E}_{t,2} = \bar{A}_{t,2} \exp(i\omega \left[t - \frac{x}{c_2}\right])$

• The reflected wave in medium 2: $\bar{E}_{r,2} = \bar{A}_{r,2} \exp(i\omega \left[t + \frac{x}{c_2}\right])$

• The transmitted wave in medium 3: $\bar{E}_{t,3} = \bar{A}_{t,3} \exp(i\omega \left[t - \frac{x}{c_3}\right])$

Where the $\bar{A}$'s are constant complex numbers.

The interfaces are respectively on $x_1$ and $x_2$. First, let us consider the interface condition at the second interface on $x_2$, we have:

$$ \begin{cases} n_2 \bar{E}_{t,2} - n_2 \bar{E}_{r,2} = n_3 \bar{E}_{t,3} \\ \bar{E}_{t,2} + \bar{E}_{r,2} = \bar{E}_{t,3} \end{cases} $$

With the complex notations, we have:

$$ \begin{cases} n_2 \bar{A}_{t,2} \exp\left(i\omega \left[t - \frac{x_2}{c_2}\right]\right) - n_2 \bar{A}_{r,2} \exp\left(i\omega \left[t + \frac{x_2}{c_2}\right]\right) = n_3 \bar{A}_{t,3} \exp\left(i\omega \left[t - \frac{x_2}{c_3}\right]\right) \\ \bar{A}_{t,2} \exp\left(i\omega \left[t - \frac{x_2}{c_2}\right]\right) + \bar{A}_{r,2} \exp\left(i\omega \left[t + \frac{x_2}{c_2}\right]\right) = \bar{A}_{t,3} \exp\left(i\omega \left[t - \frac{x_2}{c_3}\right]\right) \end{cases} $$

This leads to:

$$ \begin{cases} 2n_2 \bar{A}_{t,2} \exp\left(i\omega \left[t - \frac{x_2}{c_2}\right]\right) = (n_2 + n_3) \bar{A}_{t,3} \exp\left(i\omega \left[t - \frac{x_2}{c_3}\right]\right) \\ (n_2 - n_3) \bar{A}_{t,2} \exp\left(i\omega \left[t - \frac{x_2}{c_2}\right]\right) = (n_2 + n_3) \bar{A}_{r,2} \exp\left(i\omega \left[t + \frac{x_2}{c_2}\right]\right) \end{cases} $$

Hence, we can define $\bar{A}_{t,3}$ and $\bar{A}_{r,2}$ in function of $\bar{A}_{t,2}$:

$$ \begin{cases} \bar{A}_{t,3} = \frac{2n_2}{n_2 + n_3} \exp\left(i\omega \left[\frac{x_2}{c_3} - \frac{x_2}{c_2}\right]\right) \bar{A}_{t,2} \\ \bar{A}_{r,2} = \frac{n_2 - n_3}{n_2 + n_3} \exp\left(-i\omega \frac{2x_2}{c_2}\right) \bar{A}_{t,2} \end{cases} $$

We define the transmission and reflection coefficient as: $T_{2\to3} = \frac{2n_2}{n_2 + n_3}$ and $R_{2\to3} = \frac{n_2 - n_3}{n_2 + n_3}$

Now, let us write the interface condition on the first interface on $x_1$:

$$ \begin{cases} n_1 \bar{E}_{i,1} - n_1 \bar{E}_{r,1} = n_2 \bar{E}_{t,2} - n_2 \bar{E}_{r,2} \\ \bar{E}_{i,1} + \bar{E}_{r,1} = \bar{E}_{t,2} + \bar{E}_{r,2} \end{cases} $$

This leads to:

$$ \begin{cases} 2n_1 \bar{E}_{i,1} = (n_1 + n_2) \bar{E}_{t,2} + (n_1 - n_2) \bar{E}_{r,2} \\ 2n_1 \bar{E}_{r,1} = (n_1 - n_2) \bar{E}_{t,2} + (n_1 + n_2) \bar{E}_{r,2} \end{cases} $$
We now replace with the complex notations:

$$
\begin{cases}
2n_1 \bar{A}_{i,1} \exp\left(i\omega \left[t - \frac{x_1}{c_1}\right]\right) = (n_1 + n_2) \bar{A}_{t,2} \exp\left(i\omega \left[t - \frac{x_1}{c_2}\right]\right) + (n_1 - n_2) \bar{A}_{r,2} \exp\left(i\omega \left[t + \frac{x_1}{c_2}\right]\right) \\
2n_1 \bar{A}_{r,1} \exp\left(i\omega \left[t + \frac{x_1}{c_1}\right]\right) = (n_1 - n_2) \bar{A}_{t,2} \exp\left(i\omega \left[t - \frac{x_1}{c_2}\right]\right) + (n_1 + n_2) \bar{A}_{r,2} \exp\left(i\omega \left[t + \frac{x_1}{c_2}\right]\right)
\end{cases}
$$

We divide by $\exp(i\omega t)$ and write $\bar{A}_{r,2}$ in function of $\bar{A}_{t,2}$:

$$
\left\{
\begin{aligned}
2n_1 \bar{A}_{i,1} \exp\left(-i\omega \frac{x_1}{c_1}\right) &= (n_1 + n_2) \bar{A}_{t,2} \exp\left(-i\omega \frac{x_1}{c_2}\right) + (n_1 - n_2) R_{2\to3} \exp\left(-i\omega \frac{2x_2}{c_2}\right) \bar{A}_{t,2} \exp\left(i\omega \frac{x_1}{c_2}\right) \\
2n_1 \bar{A}_{r,1} \exp\left(i\omega \frac{x_1}{c_1}\right) &= (n_1 - n_2) \bar{A}_{t,2} \exp\left(-i\omega \frac{x_1}{c_2}\right) + (n_1 + n_2) R_{2\to3} \exp\left(-i\omega \frac{2x_2}{c_2}\right) \bar{A}_{t,2} \exp\left(i\omega \frac{x_1}{c_2}\right)
\end{aligned}
\right.
$$

Which leads to:

$$
\left\{
\begin{aligned}
2n_1 \exp\left(-i\omega \frac{x_1}{c_1}\right) \bar{A}_{i,1} &= \left[ (n_1 + n_2) + (n_1 - n_2) R_{2\to3} \exp\left(i\omega \left[\frac{2x_1}{c_2} - \frac{2x_2}{c_2}\right]\right) \right] \exp\left(-i\omega \frac{x_1}{c_2}\right) \bar{A}_{t,2} \\
2n_1 \exp\left(i\omega \frac{x_1}{c_1}\right) \bar{A}_{r,1} &= \left[ (n_1 - n_2) + (n_1 + n_2) R_{2\to3} \exp\left(i\omega \left[\frac{2x_1}{c_2} - \frac{2x_2}{c_2}\right]\right) \right] \exp\left(-i\omega \frac{x_1}{c_2}\right) \bar{A}_{t,2}
\end{aligned}
\right.
$$

We divide by $(n_1 + n_2) \exp(-i\omega \frac{x_1}{c_2})$:

$$
\left\{
\begin{aligned}
T_{1 \to 2} \exp \left( i\omega \left[ \frac{x_1}{c_2} - \frac{x_1}{c_1} \right] \right) \bar{A}_{i,1} &= \left[ 1 + R_{1 \to 2} R_{2 \to 3} \exp \left( i\omega \left[ \frac{2x_1}{c_2} - \frac{2x_2}{c_2} \right] \right) \right] \bar{A}_{t,2} \\
T_{1 \to 2} \exp \left( i\omega \left[ \frac{x_1}{c_2} + \frac{x_1}{c_1} \right] \right) \bar{A}_{r,1} &= \left[ R_{1 \to 2} + R_{2 \to 3} \exp \left( i\omega \left[ \frac{2x_1}{c_2} - \frac{2x_2}{c_2} \right] \right) \right] \bar{A}_{t,2}
\end{aligned}
\right.
$$

We can finally write $\bar{A}_{t,2}$ in function of $\bar{A}_{i,1}$ and all other $\bar{A}$'s in function of $\bar{A}_{t,2}$:

$$
\left\{
\begin{aligned}
\bar{A}_{t,2} &= \frac{T_{1 \to 2} \exp\left(i\omega \left[\frac{x_1}{c_2} - \frac{x_1}{c_1}\right]\right)}{1 + R_{1 \to 2} R_{2 \to 3} \exp\left(i\omega \left[\frac{2x_1}{c_2} - \frac{2x_2}{c_2}\right]\right)} \bar{A}_{i,1} \\
\bar{A}_{r,1} &= \frac{R_{1 \to 2} + R_{2 \to 3} \exp\left(i\omega \left[\frac{2x_1}{c_2} - \frac{2x_2}{c_2}\right]\right)}{T_{1 \to 2} \exp\left(i\omega \left[\frac{x_1}{c_2} + \frac{x_1}{c_1}\right]\right)} \bar{A}_{t,2} \\
\bar{A}_{t,3} &= T_{2 \to 3} \exp\left(i\omega \left[\frac{x_2}{c_3} - \frac{x_2}{c_2}\right]\right) \bar{A}_{t,2} \\
\bar{A}_{r,2} &= R_{2 \to 3} \exp\left(-i\omega \frac{2x_2}{c_2}\right) \bar{A}_{t,2}
\end{aligned}
\right.
$$
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