# Fourth Order Accurate Compact Scheme for First-Order Maxwell’s Equations

I. Versano  School of Mathematical Sciences, Tel-Aviv University, Tel-Aviv
6997801, Israel<EMAIL_ADDRESS>,  E. Turkel  School of
Mathematical Sciences, Tel-Aviv University, Tel-Aviv 6997801, Israel
<EMAIL_ADDRESS>and  S. Tsynkov  North Carolina State University, Box
8205, Raleigh, NC 27695, USA<EMAIL_ADDRESS>
###### Abstract.

We construct a compact fourth-order scheme, in space and time, for the time-
dependent Maxwell’s equations given as a first-order system on a staggered
(Yee) grid. At each time step, we update the fields by solving positive
definite second-order elliptic equations. We face the challenge of finding
compatible boundary conditions for these elliptic equations while maintaining
a compact stencil. The proposed scheme is compared computationally with a non-
compact scheme and data-driven approach.

Keywords: Compact finite differences, Maxwell’s equations, High order
accuracy.

Acknowledgments

The work supported by the US–Israel Bi-national Science Foundation (BSF) under
grant # 2020128.

##  1\. introduction

Time-dependent Maxwell’s equations have been established for over 125 years as
the theoretical representation of electromagnetic phenomena [10, 11]. However,
the accuracy requirement of high-frequency simulations remains a challenge due
to a pollution effect which has been discovered firstly for the convection
equation [8] and later on for the Helmholtz equation [1, 5]. The pollution
effect states that, for a $p$ order time accurate scheme, the quantity
$\omega^{p+1}h_{\tau}^{p}$ should remain constant in order for the error to
remain constant as the wave-number $\omega$ varies. A similar phenomenon
applies also to spatial errors [1, 5]. This effect motivates the need for
higher order schemes, i.e., larger $p$. Using straightforward central
differences, higher-order accuracy requires a larger stencil. This has two
disadvantages. Firstly, the larger the stencil, the more work may be needed to
invert a matrix with a larger bandwidth and more non-zero entries. Even more
serious are the difficulties near boundaries. A large stencil requires some
modification near the boundaries where all the points needed in the stencil
are not available. This raises questions about the efficiency and stability of
such schemes. In addition, higher-order compact methods can achieve the same
error while using fewer grid nodes, making them more efficient in terms of
both storage and CPU time than non-compact or low-order methods (see e.g., [3]
for the wave equation, [14] for Helmholtz equation, [18, 19] for Maxwell’s
equations, and references therein). The Yee scheme (also known as finite-
difference time-domain, or FDTD) was introduced by Yee [17] and remains the
common numerical method used for electromagnetic simulations [16]. Although it
is only a second-order method, it is still preferred for many applications
because it preserves important structural features of Maxwell’s equations that
other methods fail to capture. One of the main novelties presented by Yee was
the staggered-grid approach, where the electric field ${\bf E}$ and magnetic
field ${\bf H}$ do not live at the same discrete space or time locations but
at separate nodes on a staggered lattice. From the perspective of differential
forms in space-time, it becomes clear that the staggered-grid approach is more
faithful to the structure of Maxwell’s equations that are dual to one another
and hence ${\bf E},{\bf H}$ naturally live on two staggered, dual meshes [15].
The question of staggered versus non-staggered high-order schemes has been
studied in [7] (see also [16, pages 63-109]), where it has been shown that,
for a given order of accuracy, a staggered scheme is more accurate and
efficient than a non-staggered scheme.

In this paper, we present a new compact fourth-order accurate scheme, in both
space and time, for Maxwell’s equations in three dimensions using an equation-
based method on a staggered grid. At each time step, we update the solution by
solving uncoupled (positive definite) second-order elliptic equations using
the conjugate gradient method. This procedure involves a non-trivial treatment
at the boundaries on which boundary conditions are not given explicitly and
have to be deduced by the equation itself while maintaining the compact
stencil. While the development of the scheme is done in 3D, the simulations
are only two-dimensional for the sake of simplicity.

###  1.1. Maxwell’s equations

Let ${\bf E},{\bf H},\bf{D},\bf{B},\bf{J}$ be (real-valued) vector fields,
defined for $(t,{\bf x})\in[0,\infty)\times\Omega$, which take values in
$\mathbb{R}^{3}$. The Maxwell equations, in first-order differential vector
form, are given by

$\displaystyle\frac{\partial\bf{D}}{\partial t}$
$\displaystyle=\nabla\\!\times\\!\bf{H}-\bf{J},$
$\displaystyle\frac{\partial\bf{B}}{\partial t}$
$\displaystyle=-\nabla\times{\bf E}$ $\displaystyle\nabla\cdot\bf{H}$
$\displaystyle=0$ $\displaystyle\nabla\cdot{\bf E}$
$\displaystyle=\frac{\rho}{\varepsilon}$ (1.1) $\displaystyle\bf{B}$
$\displaystyle=\mu\bf{H}$ $\displaystyle\bf{D}$ $\displaystyle=\epsilon{\bf
E}$

in $\Omega$, where in addition

$\qquad\frac{\partial\rho}{\partial t}+\nabla\cdot\bf{J}=0.$

Here $\nabla\times$ denotes the curl operator with respect to ${\bf x}$. We
assume that smooth initial conditions are given for $\bf{D}$ and $\bf{B}$ at
$t=0$. We assume the following:

  1. (1)

$\mu,\varepsilon$ are positive constants, unless noted otherwise, satisfying
$\frac{1}{\mu\varepsilon}=c^{2}$.

  2. (2)

$\Omega=[0,1]^{3}$, $\vec{n}\times{\bf E}=0$ on $\partial\Omega$ ([12, Section
8]).

The identity $\nabla\cdot(\nabla\times)=0$ implies that
$\frac{\partial}{\partial t}\left(\nabla\cdot{\bf H}\right)=0$ and hence the
divergence conditions are only restrictions on the initial conditions. By
letting $\tau=ct$, $Z=\sqrt{\mu/\varepsilon}$, we obtain the following re-
scaled equations:

$\displaystyle\frac{\partial{\bf E}}{\partial\tau}=Z(\nabla\times\bf{H}-{\bf
J})\quad\mathrm{in}\quad(0,\infty)\times\Omega$ (1.2a)
$\displaystyle\frac{\partial{\bf
H}}{\partial\tau}=-\frac{1}{Z}\nabla\times{\bf
E}\quad\mathrm{in}\quad(0,\infty)\times\Omega$ (1.2b)
$\displaystyle\nabla\cdot{\bf H}(\tau,x,y,z)=0$ (1.2c)
$\displaystyle\nabla\cdot{\bf E}(\tau,x,y,z)=\frac{\rho}{\varepsilon}$ (1.2d)
$\displaystyle\frac{1}{c}\frac{\partial\rho}{\partial\tau}+\nabla\cdot\bf{J}=0.$
(1.2e)

This is supplemented by the initial conditions

${\bf E}(0,x,y,z)={\bf E}_{0},\quad{\bf H}(0,x,y,z)={\bf H}_{0}$

and the boundary conditions $\vec{n}\times{\bf E}=0$ on $\partial\Omega$.

######  Remark 1.1.

Differentiating (1.2) with respect to $\tau$, we see that ${\bf E}$ and ${\bf
H}$ satisfy the wave equation

$\displaystyle\frac{\partial^{2}{\bf E}}{\partial\tau^{2}}$
$\displaystyle=-Z\frac{\partial{\bf
J}}{\partial\tau}-\nabla\times\nabla\times{\bf E}$ (1.3)
$\displaystyle=-Z\frac{\partial{\bf J}}{\partial\tau}-\nabla(\nabla\cdot{\bf
E})+\Delta{\bf E}$ $\displaystyle=-Z\frac{\partial{\bf
J}}{\partial\tau}-\nabla(\nabla\cdot{\bf E})+\Delta{\bf E}$

and

$\frac{\partial^{2}{\bf H}}{\partial\tau^{2}}=-\nabla\times\nabla\times{\bf
H}-\nabla\times\partial_{\tau}{\bf J}=\Delta{\bf
H}-\nabla\times\partial_{\tau}{\bf J}.$ (1.4)

Here $\Delta$ denotes the Laplacian with respect to ${\bf x}$.

To determine the number of imposed boundary conditions, we need to examine the
eigenvalues for the reduced one-dimensional equations at that boundary. A
positive eigenvalue indicates a variable that is entering from outside the
domain and so needs to be specified. A negative eigenvalue indicates the
variable is determined from the inside and so cannot be imposed and instead is
determined by the interior equations. A zero eigenvalue is more ambiguous.

######  Remark 1.2.

An imposed boundary condition must contain only derivatives of a lower order
than the differential equation. Thus, for a first-order system, all boundary
conditions must be of Dirichlet type. A Neumann condition is not considered an
imposed boundary condition.

To derive the boundary conditions for ${\bf E}$ and ${\bf H}$, we analyze the
characteristics of the homogeneous counterpart to system (1.1), $\rho=0$ and
${\bf J}=0$, in the direction normal to the boundary (cf. [6, Section 9.1]).
The two divergence equations can be considered as constraints on the initial
data and do not affect the characteristics. Let $w$ be the vector
$(E_{x},E_{y},E_{z},H_{x},H_{y},H_{z})^{t}$. Then, considering only the $x$
space derivatives (without loss of generality) we arrive at

$\frac{\partial w}{\partial t}=Aw_{x}$

where

$A=\begin{pmatrix}0&0&0&0&0&0\\\ 0&0&0&0&0&\frac{1}{\epsilon}\\\
0&0&0&0&-\frac{1}{\epsilon}&0\\\ 0&0&0&0&0&0\\\ \frac{1}{\mu}&0&0&0&0&0\\\
0&-\frac{1}{\mu}&0&0&0&0\end{pmatrix}$

The eigenvalues of $A$ are $c\\!=\pm\\!\frac{1}{\epsilon\mu}$ (twice) and $0$
(twice). Thus, two boundary conditions need to be imposed, two are determined
from the interior and two are not clear. We shall see that we impose ${\bf
E}_{\parallel}=0$ (two conditions) and ${\bf H}_{\perp}=0$ (one condition).

The structure of the paper is as follows: In Section 2 we introduce the
background. In Section 3 we present our scheme, and Section 4 is devoted to
stability analysis. Section 5 is devoted to numerical examples.

##  2\. preliminaries

We consider a uniform discretization of $\Omega=[0,1]^{3}$, $\Delta
x\\!=\\!\Delta y\\!=\\!\Delta z\\!=\\!h$ and $h_{\tau}$ is the time-step. We
introduce the following notations:

  * •

${\bf x}=(x,y,z)\in\mathbb{R}^{3}$.

  * •

$r=\frac{h_{\tau}}{h}$ is the Courant-Friedrichs-Lewy (CFL) number.

  * •

${\bf E}(t,{\bf x})=(E_{x}(t,{\bf x}),E_{y}(t,{\bf x}),E_{z}(t,{\bf x}))$,
${\bf H}(t,{\bf x})=(H_{x}(t,{\bf x}),H_{y}(t,{\bf x}),H_{z}(t,{\bf x}))$.

  * •

$\Delta$ denotes the Laplacian with respect to ${\bf x}$.

  * •

$\Delta{\bf E}=(\Delta E_{x},\Delta E_{y},\Delta E_{z})$, $\Delta{\bf
H}=(\Delta H_{x},\Delta H_{y},\Delta H_{z})$.

  * •

${\bf E}^{n}$ (resp. ${\bf H}^{n+1/2}$) denotes ${\bf E}(t=nh_{\tau})$ (resp.
${\bf H}(t=(n+1/2)h_{\tau})$).

  * •

For $i,j,k\in\mathbb{N}\cup\\{0\\}$, $x_{i}=ih,y_{j}=jh,z_{k}=kh$.

  * •

$x_{i+\frac{1}{2}}=\left(i+\frac{1}{2}\right)h,y_{j+\frac{1}{2}}=\left(j+\frac{1}{2}\right)h,z_{k+\frac{1}{2}}=\left(k+\frac{1}{2}\right)h$.

  * •

$[N]=\\{0,1,..,N\\},\quad N\in\mathbb{N}.$

  * •

$[k,m]=\\{k,k+1,..m\\},\quad k<m\in\mathbb{N}.$

  * •

1 denotes the identity operator.

  * •

For any operator $T:\mathbb{R}^{M}\to\mathbb{R}^{N}$ ,
$\|T\|:=\sup\limits_{0\neq x\in\mathbb{R}^{N}}\frac{|Tx|}{|x|}$.

  * •

$\partial_{\tau}:=\frac{\partial}{\partial\tau}$ and
$\partial_{\tau}^{2}:=\frac{\partial^{2}}{\partial^{2}\tau}$.

  * •

$\delta_{\tau}U^{n+\frac{1}{2}}:=\frac{U^{n+1}-U^{n}}{h_{\tau}}$.

To discretize the equations, we introduce a staggered mesh in both space and
time as in the Yee scheme [17]. With this arrangement, all space derivatives
are spread over a single mesh width, and the central time and space
derivatives are centered at the same point, similar to that of the Yee scheme.
${\bf E}$ is evaluated at time $nh_{\tau}$ while $\bf{H}$ and $\bf{J}$ are
evaluated at time $(n\\!+\\!\frac{1}{2})h_{\tau}$. For the spatial
discretization, we define the following meshes:

$\displaystyle\Omega_{h}^{E_{z}}:=\\{(x_{i},y_{j},z_{k+\frac{1}{2}}),\quad(i,j,k)\in[N]^{2}\times[N-1]\\}$
(2.1)
$\displaystyle\Omega_{h}^{E_{y}}:=\\{(x_{i},y_{j+\frac{1}{2}},z_{k}),\quad(i,j,k)\in[N]\times[N-1]\times[N]\\}$
$\displaystyle\Omega_{h}^{E_{x}}:=\\{(x_{i+\frac{1}{2}},y_{j},z_{k}),\quad(i,j,k)\in[N-1]\times[N]^{2}\\}$
$\displaystyle\Omega_{h}^{H_{z}}:=\\{(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}},z_{k}),\quad(i,j,k)\in[N-1]^{2}\times[N]\\}$
$\displaystyle\Omega_{h}^{H_{y}}:=\\{(x_{i+\frac{1}{2}},y_{j},z_{k+\frac{1}{2}}),\quad(i,j,k)\in[N-1]\times[N]\times[N-1]\\}$
(2.2)
$\displaystyle\Omega_{h}^{H_{x}}:=\\{(x_{i},y_{j+\frac{1}{2}},z_{k+\frac{1}{2}}),\quad(i,j,k)\in[N]\times[N-1]^{2}\\}.$

For each $s=x,y,z$ we define the interior meshes:

$\mathrm{int}(\Omega_{h}^{E_{s}}):=(0,1)^{3}\cap\Omega_{h}^{E_{s}},\qquad\mathrm{int}(\Omega_{h}^{H_{s}}):=(0,1)^{3}\cap\Omega_{h}^{H_{s}}.$

Figure 1. The grid points of
$\Omega_{h}^{E_{z}},\Omega_{h}^{H_{x}},\Omega_{h}^{H_{y}}$ projected on the
$(x,y)$ plane in the case $N=5$.

By convention, $x_{0}=y_{0}=z_{0}=0$ and $x_{N}=y_{N}=z_{N}=1$. We denote

$\Omega_{h}^{{\bf
E}}:=\Omega_{h}^{E_{x}}\times\Omega_{h}^{E_{y}}\times\Omega_{h}^{E_{z}},\qquad\Omega_{h}^{\bf
H}:=\Omega_{h}^{H_{x}}\times\Omega_{h}^{H_{y}}\times\Omega_{h}^{H_{z}},$

and

$\mathrm{int}(\Omega_{h}^{{{\bf
E}}}):=\mathrm{int}(\Omega_{h}^{E_{x}})\times\mathrm{int}(\Omega_{h}^{E_{y}})\times\mathrm{int}(\Omega_{h}^{E_{z}}),\quad\partial\Omega_{h}^{\bf
E}:=\Omega_{h}^{\bf E}\setminus\mathrm{int}(\Omega_{h}^{{\bf E}})$
$\mathrm{int}(\Omega_{h}^{{\bf
H}}):=\mathrm{int}(\Omega_{h}^{H_{x}})\times\mathrm{int}(\Omega^{H_{y}})\times\mathrm{int}(\Omega_{h}^{H_{z}}),\quad\partial\Omega_{h}^{\bf
H}:=\Omega_{h}^{\bf H}\setminus\mathrm{int}(\Omega^{{\bf H}})$

see Figure 1 for illustration of the sets
$\Omega_{h}^{E_{z}},\Omega_{h}^{H_{x}},\Omega_{h}^{H_{y}}$ projected on the
$x,y$ plane. The discretized functions ${\bf E}_{h},{\bf H}_{h}$ are then
defined as:

${\bf E}_{h}^{n}:={\bf E}(nh_{\tau},\Omega_{h}^{\bf E}),\quad{\bf
H}_{h}^{n+\frac{1}{2}}:={\bf H}\left((n+\frac{1}{2})h_{\tau},\Omega_{h}^{\bf
H}\right).$

With this arrangement, the boundary condition $\vec{n}\times{\bf E}=0$ on
$\partial\Omega_{h}^{{\bf E}}$ implies that

$\displaystyle
E_{x}(x_{i+\frac{1}{2}},0,z_{j})=E_{x}(x_{i+\frac{1}{2}},N,z_{j})=E_{x}(x_{i+\frac{1}{2}},y_{j},0)=E_{x}(x_{i+\frac{1}{2}},y_{j},N)=0,$
$\displaystyle
E_{y}(0,y_{j+\frac{1}{2}},z_{j})=E_{y}(N,y_{j+\frac{1}{2}},z_{j})=E_{y}(x_{i},y_{j+\frac{1}{2}},0)=E_{y}(x_{i},y_{j+\frac{1}{2}},N)=0,$
(2.3) $\displaystyle
E_{z}(0,y_{j},z_{k+\frac{1}{2}})=E_{x}(N,y_{j},z_{k+\frac{1}{2}})=E_{z}(x_{i},0,z_{k+\frac{1}{2}})=E_{z}(x_{i},N,z_{k+\frac{1}{2}})=0.$

In short, $E_{x}\\!=\\!0$ when $y\\!=\\!0,1$ and $z\\!=\\!0,1$;
$E_{y}\\!=\\!0$ when $x\\!=\\!0,1$ and $z\\!=\\!0,1$;$E_{z}\\!=\\!0$ when
$x\\!=\\!0,1$ and $y\\!=\\!0,1$.

###  2.1. Boundary conditions for time derivatives

We will now derive the boundary conditions for the time derivatives of the
fields ${\bf E}$ and ${\bf H}$ rather than the fields themselves. We do so
since time marching with our compact scheme derived in Section 3 involves
solving elliptic equations for the time derivatives of the fields (see
equation (3.1a) and (3.1b)).

By (2), we have the following Dirichlet-type boundary conditions:

$\begin{gathered}\partial_{\tau}E_{x}\\!=\\!0\ \text{if}\ y\\!=\\!0,1\
\text{and}\ z\\!=\\!0,1,\\\ \partial_{\tau}E_{y}\\!=\\!0\ \text{if}\
x\\!=\\!0,1\ \text{and}\ z\\!=\\!0,1,\\\ \partial_{\tau}E_{z}\\!=\\!0\
\text{if}\ x\\!=\\!0,1\ \text{and}\ y\\!=\\!0,1.\end{gathered}$ (2.4)

Moreover, differentiating the Gauss law of electricity $\nabla\cdot{\bf
E}=\frac{\rho}{\varepsilon}$ with respect to time and using (2) we see that,
for the remaining variables the Neumann-type conditions hold:

$\displaystyle\partial_{x}\partial_{\tau}E_{x}=\dfrac{\partial_{\tau}\rho}{\varepsilon}\
\text{if}\ x=0,1,$
$\displaystyle\partial_{y}\partial_{\tau}E_{y}=\dfrac{\partial_{\tau}\rho}{\varepsilon}\
\text{if}\ y=0,1,$
$\displaystyle\partial_{z}\partial_{\tau}E_{z}=\dfrac{\partial_{\tau}\rho}{\varepsilon}\
\text{if}\ z=0,1.$

The Faraday law $\partial_{\tau}{\bf H}=-\frac{1}{Z}\nabla\times{\bf E}$
yields the following Dirichlet-type boundary conditions:

$\begin{gathered}\partial_{\tau}H_{x}=0\ \text{if}\ x=0,1,\\\
\partial_{\tau}H_{y}=0\ \text{if}\ y=0,1,\\\ \partial_{\tau}H_{z}=0\
\text{if}\ z=0,1.\end{gathered}$ (2.5)

In Table 1, we provide the values of the second-order derivatives at the
boundaries of $\Omega$ that we need for subsequent analysis and that can be
derived from equations (2.4) and (2.5) by differentiation.

|  $x=0,1$ |  $y=0,1$ |  $z=0,1$
---|---|---|---
$\partial_{\tau}^{2}E_{x}$ |  - |  0 |  0
$\partial_{\tau}^{2}E_{y}$ |  0 |  - |  0
$\partial_{\tau}^{2}E_{z}$ |  0 |  0 |  -
$\partial_{\tau}\partial_{y}H_{z}$ |  - |  - |  0
$\partial_{\tau}\partial_{z}H_{y}$ |  - |  0 |  -
$\partial_{\tau}\partial_{z}H_{x}$ |  0 |  - |  -
$\partial_{\tau}\partial_{x}H_{z}$ |  - |  - |  0
$\partial_{\tau}\partial_{x}H_{y}$ |  - |  0 |  -
$\partial_{\tau}\partial_{y}H_{x}$ |  0 |  - |  -
Table 1. Second-order derivatives in time and mixed space-time derivatives at
the boundaries of $\Omega$. The first three rows follow from (2.4). The other
six rows are an implication of (2.5).

The Ampère law $\partial_{\tau}{\bf E}=Z(\nabla\times{\bf H}-{\bf J})$ implies
that

$\displaystyle\partial_{\tau}^{2}E_{x}=Z\left(\frac{\partial
H_{z}}{\partial_{y}\partial_{\tau}}-\frac{\partial^{2}H_{y}}{\partial_{z}\partial_{\tau}}-\partial_{\tau}{\bf
J}_{x}\right),$ $\displaystyle\partial_{\tau}^{2}E_{y}=Z\left(\frac{\partial
H_{x}}{\partial_{z}\partial_{\tau}}-\frac{\partial^{2}H_{z}}{\partial_{x}\partial_{\tau}}-\partial_{\tau}{\bf
J}_{y}\right),$ $\displaystyle\partial_{\tau}^{2}E_{z}=Z\left(\frac{\partial
H_{y}}{\partial_{x}\partial_{\tau}}-\frac{\partial^{2}H_{x}}{\partial_{y}\partial_{\tau}}-\partial_{\tau}{\bf
J}_{z}\right).$

Substituting the values from Table 1 into these equations, we obtain the
following six Neumann-type boundary conditions:

$\displaystyle\partial_{y}\partial_{\tau}H_{x}=-\partial_{\tau}{\bf J}_{z}\
\text{if}\ y=0,1\quad\partial_{z}\partial_{\tau}H_{x}=\partial_{\tau}{\bf
J}_{y}\ \text{if}\ z=0,1,$
$\displaystyle\partial_{x}\partial_{\tau}H_{y}=\partial_{\tau}{\bf J}_{z}\
\text{if}\ x=0,1\quad\partial_{z}\partial_{\tau}H_{y}=-\partial_{\tau}{\bf
J}_{x}\ \text{if}\ z=0,1,$
$\displaystyle\partial_{x}\partial_{\tau}H_{z}=-\partial_{\tau}{\bf J}_{y}\
\text{if}\ x=0,1\quad\partial_{y}\partial_{\tau}H_{z}=\partial_{\tau}{\bf
J}_{x}\ \text{if}\ y=0,1.$

With these arrangements, the derivatives $\partial_{\tau}{\bf E}$ and
$\partial_{\tau}{\bf H}$ are supplied with well-defined boundary conditions,
as summarized in Table 2. In short, on $\partial\Omega$ we have
$\partial_{\tau}{\bf E}_{\parallel}=\partial_{\tau}{\bf H}_{\perp}=0$, while
$\partial_{\tau}{\bf E}_{\perp}$ and $\partial_{\tau}{\bf H}_{\parallel}$ obey
a Neumann condition. As discussed earlier, Dirichlet conditions are true
boundary conditions while Neumann conditions are derived from the interior
equations.

|  $x=0,1$ |  $y=0,1$ |  $z=0,1$
---|---|---|---
$\partial_{\tau}E_{x}$ |  Neumann |  Dirichlet |  Dirichlet
$\partial_{\tau}E_{y}$ |  Dirichlet |  Neumann |  Dirichlet
$\partial_{\tau}E_{z}$ |  Dirichlet |  Dirichlet |  Neumann
$\partial_{\tau}H_{x}$ |  Dirichlet |  Neumann |  Neumann
$\partial_{\tau}H_{y}$ |  Neumann |  Dirichlet |  Neumann
$\partial_{\tau}H_{z}$ |  Neumann |  Neumann |  Dirichlet
Table 2. Boundary conditions for the time derivatives of the Cartesian
components of ${\bf E}$ and ${\bf H}$.

##  3\. The scheme

###  3.1. Equation based derivation

We now extend the second-order accurate Yee scheme to a fourth-order accurate
scheme, in both space and time, while maintaining the compact Yee stencil. The
idea is to use a Taylor expansion in time to the next order and then replace
the resulting third-order time derivatives by space derivatives using
Maxwell’s equations.

The fourth-order Taylor expansion applied to (1.2) and combined with equations
(1.3) and (1.4) yields:

$\displaystyle Z(\nabla\times{\bf H}^{n+1/2}-{\bf
J}^{n+1/2})=\left(\frac{\partial{\bf
E}}{\partial\tau}\right)^{n+1/2}=\frac{{\bf E}^{n+1}-{\bf
E}^{n}}{h_{\tau}}-\frac{h_{\tau}^{2}}{24}\frac{\partial^{3}{\bf
E}^{n+1/2}}{\partial\tau^{3}}+O(h_{\tau}^{4})=$
$\displaystyle\delta_{\tau}{\bf
E}^{n+1/2}-\frac{h_{\tau}^{2}}{24}\partial_{\tau}\left(-Z\partial_{\tau}{\bf
J}^{n+1/2}-\nabla(\nabla\cdot{\bf E}^{n+1/2})+\Delta{\bf
E}^{n+1/2}\right)+O(h_{\tau}^{4})=$ $\displaystyle\delta_{\tau}{\bf
E}^{n+1/2}-\frac{h_{\tau}^{2}}{24}\left(-Z\partial_{\tau}^{2}{\bf
J}^{n+1/2}-\nabla(\nabla\cdot\partial_{\tau}{\bf
E}^{n+1/2})+\Delta\partial_{\tau}{\bf E}^{n+1/2}\right)+O(h_{\tau}^{4})=$
$\displaystyle\delta_{\tau}{\bf
E}^{n+1/2}-\frac{h_{\tau}^{2}}{24}\left(-Z\partial_{\tau}^{2}{\bf
J}^{n+1/2}+Z\nabla(\nabla\cdot{\bf J}^{n+1/2})+\Delta\partial_{\tau}{\bf
E}^{n+1/2}\right)+O(h_{\tau}^{4}).$

Since $\Delta\partial_{\tau}{\bf E}^{n+1/2}=\Delta\delta_{\tau}{\bf
E}^{n+1/2}+O(h_{\tau}^{2})$, we have that

$\displaystyle Z(\nabla\times{\bf H}^{n+1/2}-{\bf J}^{n+1/2})=$
$\displaystyle\delta_{\tau}{\bf
E}^{n+1/2}-\frac{h_{\tau}^{2}}{24}\left(-Z\partial_{\tau}^{2}{\bf
J}^{n+1/2}+Z\nabla(\nabla\cdot{\bf J}^{n+1/2})+\Delta\delta_{\tau}{\bf
E}^{n+1/2}\right)+O(h_{\tau}^{4}).$

Similarly, we have for ${\bf H}$:

$\displaystyle-\frac{1}{Z}\nabla\times{\bf E}^{n+1}=\frac{\partial{\bf
H}^{n+1}}{\partial\tau}=$ $\displaystyle\delta_{\tau}{\bf
H}^{n+1}-\frac{h_{\tau}^{2}}{24}\partial_{\tau}\left(\Delta{\bf
H}^{n+1}-\nabla\times\partial_{\tau}{\bf J}^{n+1}\right)+O(h_{\tau}^{4})=$
$\displaystyle\delta_{\tau}{\bf
H}^{n+1}-\frac{h_{\tau}^{2}}{24}\left(\delta_{\tau}\Delta{\bf
H}^{n+1}-\nabla\times\partial_{\tau}^{2}{\bf J}^{n+1}\right)+O(h_{\tau}^{4}).$

Thus, at each time step we obtain several uncoupled modified Helmholtz
equations which are positive definite. Note, that for more complicated
situations the equations might be coupled at the boundary.

$\displaystyle-\Delta\delta_{\tau}{\bf
E}^{n+\frac{1}{2}}+\frac{24}{h_{\tau}^{2}}\delta_{\tau}{\bf
E}^{n+\frac{1}{2}}$ $\displaystyle=Z\frac{24}{h_{\tau}^{2}}\nabla\times{\bf
H}^{n+\frac{1}{2}}+P({\bf
J})^{n+\frac{1}{2}}+O(h_{\tau}^{4}+h^{4})\quad\mathrm{in}\quad(0,1)^{3},$
(3.1a) $\displaystyle-\Delta\delta_{\tau}{\bf
H}^{n+1}+\frac{24}{h_{\tau}^{2}}\delta_{\tau}{\bf H}^{n+1}$
$\displaystyle=-\frac{1}{Z}\frac{24}{h_{\tau}^{2}}\nabla\times{\bf
E}^{n+1}-\nabla\times\partial_{\tau}^{2}{\bf
J}^{n+1}+O(h_{\tau}^{4}+h^{4})\quad\mathrm{in}\quad(0,1)^{3},$ (3.1b)

where

$P({\bf
J})^{n+\frac{1}{2}}\stackrel{{\scriptstyle\text{def}}}{{=}}Z\left(-{\bf
J}^{n+\frac{1}{2}}-\nabla(\nabla\cdot{\bf
J}^{n+\frac{1}{2}})+\frac{\partial^{2}{\bf
J}^{n+\frac{1}{2}}}{\partial\tau^{2}}\right).$ (3.2)

Equations (3.1) are supplemented by the boundary conditions for the functions
$\delta_{\tau}{\bf E}$ and $\delta_{\tau}{\bf H}$, as discussed in Section 2.

###  3.2. Modified Helmholtz equation

In a Cartesian coordinate system, each of the equations (3.1a) and (3.1b) gets
split into three scalar modified Helmholtz equations:

$-\Delta\phi+\kappa^{2}\phi=\kappa^{2}F,\ \kappa^{2}>0.$ (3.3)

We use the fourth-order accurate compact finite-difference scheme [14] for
solving elliptic equations of the type (3.3).

Consider an equally-spaced mesh of dimension $N_{x}\times N_{y}\times N_{z}$
and size $h$ in each direction. We denote by $D_{xx}$, $D_{yy}$, and $D_{zz}$
the standard second-order central-difference operators and define

$\Delta_{h}\stackrel{{\scriptstyle\text{def}}}{{=}}D_{xx}+D_{yy}+D_{zz},\quad\varUpsilon_{h}\stackrel{{\scriptstyle\text{def}}}{{=}}D_{zz}D_{xx}+D_{yy}D_{zz}+D_{xx}D_{yy}.$

Then, the fourth-order accurate scheme for (3.3) is given by [14]:

$-\left(\Delta_{h}+\frac{h^{2}}{6}\varUpsilon_{h}\right)\phi+\kappa^{2}\left(1+\frac{\kappa^{2}h^{2}}{12}\right)\phi=\kappa^{2}\left(1+\frac{\kappa^{2}h^{2}}{12}+\frac{h^{2}}{12}\Delta_{h}\right)F.$
(3.4)

If $\Delta F$ is known at the grid nodes to fourth-order accuracy, then (3.3)
can be simplified further to

$-\left(\Delta_{h}+\frac{h^{2}}{6}\varUpsilon_{h}\right)\phi+\kappa^{2}\left(1+\frac{\kappa^{2}h^{2}}{12}\right)\phi=\kappa^{2}\left(1+\frac{\kappa^{2}h^{2}}{12}\right)F+\frac{\kappa^{2}h^{2}}{12}\Delta
F.$ (3.5)

Letting $\kappa^{2}\\!=\\!\dfrac{24}{h_{\tau}^{2}}$ in (3.1) and (3.4) gives
rise to six elliptic equations of type (3.3) with $\phi$ and $F$ specified in
Table 3. Once the discretization in space has been applied (Section 3.4),
equations (3.1) lead to six positive definite linear systems to be solved
independently by conjugate gradients. Note that, it is difficult to use
multigrid since the grid dimension for different various equations is not the
same and hence not all equations can have $2^{n}$ grid nodes in every
direction.

Eq. |  $\phi$ |  F
---|---|---
1 |  $\delta_{\tau}E_{x}^{n+\frac{1}{2}}$ |  $Z(\partial_{y}H_{z}^{n+\frac{1}{2}}-\partial_{z}H_{y}^{n+\frac{1}{2}})+\frac{h_{\tau}^{2}}{24}P({\bf J}^{n+\frac{1}{2}})_{x}$
2 |  $\delta_{\tau}E_{y}^{n+\frac{1}{2}}$ |  $Z(\partial_{z}H_{x}^{n+\frac{1}{2}}-\partial_{x}H_{z}^{n+\frac{1}{2}})+\frac{h_{\tau}^{2}}{24}P({\bf J}^{n+\frac{1}{2}})_{y}$
3 |  $\delta_{\tau}E_{z}^{n+\frac{1}{2}}$ |  $Z(\partial_{x}H_{y}^{n+\frac{1}{2}}-\partial_{y}H_{x}^{n+\frac{1}{2}})+\frac{h_{\tau}^{2}}{24}P({\bf J}^{n+\frac{1}{2}})_{z}$
4 |  $\delta_{\tau}H_{x}^{n+\frac{3}{2}}$ |  $-\frac{1}{Z}(\partial_{y}E_{z}^{n+1}-\partial_{z}E_{y}^{n+1})-\frac{h_{\tau}^{2}}{24}(\nabla\times\partial_{\tau}^{2}{\bf J}^{n+1})_{x}$
5 |  $\delta_{\tau}H_{y}^{n+\frac{3}{2}}$ |  $-\frac{1}{Z}(\partial_{z}E_{x}^{n+1}-\partial_{x}E_{z}^{n+1})-\frac{h_{\tau}^{2}}{24}(\nabla\times\partial_{\tau}^{2}{\bf J}^{n+1})_{y}$
6 |  $\delta_{\tau}H_{z}^{n+\frac{3}{2}}$ |  $-\frac{1}{Z}(\partial_{x}E_{y}^{n+1}-\partial_{y}E_{x}^{n+1})-\frac{h_{\tau}^{2}}{24}(\nabla\times\partial_{\tau}^{2}{\bf J}^{n+1})_{z}$
Table 3. The expressions for $\phi$ and $F$ in (3.3) ($P({\bf J})$ is given by
(3.2)).

###  3.3. Neumann boundary conditions

We reemphasize that Neumann conditions for a first-order system are not
imposed boundary conditions. They are rather an implication of the PDEs in the
interior. Let $\Omega=[0,1]^{3}$. We define ghost nodes outside the numerical
grid so that the Neumann conditions are satisfied to fourth-order accuracy
(see Figure 2).

Figure 2.  Left: The grid nodes for $H_{x}$ projected onto the $(x,y)$ plane
in the case $N_{x}=N_{y}=5$. The unit square is bounded by the dashed line.
Solid circles show the mesh points of $\Omega_{h}^{H_{x}}$. Hollow circles
denote the ghost points induced by the Neumann boundary conditions at $y=0,1$.
Right: The grid nodes for $H_{y}$ projected onto the $(x,y)$ plane in the case
$N_{x}=N_{y}=5$. The unit square is bounded by the dashed frame. Solid stars
show the mesh points of $\Omega_{h}^{H_{y}}$. Hollow stars denote ghost points
induced by the Neumann boundary conditions at $x=0,1$.

We propose two methods for the implementation of the Neumann boundary
conditions.

#### Method 1 (equation based)

Consider, without loss of generality, equation (3.3) with the Neumann boundary
condition $\partial_{x}\phi(0,y,z)=g(y,z)$. Differentiating (3.3) with respect
to $x$ (cf. [14]), we have:

$\phi_{x}(0,x,z)=\frac{\phi(\frac{1}{2}h,y,z)-\phi(-\frac{1}{2}h,y,z)}{h}+\frac{h^{2}}{24}\left(\phi_{xyy}+\phi_{xzz}-\kappa^{2}\phi_{x}-\kappa^{2}F_{x}\right)+O(h^{4}).$
(3.6)

Assume that $\phi_{xyy}+\phi_{xzz}$ is known for $x=0$. Setting $x\\!=\\!0$ in
(3.6), we define the ghost value $\phi(-\frac{1}{2}h,y,z)$ such that the
Neumann boundary condition $\partial_{x}\phi\\!=\\!=g(y,z)$ is satisfied to
fourth-order.

#### Method 2 (Taylor based)

By Taylor’s expansion

$\phi_{x}(0,y,z)=\frac{\phi(\frac{1}{2}h,y,z)-\phi(-\frac{1}{2}h,y,z)}{h}-\frac{h^{2}}{24}\partial_{xxx}\phi(0,y,z)+O(h^{4}).$
(3.7)

Hence, knowing $\partial_{x}\phi$ and $\partial_{xxx}\phi$ at $x\\!=\\!0$
enables one to obtain the ghost variable $\phi(-\frac{1}{2}h,y,z)$ using
(3.7).

Next, we show how to approximate the Neumann-type boundary conditions to
fourth-order accuracy for $\partial_{\tau}{\bf E}$ using Method 1 and for
$\partial_{\tau}{\bf H}$ using Method 2.

######  Lemma 3.1.

Let $\phi(x,y,z)=\partial_{\tau}E_{x}$ be evaluated at some fixed $\tau_{0}$.
Then, the quantities

$\phi_{xyy}(0,y,z),\ \phi_{xzz}(0,y,z),\ \phi_{x}(0,y,z),\ \text{and}\ F_{x}$

are known to fourth-order accuracy. (Recall that, $F$ contains the terms
involving only $\partial_{z}H_{y}-\partial_{y}H_{z}$ and ${\bf J}$.) If, in
addition, ${\bf J}=0$ and $\rho=0$, then
$\phi_{xyy}(0,y,z)=\phi_{xzz}(0,y,z)=\phi_{x}(0,y,z)=F_{x}(0,y,z)=0$.

###### Proof.

$\phi_{x}(0,y,z)$ equals to $\partial_{\tau}E_{x}$ at $x=0$, which has already
been specified. Moreover, we have

$\partial_{x}E_{x}+\partial_{y}E_{y}+\partial_{z}E_{z}=\frac{\rho}{\varepsilon}.$
(3.8)

Differentiating (3.8) twice with respect to $y$ and once with respect to
$\tau$, we arrive at

$-\partial_{xyy}\phi+\partial_{\tau
yy}(\rho/\varepsilon)=\partial_{yyy}\phi+\partial_{zyy}\phi.$

If $x=0$, then $\partial_{\tau yyy}E_{x}=\partial_{\tau zyy}E_{x}=0$, and
hence if $x=0$, then $\partial_{xyy}\phi=\partial_{\tau
yy}(\rho/\varepsilon)$. Differentiating (3.8) twice with respect to $z$ and
once with respect to $\tau$, we get

$-\partial_{xzz}\phi+\partial_{\tau
zz}(\rho/\varepsilon)=\partial_{yzz}\phi+\partial_{zzz}\phi.$

If $x=0$, then $\partial_{\tau zzz}E_{x}=\partial_{\tau yzz}E_{x}=0$, and
hence at $x=0$ we have $\partial_{xzz}\phi=\partial_{\tau
zz}(\rho/\varepsilon)$.

To obtain $F_{x}$ at $x=0$, it is sufficient to evaluate
$\partial_{x}(\partial_{z}H_{y}-\partial_{y}H_{z})$. This quantity, in turn,
can be derived via the Ampère law $\partial_{\tau}{\bf E}=Z(\nabla\times{\bf
H}-{\bf J})$ if $\partial_{x\tau}E_{x}$ is available. And the latter is known
through the Neumann boundary condition for $\partial_{\tau}{\bf E}$ at $x=0$.

The second assertion follows immediately by virtue of the same argument. ∎

By repeating the previous proof for the remaining faces of the cube
$\Omega=[0,\,1]^{3}$ and field components $E_{y}$ and $E_{z}$, we obtain the
following

######  Corollary 3.2.

Neumann boundary conditions for equation (3.1a) admit a compact fourth-order
accurate discretization by means of Method 1.

Next, we will show how one can use Method 2 to build a compact discretization
of the Neumann boundary conditions for $\partial_{\tau}{\bf H}$.

######  Lemma 3.3.

Let $\phi(x,y,z)=\partial_{\tau}H_{y}$ be evaluated at some fixed $\tau_{0}$.
Then, the term $\phi_{xxx}$ at $x=0,1$ can be approximated to fourth-order
accuracy.

###### Proof.

Assume without loss of generality that $Z=1$. The equation

$\partial_{\tau}E_{z}=\partial_{x}H_{y}-\partial_{y}H_{x}-{\bf J}_{z}$ (3.9)

(cf. (1.2a)) and boundary conditions $\partial_{\tau}E_{z}=0$ at $x=0,1$ imply
that at $x=0,1$

$\partial_{\tau xxx}H_{y}=\partial_{\tau yxx}H_{x}+\partial_{\tau xx}{\bf
J}_{z}.$

Then, we use $\partial_{\tau}\nabla\cdot{\bf H}=0$ and replace
$\partial_{x}H_{x}$ with $-\partial_{\tau y}H_{y}-\partial_{\tau z}H_{z}$.
This yields:

$\partial_{\tau yxx}H_{x}=-\partial_{\tau yyx}H_{y}-\partial_{\tau xyz}H_{z}.$

To complete the proof, we have to approximate $\partial_{\tau yyx}H_{y}$ and
$\partial_{\tau xyz}H_{z}$. This is done as follows. We use (3.9) again to
obtain

$\partial_{\tau yyx}H_{y}=\partial_{\tau yyy}H_{x}+\partial_{\tau yy}{\bf
J}_{z}.$

At $x=0,1$, the first equation of (2.5) implies that $\partial_{\tau
yyy}H_{x}=0$.

To approximate $\partial_{\tau xyz}H_{z}$, we use the $y$-component of the
Ampère law (1.2a):

$\partial_{\tau}E_{y}=\partial_{z}H_{x}-\partial_{x}H_{z}-{\bf J}_{y}.$

Therefore, at $x=0,1$ we have:

$\partial_{x}H_{z}={\bf J}_{y}+\partial_{z}H_{x}$

and

$\partial_{\tau xyz}H_{z}=\partial_{\tau yz}{\bf J}_{y}+\partial_{\tau
yzz}H_{x}.$

Since $\partial_{\tau}H_{x}=0$ at $x=0,1$, we finally derive:

$\partial_{\tau yzz}H_{x}=0.$

This completes the proof. ∎

By repeating the previous procedure for $H_{x}$ and $H_{z}$, we obtain the
following corollaries.

######  Corollary 3.4.

Neumann boundary conditions for equation (3.1b) admit a compact fourth-order
accurate discretization by means of Method 2.

######  Corollary 3.5.

Assume that $\rho=0$ and ${\bf J}=0$. Then, all Neumann boundary conditions
obtained for $\partial_{\tau}{\bf E}$ and $\partial_{\tau}{\bf H}$ are
homogeneous and all the derivatives in (3.6), (3.7) vanish. In particular,
$\phi(\tilde{x}-\frac{1}{2})=\phi(\tilde{x}+\frac{1}{2})$ is a fourth-order
accurate approximation to the homogeneous Neumann condition at the boundary
point $\tilde{x}$ (see also Section 5).

We emphasize that, the right hand side, $F$, needs to be approximated with
fourth-order accuracy. This is done using a fourth-order Padé approximation
for the curl operator (see details in Appendix A.1). Then, we solve equations
1, 2, and 3 from Table 3 using the scheme (3.4), and subsequently solve
equations 4, 5, and 6 from Table 3 using the scheme (3.5).

In view of equations (3.4) and (3.5), we also define the following operators.

######  Definition 3.6.

For any $s=x,y,z$ and $G=H_{s},E_{s},$ let

$P_{1}^{G},P_{2}^{G}:\mathrm{int}(\Omega_{h})^{G}\to\mathrm{int}(\Omega_{h})^{G}$

be the symmetric operators given by

$P_{1}^{G}:=-\left(\Delta_{h}+\frac{h^{2}}{6}\varUpsilon_{h}\right)+\frac{24}{h_{\tau}^{2}}\left(1+\frac{2}{r^{2}}\right),$
$P_{2}^{G}:=\frac{24}{h_{\tau}^{2}}\left(-1-\frac{2}{r^{2}}-\frac{h^{2}}{12}\Delta_{h}\right).$

We omit the superscript $G$ whenever it does not lead to any ambiguity.

Hereafter, $r=\frac{h_{\tau}}{h}$ will denote the CFL number. As shown in
Appendix A.2, if $r<\sqrt{3+\sqrt{21}}$, then $P_{2}^{G}$ and $P_{1}^{G}$ are
symmetric positive definite matrices and the following estimates hold:

$\displaystyle 0<\frac{2}{r^{2}}\leq\|\frac{h_{\tau}^{2}}{24}P_{2}^{G}\|\leq
1+\frac{2}{r^{2}},$ (3.10) $\displaystyle
1<1+\frac{2}{r^{2}}-\frac{r^{2}}{6}\leq\|\frac{h_{\tau}^{2}}{24}P_{1}^{G}\|\leq
1+\frac{2}{r^{2}}+\frac{r^{2}}{2}.$ (3.11)

In particular,

$\|\left(P_{1}^{G}\right)^{-1}P_{2}^{G}\|\geq\frac{\frac{2}{r^{2}}}{1+\frac{2}{r^{2}}+\frac{r^{2}}{2}}.$

###  3.4. The numerical scheme

Let

${\bf
curl}_{h}\stackrel{{\scriptstyle\text{def}}}{{=}}\begin{pmatrix}0&-\delta_{z}&\delta_{y}\\\
\delta_{z}&0&-\delta_{x}\\\ -\delta_{y}&\delta_{x}&0\\\ \end{pmatrix}$

be the block matrix representing a spatial Padé fourth-order finite-difference
approximation of the curl operator. See Appendix A.1 for the full definition
of the matrix ${\bf curl}_{h}$.

Let ${\bf E}^{n}_{h}$ and $\Delta{\bf E}^{n}_{h}$ be given on
$\Omega_{h}^{{\bf E}}$ and let ${\bf H}^{n+\frac{1}{2}}_{h}$ be given on
$\Omega_{h}^{{\bf H}}$. To advance in time, we take two steps.
Step 1:

${{\bf E}}_{h}^{n+1}=\begin{pmatrix}E_{x}^{n}+h_{\tau}E_{x}^{*}\\\
E_{y}^{n}+h_{\tau}E_{y}^{*}\\\ E_{y}^{n}+h_{\tau}E_{y}^{*}\end{pmatrix}$

where

$\displaystyle\underbrace{\begin{pmatrix}P_{1}^{E_{x}}&0&0\\\
0&P_{1}^{E_{y}}&\\\ 0&0&P_{1}^{E_{z}}\end{pmatrix}}_{P_{1}^{\bf
E}}\begin{pmatrix}E_{x}^{*}\\\ E_{y}^{*}\\\ E_{z}^{*}\end{pmatrix}=ZP({\bf
J}^{n+1/2})+$
$\displaystyle\kappa^{2}Z\begin{pmatrix}0&-\delta_{z}&\delta_{y}\\\
\delta_{z}&0&-\delta_{x}\\\
-\delta_{y}&\delta_{x}&0\end{pmatrix}\left[\left(1+\frac{\kappa^{2}h^{2}}{12}\right)\begin{pmatrix}H_{x}^{n+\frac{1}{2}}\\\
H_{y}^{n+\frac{1}{2}}\\\
H_{z}^{n+\frac{1}{2}}\end{pmatrix}+\frac{h^{2}}{12}\begin{pmatrix}\Delta_{h}H_{x}^{n+\frac{1}{2}}\\\
\Delta_{h}H_{y}^{n+\frac{1}{2}}\\\
\Delta_{h}H_{z}^{n+\frac{1}{2}}\end{pmatrix}\right]$

(cf. Section 3.2). Extend ${{\bf E}}_{h}^{n+1}$ to $\Omega_{h}^{\bf E}$ using
the boundary conditions. Define

$\Delta{\bf E}_{h}^{n+1}=\Delta{\bf E}_{h}^{n}+\kappa^{2}({\bf
E}_{h}^{n+1}-{\bf
E}_{h}^{n})-Z\frac{24}{h_{\tau}}\begin{pmatrix}0&-\delta_{z}&\delta_{y}\\\
\delta_{z}&0&-\delta_{x}\\\
-\delta_{y}&\delta_{x}&0\end{pmatrix}\begin{pmatrix}H_{x}^{n+\frac{1}{2}}\\\
H_{y}^{n+\frac{1}{2}}\\\ H_{z}^{n+\frac{1}{2}}\end{pmatrix}-h_{\tau}P({\bf
J})^{n+\frac{1}{2}}$

on $\Omega_{h}^{\bf E}$ (see (3.1))
Step 2:

${{\bf
H}}_{h}^{n+\frac{3}{2}}=\begin{pmatrix}H_{x}^{n+\frac{1}{2}}+h_{\tau}H_{x}^{*}\\\
H_{y}^{n+\frac{1}{2}}+h_{\tau}H_{y}^{*}\\\
H_{y}^{n+\frac{1}{2}}+h_{\tau}H_{y}^{*}\end{pmatrix}$

where

$\displaystyle\underbrace{\begin{pmatrix}P_{1}^{H_{x}}&0&0\\\
0&P_{1}^{H_{y}}&\\\ 0&0&P_{1}^{H_{z}}\end{pmatrix}}_{P_{1}^{\bf
H}}\begin{pmatrix}H_{x}^{*}\\\ H_{y}^{*}\\\
H_{z}^{*}\end{pmatrix}-\frac{1}{Z}\nabla\times\partial_{\tau}^{2}{\bf
J}^{n+1}=$
$\displaystyle-\frac{\kappa^{2}}{Z}\begin{pmatrix}0&-\delta_{z}&\delta_{y}\\\
\delta_{z}&0&-\delta_{x}\\\
-\delta_{y}&\delta_{x}&0\end{pmatrix}\left[\left(1+\frac{\kappa^{2}h^{2}}{12}\right)\begin{pmatrix}E_{x}^{n+1}\\\
E_{y}^{n+1}\\\ E_{z}^{n+1}\end{pmatrix}+\frac{h^{2}}{12}\begin{pmatrix}\Delta
E_{x}^{n+1}\\\ \Delta E_{y}^{n+1}\\\ \Delta E_{z}^{n+1}\end{pmatrix}\right],$

(cf. Section 3.2). Extend ${\bf H}_{h}^{n+3/2}$ to $\Omega_{h}^{\bf H}$ using
the appropriate boundary condition.

##  4\. stability analysis

We assume, with no loss of generality, that $\rho\\!=\\!0$, ${\bf J}\\!=\\!0$,
and $Z\\!=\\!1$. The scheme can therefore be written in a compact way:

$\displaystyle{\bf E}_{h}^{n+1}-{\bf E}_{h}^{n}$
$\displaystyle=h_{\tau}P_{1}^{-1}{\bf
curl}_{h}\left(\kappa^{2}\left(1+\frac{\kappa^{2}h^{2}}{12}\right){\bf
H}_{h}^{n+1/2}+\frac{\kappa^{2}h^{2}}{12}\Delta{\bf H}_{h}^{n+1/2}\right),$
(4.1a) $\displaystyle{\bf H}_{h}^{n+3/2}-{\bf H}_{h}^{n+1/2}$
$\displaystyle=-h_{\tau}P_{1}^{-1}{\bf
curl}_{h}\left(\kappa^{2}\left(1+\frac{\kappa^{2}h^{2}}{12}\right){\bf
E}_{h}^{n+1}+\frac{\kappa^{2}h^{2}}{12}\Delta{\bf E}_{h}^{n+1}\right).$ (4.1b)

We approximate $\Delta$ in equations (4.1a) and (4.1b) using the standard
second-order difference Laplacian $\Delta_{h}$ (same as in (3.4)). Then, we
recast (4.1a), (4.1b) as:

$\displaystyle\left(\bf{1}-h_{\tau}\begin{pmatrix}0&0\\\ -P_{1}^{-1}P_{2}{\bf
curl}_{h}&0\end{pmatrix}\right)\begin{pmatrix}{\bf E}^{n+1}\\\ {\bf
H}^{n+3/2}\end{pmatrix}=$
$\displaystyle\left(\bf{1}-h_{\tau}\begin{pmatrix}0&-P_{1}^{-1}P_{2}{\bf
curl}_{h}\\\ 0&0\end{pmatrix}\right)\begin{pmatrix}{\bf E}^{n}\\\ {\bf
H}^{n+1/2}\end{pmatrix},$ (4.2)

where the operators $P_{1}$ and $P_{2}$ are introduced in Definition 3.6.
Next, assume the solution is in the form of a plane wave:

$\begin{pmatrix}{\bf E}^{n}\\\ {\bf
H}^{n+1/2}\end{pmatrix}_{i,j,k}=\sigma^{n}(\exp(-\sqrt{-1}(\xi_{x}ih+\xi_{y}jh+\xi_{z}kh)))\begin{pmatrix}{\bf
E}_{0}\\\ {\bf H}_{0}\end{pmatrix}.$

Let $Q=-P_{1}^{-1}P_{2}{\bf curl}_{h}$ and let $\lambda$ denote an eigenvalue
of $Q$. We substitute the plane wave solution into (4.2) and derive

$\begin{pmatrix}\sigma-1&\lambda h_{\tau}\\\ -\sigma
h_{\tau}\lambda&\sigma-1\end{pmatrix}\begin{pmatrix}{\bf E}_{0}\\\ {\bf
H}_{0}\end{pmatrix}=\begin{pmatrix}0\\\ 0\end{pmatrix}.$

Therefore,

$\sigma^{2}-\sigma(2-\lambda^{2}h_{\tau}^{2})+1=0.$

In particular, $|h_{\tau}\lambda|\leq 2$ implies that $|\sigma|\leq 1$. Then,
using the definition of the matrix $A$ given by equation (A.2) and Remark A.1
(see Appendix A.1), we can derive the stability condition in the following
form:

$r2\sqrt{3}\|A^{-1}\|\|-P_{1}^{-1}P_{2}\|\leq 2\Leftrightarrow
r\leq\frac{1}{\sqrt{3}\|A^{-1}\|\cdot\|-P_{1}^{-1}P_{2}\|}.$ (4.3)

By (3.10), $r<\sqrt{3+\sqrt{21}}$ implies that
$\|-P_{1}^{-1}P_{2}\|\geq\frac{\frac{2}{r^{2}}}{1+\frac{2}{r^{2}}+\frac{r^{2}}{2}}$
(see Appendix A.2 for detail). Therefore, the scheme is stable provided that

$r\leq\min\left\\{\sqrt{3+\sqrt{21}},\frac{1+\frac{2}{r^{2}}+\frac{r^{2}}{2}}{\sqrt{3}\|A^{-1}\|\frac{2}{r^{2}}}\right\\}.$

The inequality

$\frac{1+\frac{2}{r^{2}}+\frac{r^{2}}{2}}{\frac{2}{r^{2}}}\geq 1,\quad
r\in(0,\infty)$

gives rise to a sufficient stability condition

$0<r\leq\frac{1}{\sqrt{3}\|A^{-1}\|}\sim\frac{5}{6\sqrt{3}}.$ (4.4)

In the following section we will examine a two-dimensional reduction of
Maxwell’s equation. A modification of Remark A.1 readily implies the stability
condition

$0<r\leq\frac{1}{\sqrt{3}\|A^{-1}\|}\sim\frac{5}{6\sqrt{2}}$ (4.5)

(cf. [18]). The numerical results summarized in Table 5 corroborate the
stability estimate (4.5).

##  5\. Example: Transverse magnetic waves in $[0,1]^{2}$

###  5.1. Numerical simulations data

For the computations, we consider the scaled two-dimensional TM system without
any current or charges. Thus, the equations reduce to

$\displaystyle\frac{\partial E_{z}}{\partial\tau}$
$\displaystyle=Z\left(-\frac{\partial H_{x}}{\partial y}+\frac{\partial
H_{y}}{\partial x}\right),$ (5.1a) $\displaystyle\frac{\partial
H_{x}}{\partial\tau}$ $\displaystyle=-\frac{1}{Z}\frac{\partial
E_{z}}{\partial y},$ (5.1b) $\displaystyle\frac{\partial H_{y}}{\partial\tau}$
$\displaystyle=\frac{1}{Z}\frac{\partial E_{z}}{\partial x}.$ (5.1c)

Since ${\bf H}$ is given at time moments $n+1/2$ we need to specify its
initial condition at the time $t\\!=\\!\frac{\Delta_{\tau}}{2}$. By Taylor
series,

${\bf H}\left(\frac{h_{\tau}}{2}\right)\\!=\\!{\bf
H}(0)\\!+\\!\left(\frac{h_{\tau}}{2}\right){\bf
H}_{\tau}(0)\\!+\\!\frac{1}{2}\left(\frac{h_{\tau}}{2}\right)^{2}{\bf
H}_{\tau\tau}(0)\\!+\\!\frac{1}{6}\left(\frac{h_{\tau}}{2}\right)^{3}{\bf
H}_{\tau\tau\tau}(0)\\!+\\!\frac{1}{24}\left(\frac{h_{\tau}}{2}\right)^{4}{\bf
H}_{\tau\tau\tau\tau}(0).$ (5.2)

${\bf H}_{\tau}(0)$ is given by (1.2b) and ${\bf H}_{\tau\tau}(0)$ is given by
(1.4). Differentiating (1.4) yields:

$\frac{\partial^{3}{\bf H}}{\partial\tau^{3}}=-\nabla\times\nabla\times{\bf
H}_{\tau}-\nabla\times\partial_{\tau\tau}{\bf J}=\Delta{\bf
H}_{\tau}-\nabla\times\partial_{\tau\tau}{\bf J},$

where, again, ${\bf H}_{\tau}$ is given by (1.2b). Similarly,

$\frac{\partial^{4}{\bf H}}{\partial\tau^{4}}=-\nabla\times\nabla\times{\bf
H}_{\tau\tau}-\nabla\times\partial_{\tau^{3}}{\bf J}=\Delta{\bf
H}_{\tau\tau}-\nabla\times\partial_{\tau\tau\tau}{\bf J},$

where ${\bf H}_{\tau\tau}$ is given by (1.4). At the nodes on or next to a
boundary, we have either Dirichlet or Neumann conditions as given in Table 2.
For (5.1), we have assumed ${\bf J}\\!=\\!0$. For the computations, we further
assume $Z\\!=\\!1$.

Let $k_{x},k_{y}\in\mathbb{N}$ and let
$\omega:=Z\pi\sqrt{k_{x}^{2}+k_{y}^{2}}$. We consider the analytical solutions

$\displaystyle E_{z}(x,y)$ $\displaystyle=\cos(\omega\tau)\sin{\left(Z\pi
k_{x}x\right)}\sin{\left(Z\pi k_{y}y\right)},$ (5.3a) $\displaystyle
H_{x}(x,y)$ $\displaystyle=-\frac{\sin(\omega\tau)}{\omega}\pi k_{y}\sin(Z\pi
k_{x}x)\cos{\left(Z\pi k_{y}y\right)},$ (5.3b) $\displaystyle H_{y}(x,y)$
$\displaystyle=\frac{\sin\left(\omega\tau\right)}{\omega}\pi
k_{x}\cos{\left(Z\pi k_{x}x\right)}\sin{\left(Z\pi k_{y}y\right)}.$ (5.3c)

We use the mean absolute error

$\displaystyle\frac{1}{3N_{\tau}\cdot N^{2}}\sum_{i,j,n=0}^{N,N,N_{\tau}}$
$\displaystyle\left|{E_{z}}^{n}_{\mathrm{numer.}}(i,j)-{E_{z}}_{\mathrm{true}}^{n}(i,j)\right|+\left|{H_{x}}^{n}_{\mathrm{numer.}}(i,j)-{H_{x}}_{\mathrm{true}}^{n}(i,j)\right|+$
$\displaystyle\left|{H_{y}}^{n}_{\mathrm{numer.}}(i,j)-{H_{y}}_{\mathrm{true}}^{n}(i,j)\right|,$

where the subscripts “$\mathrm{numer.}$” and “$\mathrm{true}$” refer to the
numerical and analytical solution, respectively. The ghost nodes are easily
defined using Corollary 3.5. The ghost points for $H_{x}$ (see Figure 2) are
defined similarly.

###  5.2. Comparable schemes: $C4,NC,AI^{h}$

We compare our proposed scheme, denoted by $C4$, with the following two
schemes, which are second-order accurate in time. These schemes are obtained
according to the Yee updating rules:

$\displaystyle
E_{z}^{n+1}=E_{z}^{n}+h_{\tau}(D_{x}H_{y}^{n+1/2}-D_{y}H_{x}^{n+1/2}),$
$\displaystyle H_{x}^{n+3/2}=H_{x}^{n+1/2}-h_{\tau}D_{y}E_{z}^{n+1},$
$\displaystyle H_{y}^{n+3/2}=H_{y}^{n+1/2}+h_{\tau}D_{x}E_{z}^{n+1},$

where $D_{x}$ and $D_{y}$ are finite difference operators that approximate the
first derivative on different stencils. For the derivatives in the $x$
direction, we consider a general stencil

$S=\frac{1}{h}\begin{pmatrix}-b&-a&a&b\\\ -d&-c&c&d\\\
-b&-a&a&b\end{pmatrix},$

whereas for the derivatives in the $y$ direction we use the transposed stencil
$S^{t}$. In particular,

$\displaystyle\partial_{x}u_{ij}\sim\frac{c(u_{i+\frac{1}{2},j}-u_{i-\frac{1}{2},j})+d(u_{i+\frac{3}{2},j}-u_{i-\frac{3}{2},j})}{h}+$
$\displaystyle\frac{a(u_{i+\frac{1}{2},j+1}+u_{i+\frac{1}{2},j-1}-u_{i-\frac{1}{2},j+1}-u_{i-\frac{1}{2},j-1})+b(u_{i+\frac{3}{2},j+1}+u_{i+\frac{3}{2},j-1}-u_{i-\frac{3}{2},j+1}-u_{i-\frac{3}{2},j-1})}{h}.$

For grid nodes near the boundary, we use the standard fourth-order accurate
one-sided finite-difference approximation of the first derivatives.

By a Taylor expansion, we obtain second-order accuracy provided that

$c+3d+2a+6b=1,$

and fourth-order accuracy if the additional constraints hold:

$\displaystyle c+27d+2a+54b$ $\displaystyle=0,$ (5.4a) $\displaystyle a+3b$
$\displaystyle=0.$ (5.4b)

We define the following stencils:

$K_{2}=\frac{1}{\Delta x}\begin{pmatrix}-b&-a&a&b\\\
-d&-1+3d+2a+6b&1-3d-2a-6b&d\\\ -b&-a&a&b\end{pmatrix},$ (5.5)
$K_{4}=\frac{1}{\Delta x}\begin{pmatrix}\frac{a}{3}&-a&a&-\frac{a}{3}\\\
-\frac{16a-1}{24}&-\frac{9-16a}{8}&\frac{9-16a}{8}&\frac{16a-1}{24}\\\
\frac{a}{3}&-a&a&-\frac{a}{3}\end{pmatrix}.$ (5.6)

The non-compact fourth-order accurate scheme that we denote NC exploits the
stencil $K_{4}$ with $a\\!=\\!0$. Its order of accuracy is
$O(h^{4})+O(h_{\tau}^{2})$.

The data-driven scheme $AI^{h}$ (see Appendix A.3) uses the stencil $K_{2}$,
where the free parameters are obtained by means of a minimization process over
the given training data. This scheme is of order $O(h^{2})+O(h_{\tau}^{2})$.
We use this approach since it has been recently shown to reduce the numerical
dispersion for the wave equation [9]. The details of the minimization process
and the selected training data are detailed in Appendix A.3.

$N$ |  NC |  C4 |  $AI^{h}$
---|---|---|---
$k_{x}=k_{y}=2$ (PPW 64), $r=\frac{5}{6\sqrt{2}}$
32 |  1.94 |  5.06 |  2.69
64 |  1.98 |  4.98 |  1.89
128 |  1.98 |  4.84 |  1.88
256 |  1.99 |  4.48 |  1.93
512 |  2.00 |  4.09 |  1.96
$k_{x}=k_{y}=21$ (PPW $\sim$6), $r=\frac{5}{6\sqrt{2}}$
32 |  0.13 |  0.51 |  0.86
64 |  0.81 |  4.80 |  3.32
128 |  1.74 |  4.78 |  1.93
256 |  1.96 |  5.09 |  1.35
512 |  1.99 |  4.18 |  1.85
$k_{x}=k_{y}=2$, $r=\frac{1}{6\sqrt{2}}$
32 |  1.31 |  5.21 |  1.88
64 |  1.86 |  4.40 |  1.97
128 |  1.97 |  3.96 |  1.99
256 |  2.00 |  3.92 |  1.99
512 |  2.00 |  3.94 |  2.00
$k_{x}=k_{y}=21$, $r=\frac{1}{6\sqrt{2}}$
32 |  0.02 |  0.30 |  0.43
64 |  3.17 |  4.28 |  0.55
128 |  7.46 |  4.00 |  1.61
256 |  -0.97 |  3.88 |  1.94
512 |  1.66 |  3.93 |  1.99
Table 4. Grid convergence for the three schemes we have chosen. $N$ is the grid dimension. PPW is the points-per-wavelength ratio. $T=\frac{1}{\sqrt{2}}$, $Z=1$.  $r$ | NC | C4 | $AI^{h}$
---|---|---|---
$k_{x}=k_{y}=2$ (PPW 64)
$\frac{1}{6\sqrt{2}}$ | 2.86e-05 | 4.06e-07 | 7.92e-04
$\frac{2}{6\sqrt{2}}$ | 1.18e-04 | 3.38e-07 | 3.86e-03
$\frac{3}{6\sqrt{2}}$ | 2.66e-04 | 2.26e-07 | 3.05e-03
$\frac{4}{6\sqrt{2}}$ | 4.73e-04 | 1.01e-07 | 1.68e-03
$\frac{5}{6\sqrt{2}}$ | 7.38e-04 | 2.49e-07 | 9.23e-05
$\frac{1}{\sqrt{2}}$ | $\infty$ | $\infty$ | $-$
$k_{x}=k_{y}=21$ (PPW $\sim$6)
$\frac{1}{6\sqrt{2}}$ | 1.30e-01 | 5.72e-02 | 2.57e-01
$\frac{2}{6\sqrt{2}}$ | 3.39e-02 | 4.89e-02 | 2.95e-01
$\frac{3}{6\sqrt{2}}$ | 1.39e-01 | 3.49e-02 | 2.79e-01
$\frac{4}{6\sqrt{2}}$ | 2.93e-01 | 1.58e-02 | 2.55e-01
$\frac{5}{6\sqrt{2}}$ | 2.48e-01 | 2.60e-02 | 6.48e-02
$\frac{1}{\sqrt{2}}$ | $\infty$ | $\infty$ | $-$
Table 5.  The error for different CFL numbers $r$.
$T\\!=\\!\frac{4}{\sqrt{2}}$, $Z=1$, $h\\!=\\!\frac{1}{64}$.

###  5.3. Observations

In Table 4, we verify the fourth-order accuracy of our scheme C4 and provide
the rates of grid convergence for the other two (NC, $AI^{h}$) schemes as
well.

In Table 5, we examine the effect of the CFL number on the mean error and
verify the results of our stability analysis (Section 4). For the points-per-
wavelength (PPW) ratio of 64 ($k_{x}=k_{y}=2$), the mean error is smaller for
scheme C4, as expected. Moreover, the NC scheme is more accurate than AI since
it is of a higher spatial order, For the PPW ratio of approximately 6.4, the
AI scheme shows a similar error to that of C4 even though it is of only second
order. This is due to the fact that the AI scheme was trained on coarse grids
with low number of points of points per wavelength.

In Figure 3, we examine the mean error as a function of the wave number. As
expected, as the PPW ratio decreases, the C4 scheme no longer outperforms the
other schemes. This is consistent with the Taylor approximation where the
local truncation error depends on higher-order derivatives that increase for
shorter wavelengths (larger wavenumbers $\sqrt{k_{x}^{2}+k_{y}^{2}}$).

Figure 3.  In both panels $k_{y}=k_{x}$, $T=\frac{4}{\sqrt{2}}$, $Z=1$, and
$h=\frac{1}{64}$. The points-per-wavelength ratio ranges from 64
($\sqrt{k_{x}^{2}+k_{y}^{2}}=2$) to $\sim 2.5$
($\sqrt{k_{x}^{2}+k_{y}^{2}}=50$).

In Figures 4 and 5, we examine how the error behaves as a function of the time
step. As expected, scheme C4 is more accurate when the wave number is low
($\sqrt{k_{x}^{2}+k_{y}^{2}}=2$) and the local truncation errors are
relatively small. However, we see that for $k_{x}=k_{y}=11$ the non-compact
scheme is more accurate for $r=\frac{1}{6\sqrt{2}}$ while C4 is more accurate
for $r=\frac{5}{6\sqrt{2}}$ when the time step increases and as a result the
temporal error is larger.

Figure 4.  $h=\frac{1}{64},r=\frac{1}{6\sqrt{2}}$, $Z=1$.

Figure 5. $h=\frac{1}{64},r=\frac{5}{6\sqrt{2}}$, $Z=1$.

Moreover, we see that for $r=\frac{1}{6\sqrt{2}}$, the scheme NC is more
accurate than AI while the situation is the opposite for
$r=\frac{5}{6\sqrt{2}}$ . This is not a surprise since the dispersion effects
may grow as the time step increases [2].

For higher wave numbers ($\sqrt{k_{x}^{2}+k_{y}^{2}}=29$ and
$\sqrt{k_{x}^{2}+k_{y}^{2}}=51$), the fourth-order accuracy of C4 and the
spatial fourth-order accuracy of NC lose their advantage due to the pollution
effect, and as the wave number increases the AI scheme demonstrates a similar
accuracy even though it is only second-order accurate in both time and space.

##  6\. Conclusions

We have constructed a compact implicit scheme for the 3D Maxwell’s equations.
The scheme is fourth-order accurate in both space and time and can therefore
also be useful for long-time simulations, including scenarios where the PPW
ratio is not very high. The scheme is compact and maintains its accuracy near
the boundaries using equation-based approximations. The temporal grid is
staggered, which allows us to solve 3 scalar uncoupled elliptic equations at
each half-time step. That can be done in parallel with the help of the
conjugate gradient method. The elliptic equations are strictly positive
definite, which enables rapid convergence of the iterative scheme. With fixed
CFL number, approximately three iterations of conjugate gradient are required
for convergence regardless the grid size. This is likely explained by the fact
that the amount of positivity in equations (3.1) increases as the grid length
decreases because the quantity $\frac{24}{h_{\tau}^{2}}$ becomes larger as $h$
decrease and the CFL number remains fixed. Overall, the method is efficient in
both CPU and memory requirements. Although we have tested numerically only 2D
examples, the implementation in 3D is straightforward.

As the method we have presented is finite-difference, the pollution effect
cannot be avoided completely. For higher wavenumbers, lower-order methods
obtained by minimizing a certain loss function over the stencil parameters
(such as data-driven) can outperform higher-order schemes.

## Appendix A

###  A.1. Operator details

Let $\otimes$ denote the Kronecker product of matrices and let $I,J,K$ be
index sets. Let

$D_{1}=\frac{1}{h}\begin{pmatrix}-1&1&0&0&0&\dots&0\\\ 0&-1&1&0&0&\dots&0\\\
0&0&-1&1&0&\dots&0\\\ \vdots&\ddots&\ddots&\ddots&\ddots&\ddots&\vdots\\\
0&0&\dots&0&-1&1&0\\\ 0&0&\dots&0&0&-1&1\end{pmatrix}$ (A.1)

denote the standard central finite-difference matrix with grid size $h/2$
whose dimension equals $(m-1)\times m$, where $m\in\\{|I|,|J|,|K|\\}$. We
define the following operators from $I\times J\times K$ to $I\times J\times
K$.

$\displaystyle D_{x}:=D_{1}\otimes\mathbf{1}_{J}\otimes\mathbf{1}_{K},$
$\displaystyle D_{y}:=\mathbf{1}_{I}\otimes D_{1}\otimes\mathbf{1}_{K},$
$\displaystyle D_{z}:=\mathbf{1}_{I}\otimes\mathbf{1}_{J}\otimes D_{1}.$

Next, we consider a fourth-order accurate approximation for first derivatives:
Assume that $f(x)$ is known at $N$ points $x_{0},x_{1},..,x_{N-1}$. Then, one
estimates $f^{\prime}(x+h/2)$ at the $N-1$ points
$x_{\frac{1}{2}},..,x_{N-3/2}$ to fourth order as:

$\underbrace{\frac{1}{24}\begin{pmatrix}26&-5&4&-1&0&\dots&0\\\
1&22&1&0&0&\dots&0\\\ 0&1&22&1&0&\dots&0\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\ddots&\vdots\\\ 0&0&\dots&1&22&1&0\\\
0&0&\dots&0&1&22&1\\\
0&0&\dots&-1&4&-5&26\end{pmatrix}}_{A}\begin{pmatrix}f^{\prime}(x_{\frac{1}{2}})\\\
f^{\prime}(x_{3/2})\\\ \cdots\\\ \cdots\\\ \cdots\\\ \cdots\\\
f^{\prime}(x_{N-3/2})\end{pmatrix}=\frac{1}{h}\begin{pmatrix}f(x_{1})-f(x_{0})\\\
\cdots\\\ \cdots\\\ \cdots\\\ \cdots\\\ \cdots\\\
f(x_{N-1)}-f(x_{N-2})\end{pmatrix}+O(h^{4}).$ (A.2)

We define the following operators from $I\times J\times K$ to $I\times J\times
K$.

$\displaystyle a_{x}^{-1}:=A^{-1}\otimes\mathbf{1}_{J}\otimes\mathbf{1}_{K},$
$\displaystyle a_{y}^{-1}:=\mathbf{1}_{I}\otimes A^{-1}\otimes\mathbf{1}_{K},$
$\displaystyle a_{z}^{-1}:=\mathbf{1}_{I}\otimes\mathbf{1}_{J}\otimes A^{-1},$

and

$\displaystyle\delta_{x}:=a_{x}^{-1}D_{x},$
$\displaystyle\delta_{y}:=a_{y}^{-1}D_{y},$
$\displaystyle\delta_{z}:=a_{z}^{-1}D_{z}.$

Finally, we define a fourth-order approximation of the curl operator as a
matrix that operates on the vectorized tensor of dimension
$|I|\cdot|J|\cdot|K|\times 1$ :

${\bf curl}_{h}:=\begin{pmatrix}0&-\delta_{z}&\delta_{y}\\\
\delta_{z}&0&-\delta_{x}\\\ -\delta_{y}&\delta_{x}&0\\\ \end{pmatrix}.$

######  Remark A.1.

By [6, section 11], $\|{\bf curl}_{h}\|\leq\frac{2\|A^{-1}\|\sqrt{3}}{h}$ and
by [18], $\|A^{-1}\|\sim\frac{6}{5}$.

We recall the standard second-order finite difference matrices for the second
derivative. We define the operators $D_{xx},D_{yy},D_{zz}$ from $I\times
J\times K$ to $I\times J\times K$ using the square matrix

$D_{2}=\frac{1}{h^{2}}\begin{pmatrix}-2&1&0&0&0&\dots&0\\\
1&-2&1&0&0&\dots&0\\\ 0&1&-2&1&0&\dots&0\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\ddots&\vdots\\\ 0&0&\dots&1&-2&1&0\\\
0&0&\dots&0&1&-2&1\\\ 0&0&\dots&0&0&1&-2\end{pmatrix},$ $\displaystyle
D_{xx}:=D_{2}\otimes\mathbf{1}_{J}\otimes\mathbf{1}_{K},$ $\displaystyle
D_{yy}:=\mathbf{1}_{I}\otimes D_{2}\otimes\mathbf{1}_{K},$ $\displaystyle
D_{zz}:=\mathbf{1}_{I}\otimes\mathbf{1}_{J}\otimes D_{2}.$

###  A.2. Operator estimates

Let $r=\frac{h_{\tau}}{h}$ be the CFL number. Consider the finite-difference
operators

$P_{1}:=-\left(\Delta_{h}+\frac{h^{2}}{6}\varUpsilon_{h}\right)+\frac{24}{h_{\tau}^{2}}\left(1+\frac{2}{r^{2}}\right),$

and

$P_{2}:=\frac{24}{h_{\tau}^{2}}\left(-1-\frac{2}{r^{2}}-\frac{h^{2}}{12}\Delta_{h}\right),$

where

$\Delta_{h}\stackrel{{\scriptstyle\text{def}}}{{=}}D_{xx}+D_{yy}+D_{zz},\quad\varUpsilon_{h}\stackrel{{\scriptstyle\text{def}}}{{=}}D_{zz}D_{xx}+D_{yy}D_{zz}+D_{xx}D_{yy}.$

Let $\sigma(\cdot)$ denote the spectrum of a given operator. The inclusion

$\sigma(h^{2}\Delta_{h})\subset[-12,0]$

implies that

$\frac{h_{\tau}^{2}}{24}\sigma(-P_{2})\subset[1+2/r^{2}-1,1+2/r^{2}].$

The operator
$\frac{h_{\tau}^{2}}{24}\left(-\Delta_{h}-\frac{h^{2}}{6}\varUpsilon_{h}\right)+1+\frac{2}{r^{2}}$
operates on a Fourier ansatz
$\exp(\sqrt{-1}(i\theta_{x}+j\theta_{y}+k\theta_{z})$ by

$\displaystyle
1+\frac{2}{r^{2}}+\frac{r^{2}}{24}\underbrace{\left(6-2(\cos\theta_{x}+\cos\theta_{y}+\cos\theta_{z})-\frac{4T}{6}\right)}_{S},$

where

$T=(\cos(\theta_{z})-1)(\cos(\theta_{x})-1)+(\cos(\theta_{x})-1)(\cos(\theta_{y})-1)+(\cos(\theta_{z})-1)(\cos(\theta_{y})-1),$

and $-4\leq S\leq 12$. Hence,

$\displaystyle\begin{split}\frac{h_{\tau}^{2}}{24}\sigma(P_{1})\subset[1+2/r^{2}-r^{2}/6,1+2/r^{2}+r^{2}/2],\\\
\frac{h_{\tau}^{2}}{24}\sigma(-P_{2})\in[1+\frac{2}{r^{2}}-1,1+\frac{2}{r^{2}}].\end{split}$
(A.3)

As a result, $r<\sqrt{3+\sqrt{21}}$ implies that $-P_{2},P_{1}$ are positive
definite and

$0<\frac{2}{r^{2}}\leq\|\frac{h_{\tau}^{2}}{24}P_{2}\|\leq 1+\frac{2}{r^{2}},$

and

$0<1+\frac{2}{r^{2}}-\frac{r^{2}}{6}\leq\|\frac{h_{\tau}^{2}}{24}P_{1}\|\leq
1+\frac{2}{r^{2}}+\frac{r^{2}}{2}.$

###  A.3. The data-driven scheme

We follow the approach of [9]. The general framework for a data-driven scheme
can be described as follows. Let $\\{(X^{n},Y^{n})\\}_{n}$ be $n$-tuples of
data-points such that for any $n$, $X^{n}\in\mathbb{R}^{m}$ and
$Y^{n}\in\mathbb{R}^{l}$. We define a network $\mathcal{N}_{\vec{a}}$ as a
function from $\mathbb{R}^{m}\to\mathbb{R}^{l}$ which depend on parameters
$\vec{a}$. We define a loss function of the form
$\mathcal{L}_{\vec{a}}:=\sum_{n}|\mathcal{N}_{\vec{a}}(X^{n})-Y^{n}|$ where
$|\cdot|$ is a given norm between two vectors in $\mathbb{R}^{l}$. We strive
to solve the problem

$\underset{\vec{a}}{\mathrm{argmin}}\ \mathcal{L}_{\vec{a}}$

and then obtain the optimal parameters $\vec{a}$ and the corresponding optimal
network $\mathcal{N}_{\vec{a}}$. The minimization process can be done using
several variations of the gradient descent method.

In our framework, we wish to find the optimal parameters $a,b,d$ for
evaluation of the first derivative using the stencil

$K_{2}=\frac{1}{\Delta x}\begin{pmatrix}-b&-a&a&b\\\
-d&-1+3d+2a+6b&1-3d-2a-6b&d\\\ -b&-a&a&b\end{pmatrix}$

in Yee updating rules (see subsection 5.2). The network will be a function
that takes a solution ${\bf E}^{n},{\bf H}^{n+1/2}$ and returns ${\bf
E}^{n+1},{\bf H}^{n+3/2}$ using Yee updating rules. We generate a data set of
analytical solutions using (5.3a) as follows. We fix $h=\frac{1}{16}$, CFL
number $r$ and the corresponding $h_{\tau}$. We also fix the final time $T$.
The number of spatial grid points is then denoted by $N$ and the number of
times steps is denoted by $N_{\tau}$. For any $(k_{x},k_{y})$ in
$\\{12,13,14,15\\}^{2}$ we generate analytical solutions $E_{z},H_{x},H_{y}$
defined by (5.3a). For any $0\leq n\leq N_{\tau}$ we let
${E_{z}}^{n}_{k_{x},k_{y},\mathrm{true}},{H_{x}}^{n}_{k_{x},k_{y}\mathrm{true}},{H_{y}}^{n}_{k_{x},k_{y},\mathrm{true}}$
be the corresponding values of the analytical solutions evaluated at the time
$t=nh_{\tau}$. We let

$X^{n,k_{x},k_{y}}:=({E_{z}}^{n}_{k_{x},k_{y},\mathrm{true}},{H_{x}}^{n}_{k_{x},k_{y}\mathrm{true}},{H_{y}}^{n}_{k_{x},k_{y},\mathrm{true}}).$

Next, we build the network, $\mathcal{N}_{a,b,d}$, which takes a numerical
solution at time step $n$, and outputs the solution at time $n+1$ evaluated
using Yee updating rules with the stencil $K_{2}$.

A single data point in our data set is then defined by
$(\mathrm{input}^{n,k_{x}.k_{y}},\mathrm{output}^{n,k_{x},k_{y}})$ where

$\mathrm{input}^{n,k_{x},k_{y}}=X^{n,k_{x},k_{y}}$

and

$\mathrm{output}^{n,k_{x},k_{y}}=(X^{n+1,k_{x},k_{y}},X^{n+2,k_{x},k_{y}},X^{n+3,k_{x},k_{y}}).$

We omit the superscripts $k_{x},k_{y}$ from now on.

The loss function is then defined over all data points as follows:

$\mathcal{L}_{a,b,d}=\sum_{n}|X^{n+1}-\mathcal{N}(X^{n})|+|X^{n+2}-\mathcal{N}(\mathcal{N}(X^{n}))|+|X^{n+3}-\mathcal{N}(\mathcal{N}(\mathcal{N}(X^{n})))|$
(A.4)

where $|\cdot|$ in (A.4) is taken to be the mean absolute error between two
vectors. The loss function is minimized using Adam optimizer and Keras [4]
over the data points with the usual splitting routine of train, test, and
validation sets.

## References

  * [1] Bayliss, A., Goldstein, C.I. and Turkel, E., 1985. On accuracy conditions for the numerical computation of waves. Journal of Computational Physics, 59(3), pp. 396-404.
  * [2] Blinne, A., Schinkel, D., Kuschel, S., Elkina, N., Rykovanov, S.G. and Zepf, M., 2018. A systematic approach to numerical dispersion in Maxwell solvers. Computer Physics Communications, 224, pp. 273-281.
  * [3] Britt, S., Tsynkov, S. and Turkel, E., 2011. Numerical simulation of time-harmonic waves in inhomogeneous media using compact high order schemes. Communications in Computational Physics, 9(3), pp. 520-541.
  * [4] Chollet, F. & others, 2015. Keras. Available at: https://github.com/fchollet/keras.
  * [5] Deraemaeker, A., Babuška, I. and Bouillard, P., 1999. Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions. International journal for numerical methods in engineering, 46(4), pp. 471-499.
  * [6] Gustafsson, B., Kreiss, H.O. and Oliger, J., 1995. Time dependent problems and difference methods (Vol. 24). John Wiley & Sons.
  * [7] Gottlieb, D. and Yang, B., 1996, March. Comparisons of staggered and non-staggered schemes for Maxwell’s equations. In 12th Annual Review of Progress in Applied Computational Electromagnetics (Vol. 2, pp. 1122-1131).
  * [8] Kreiss, H.O. and Oliger, J., 1972. Comparison of accurate methods for the integration of hyperbolic equations. Tellus, 24(3), pp. 199-215.
  * [9] Ovadia, O., Kahana, A. and Turkel, E., 2022. A Convolutional Dispersion Relation Preserving Scheme for the Acoustic Wave Equation. arXiv preprint arXiv:2205.10825.
  * [10] L. D. Landau and E. M. Lifshitz, 1975. Course of Theoretical Physics, Vol. 2, The Classical Theory of Fields, Fourth ed., Pergamon Press, Oxford.
  * [11] L. D. Landau and E. M. Lifshitz, 1984. Course of Theoretical Physics. Vol. 8, Electrodynamics of Continuous Media, Pergamon International Library of Science, Technology, Engineering and Social Studies, Pergamon Press, Oxford.
  * [12] Leis, R., 2013. Initial boundary value problems in mathematical physics. Courier Corporation.
  * [13] Sakkaplangkul, P. and Bokil, V., 2021. Convergence analysis of Yee-FDTD schemes for 3D Maxwell’s equations in linear dispersive media. International journal of numerical analysis and modeling, 18(4).
  * [14] Singer, I. and Turkel, E., 1998. High-order finite difference methods for the Helmholtz equation. Computer methods in applied mechanics and engineering, 163(1-4), pp. 343-358.
  * [15] Stern, A., Tong, Y., Desbrun, M. and Marsden, J.E., 2015. Geometric computational electrodynamics with variational integrators and discrete differential forms. Geometry, Mechanics, and Dynamics: The Legacy of Jerry Marsden, pp. 437-475.
  * [16] Taflove, A., Hagness, S.C. and Piket-May, M., 2005. Computational electromagnetics: the finite-difference time-domain method, Artech house, Boston 1998.
  * [17] Yee, K., 1966. Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media. IEEE Transactions on antennas and propagation, 14(3), pp. 302-307.
  * [18] Yefet, A. and Turkel, E., 2000. Fourth order compact implicit method for the Maxwell equations with discontinuous coefficients. Applied Numerical Mathematics, 33(1-4), pp. 125-134.
  * [19] Yefet, A. and Petropoulos, P.G., 2001. A staggered fourth-order accurate explicit finite difference scheme for the time-domain Maxwell’s equations. Journal of Computational Physics, 168(2), pp. 286-315.


# GymD2D: A Device-to-Device Underlay Cellular Offload Evaluation Platform

David Cotton  0000-0002-8817-3736 School of Electrical and Data Engineering
University of Technology Sydney
Sydney, Australia
<EMAIL_ADDRESS>Zenon Chaczko  School of Electrical and
Data Engineering
University of Technology Sydney
Sydney, Australia
<EMAIL_ADDRESS>
###### Abstract

Cellular offloading in device-to-device communication is a challenging
optimisation problem in which the improved allocation of radio resources can
increase spectral efficiency, energy efficiency, throughout and reduce
latency. The academic community have explored different optimisation methods
on these problems and initial results are encouraging. However, there exists
significant friction in the lack of a simple, configurable, open-source
framework for cellular offload research. Prior research utilises a variety of
network simulators and system models, making it difficult to compare results.
In this paper we present GymD2D, a framework for experimentation with physical
layer resource allocation problems in device-to-device communication. GymD2D
allows users to simulate a variety of cellular offload scenarios and to extend
its behaviour to meet their research needs. GymD2D provides researchers an
evaluation platform to compare, share and build upon previous research. We
evaluated GymD2D with state-of-the-art deep reinforcement learning and
demonstrate these algorithms provide significant efficiency improvements.

###### Index Terms:

device-to-device (D2D) communication, cellular offload, resource allocation,
radio resource management, network simulator, deep reinforcement learning,
OpenAI Gym

††publicationid: pubid:  ©2021 IEEE. Personal use of this material is
permitted. Permission from IEEE must be obtained for all other uses, in any
current or future media, including reprinting/republishing this material for
advertising or promotional purposes, creating new collective works, for resale
or redistribution to servers or lists, or reuse of any copyrighted component
of this work in other works.

##  I Introduction

Wireless data use is increasing rapidly presenting a challenge for cellular
service providers. The explosion of video and smart device traffic has
intensified demand for limited cellular resources [1]. A multi-faceted
approach to developing next generation networks include: unlocking millimeter-
wave frequencies, increasing cell density, multiple-input multiple-output,
improved information encoding, and the focus of this research—smarter
coordination [2].

Device-to-device (D2D) communication is a broad set of protocols for ad-hoc,
peer-to-peer wireless communication for cellular or Internet connected
wireless devices. In contrast to normal cellular operation, user equipment
(UE) utilising D2D mode communicate directly with one another instead of
through base stations and connected networks [3]. D2D has been proposed for a
variety of use cases such as: public safety communications, communications
relay, localised services, Internet of things and cellular offload [4].

Figure 1: Example of D2D cellular offloading with underlay networks. This
example demonstrates the cellular offload optimisation problem. Depicted is a
BS, CUE and 2 DUE pairs. The BS has 2 RBs available to allocate. The BS has
assigned the CUE RB1 for downlink and has the option of either assigning pair
A RB1 and pair B RB2 or vice versa. If the BS assigns RB1 to pair A it is
likely that due to their proximity to the CUE, significant interference would
occur. As pair B is situated relatively distant to the CUE, if pair B
transmitted with a lower power level, their interference on the CUE would be
negligible. Therefore the optimal solution is to assign RB1 to pair B and RB2
to pair A.

In this research we are interested in cellular offload, a mechanism to divert
cellular traffic through alternative channels to improve network efficiency.
This can be achieved by communicating out of band, such as over WiFi, or
inband via overlay or underlay networking. Underlay cellular offload is a form
of opportunistic spectrum access in which D2D UE (DUE) act as secondary users
sharing radio resources with primary cellular UE (CUE) [5]. In underlay
networking, DUE are responsible for managing their interference with CUE so as
to avoid degrading primary network performance.

Radio resource management (RRM) is the system wide management of radio
resources, such as transmit power and resource blocks (RB), across a wireless
network to manage interference and utilise resources as efficiently as
possible. In cellular systems, which are limited by co-channel interference,
improved resource allocation can increase spectral efficiency, energy
efficiency, throughput and reduce latency. We show a simplified example of the
resource allocation problem in Figure 1.

A challenge facing researchers developing D2D algorithms is the need for
improved tooling, reliable measurements and established benchmarks. An open-
source evaluation platform allows researchers to compare results, share
algorithms and build upon previous research.

In this paper we present GymD2D, a network simulator and evaluation platform
for RRM in D2D underlay cellular offload. GymD2D provides convenient
abstractions to aid researchers in quick prototyping of resource allocation
algorithms. The toolkit allows users to programmatically configure the
environment to simulate a wide variety of scenarios. GymD2D has been designed
with extensibility as a core design principle, allowing users to override and
extend its behaviour to meet their research needs. It has been developed in
the Python programming language to allow users to leverage its extensive
ecosystem of scientific computing packages. The open-source nature of GymD2D
centralises development effort and avoids the redundant work of individual
researchers creating their own simulators. This puts more eyes on bug fixing,
provides a more stable platform and increases confidence in reported empirical
results. GymD2D reduces entry barriers for junior researchers, helps
researchers from other disciplines to cross-pollinate ideas easier and more
generally increases participation. Our software package is provided to the
community under a MIT licence at https://github.com/davidcotton/gym-d2d.

##  II Background

###  II-A Device-to-device communication

In this section we provide an overview of the D2D RRM literature to situate
the reader as to the requirements of the platform. Firstly, we highlight the
most common optimisation problems. Secondly, we analyse key differences across
simulators, paying special attention to simplifying assumptions frequently
observed. Thirdly, we survey the optimisation algorithms used for resource
allocation. Finally, we outline limitations of existing research, providing
direction for the simulation requirements of future work.

####  II-A1 Optimisation Problems

RRM is an optimisation problem where the objective is to utilise radio
resources as efficiently as possible. D2D RRM has been proposed for improving
spectral efficiency, energy efficiency and quality of service. In these uses
cases, the objective can be to optimise data rates, throughput, capacity,
signal to interference noise ratio (SINR), power consumption, energy
efficiency or latency [4, 6]. D2D systems can be centrally managed by the
network operator, manage interference autonomously or use a hybrid control
mode which aims to combine the benefits of both. The choice of control mode
limits the applicability of certain algorithms which may only be feasible in
centrally managed paradigms.

####  II-A2 Simulation models

D2D cellular offload typically investigates networks using orthogonal
frequency division multiple access (OFDMA), communicating on licensed bands
using underlay networking. The most common scenario is a single macro base
stations (MBS) surrounded by many randomly positioned CUEs and DUEs. It is
generally assumed that cellular systems are under full load and each RB is
allocated to a CUE. In the literature, simulations vary in scope from 2–30
RBs, 2–30 CUEs and 2–60 DUEs, while MBS operate with a cell radius of 20–500m.
It is frequently assumed that DUE are already paired and operate in a range
between 10–30m apart. Typically, omni-directional antenna and isotropic
propagation are utilised. Path loss is commonly modeled using log-distance
models, with or without shadowing.

####  II-A3 Optimisation Algorithms

A wide variety of optimisation methods have been investigated on a range of
D2D RRM problems. Initially, D2D radio resources were proposed to be managed
using existing cellular uplink power control mechanisms [7]. Consequently, it
was identified that resources could be more efficiently allocated with the use
of mathematical optimisation [8]. However, due to the computational complexity
of these methods and the millisecond timescales involved, it may not feasible
to solve to optimality. This can be addressed with the use of greedy heuristic
algorithms which reduce computational complexity at the cost of global
optimality [9]. Alternatively, resource allocation can be optimised graph
theoretically [10], game theoretically [11], with evolutionary algorithms [12]
or reinforcement learning (RL) [13]. More recently, deep reinforcement
learning (DRL), a subfield of RL which uses deep neural networks to represent
policy and/or value functions has demonstrated promising results [14, 15]. DRL
is well suited for many D2D RRM problems as neural networks provide rich
approximations, scale well and generalise to unseen data.

####  II-A4 Research limitations

A common limitation observed in D2D cellular offload research is BSs not
enforcing uplink power control for CUE, who transmit at maximum power, a very
energy inefficient approach. Another research challenge is accounting for
large SINR increases on the primary network, such as could significantly
impact primary network throughput or drive up CUE transmit power levels.
Resource allocations algorithms need to demonstrate their effectiveness in
larger search spaces that more closely reflect real world demands. Iterative
learning algorithms need to be capable of generalising to out of training
distribution data and be robust under diverse propagation conditions. Lastly,
in our opinion, one of the greatest limitations of existing research is the
lack of established benchmarks and comparison with other algorithms.

###  II-B OpenAI Gym

OpenAI Gym is an open-source software toolkit for reinforcement learning (RL)
[16]. Gym provides an abstraction layer that enables a variety of tasks, known
as environments in RL parlance, to be wrapped to present a consistent
interface. The abstraction provided by Gym allows the easy interchange of
algorithms and environments. This makes it is easy to test how a single
algorithm generalises across a diverse set of environments or to benchmark
different algorithms on a given environment. The simplicity and flexibility
Gym offers has proved very popular and has lead to it becoming the de facto
environment format in RL. While Gym was designed for RL research, the
application programming interface (API) it provides makes it easy to apply
many other algorithms types.

###  II-C Network simulation

One of the most widely used network simulators in education and research is
ns-3. Ns-3 is an open-source, modular, discrete-event simulator for wired and
wireless networks. It provides the full TCP/IP stack and wireless propagation
modelling. Another popular alternative with comparable features is OMNeT++,
while there exists similar commercial tools such as NetSim and MATLAB. Ns-3
has been incorporated into an OpenAI Gym environment under the ns3-gym project
[17].

##  III Design Principles

The design of GymD2D has been inspired by the authors experience developing
and comparing reinforcement learning algorithms. In our experience the
following design principles stimulate experimentation and the sharing of
ideas.

  * •

Simple: Easy to get started with, the framework should allow researchers to be
productive quickly.

  * •

Configurable: The framework should be easily configured to meet the broad
range of D2D cellular offload use cases. Configurability allows researchers to
programmatically test algorithm generalisation and scalability.

  * •

Extensible: The framework should allow users to extend the system’s behaviour
to meet their needs. The nature of research dictates a stream of new ideas we
can’t anticipate but we can provide researchers the flexibility to adapt.

  * •

Scalable: The framework should be performant and easily parallelisable.
Developing new algorithms requires significant experimentation and reducing
the time spent waiting for results is important for productivity. Some
algorithms, such as policy gradient DRL, require parallel environments to
function. Real world solutions are often a combination of both algorithmic and
architectural components.

  * •

Reproducible: Experiments should be easily repeatable. To build confidence in
our deductions, it is important that we can reperform experiments to ensure
the observed outcomes were not statistical anomalies. Reproducibility allows
researchers to share their contributions with community more easily.

##  IV System Design

###  IV-A System model

GymD2D is designed to simulate physical layer resource allocation problems in
D2D underlay cellular offload. The framework abstracts away data link and
above layers, D2D session establishment and management concerns. The system
model contains a single MBS $b$, a set of $M$ CUEs
$\mathcal{M}=\\{1,\dots,M\\}$ and a set of $N$ DUE pairs
$\mathcal{N}=\\{1,\dots,N\\}$, that reside within the coverage area of the
cell. We denote the $m^{th}$ CUE $C_{m}$, the $n^{th}$ DUE pair $D_{n}$ and
the transmitter and receiver of pair $D_{n}$ by $D_{n}^{t}$ and $D_{n}^{r}$
respectively.

The system employs OFDMA, with a set of $K$ RBs $k\in\mathcal{K}$ are
available for allocation. An assumption is made that all devices are equipped
with omni-directional antenna and transmit isotropically. Accordingly, the
network resides within a circular cell of radius $R$, with the MBS located in
the centre at position $(0,0)$. The simulation environment contains no
obstructions or outside interference. D2D communicate one-to-one and D2D relay
is not supported.

We denote the effective isotropic radiated power (EIRP) $P$ of BSs, CUEs and
DUEs as $P^{b}$, $P^{c}$, $P^{d}$ respectively. The EIRP of a BS is
calculated,

$P^{b}=P_{tx}-10log_{10}s+g_{ant}-l_{ix}-l_{cb}+g_{amp}$ (1)

and the EIRP of CUE and DUE,

$P^{c}=P^{d}=P_{tx}-10log_{10}s+g_{ant}-l_{ix}-l_{bd}$ (2)

where $P_{tx}$ is the transmission power level, $s$ is the number of
subcarriers, $g_{ant}$ is the transmitting antenna gain, $l_{ix}$ is the
interference margin loss to approximate noise from surrounding cells, $l_{bd}$
is body loss to approximate attenuation caused by the user, $l_{cb}$ is cable
loss, and $g_{amp}$ is amplifier gain.

We denote the received signal level $R$ from transmitter $i$ at receiver $j$
of BS, CUEs and DUEs as $R^{b}_{i,j}$, $R^{c}_{i,j}$, $R^{d}_{i,j}$. The
received signal level of BS as,

$R^{b}_{i,j}=P_{i}-PL_{i,j}+g_{ant}-l_{cb}+g_{amp}$ (3)

and the received signal level of CUE or DUE,

$R^{c}_{i,j}=R^{d}_{i,j}=P_{i}-PL_{i,j}+g_{ant}-l_{bd}$ (4)

where $P_{i}$ is the EIRP from transmitter $i$ and $PL_{i,j}$ is the path loss
of the chosen path loss model between $i$ and $j$.

We assume D2D transmissions are synchronised to cellular transmissions and
occupy the same $K$ orthogonal resources. During both uplink and downlink, co-
channel interference is calculated for each receiver sharing RB $k$. GymD2D
considers co-channel interference between:

  * •

D2D to cellular, interference from secondary DUE on the primary cellular
network,

  * •

cellular to D2D, interference from CUE or BS to DUE, and

  * •

D2D to D2D, the interference between DUE pairs sharing a RB.

Accordingly, we model the instantaneous SINR $\xi$ of receiver $j$ from
transmitter $i$ on RB $k$,

$\xi_{i,j,k}=\frac{R_{i,j}}{\sum_{n\in\mathcal{T}_{k},n\neq
i}R_{n,j}+\sigma^{2}}$ (5)

where $\mathcal{T}_{k}$ is the set of transmitters allocated to RB $k$ and
$\sigma^{2}$ is additive white Gaussian noise (AWGN).

The capacity of channel $C_{i,j}$ can be calculated using the SINR
$\xi_{i,j}$,

$C_{i,j}[Mbps]=Blog_{2}(1+\xi_{i,j})$ (6)

where $B$ is the channel bandwidth in MHz.

###  IV-B Path loss models

GymD2D contains several of the most common path loss models and makes it easy
for users to implement their own custom models. By default, GymD2D uses the
simplest model, free space path loss (FSPL),

$FSPL(f,d)[dB]=10nlog_{10}\Big{(}\frac{4\pi fd}{c}\bigg{)}$ (7)

where $n=2$ is the path loss exponent (PLE) in free space, $f$ is the carrier
frequency in Hz, $d$ is the distance between the transmitter and receiver and
$c$ is the speed of light in m/s.

To simulate obstructed propagation environments it can be useful to model
fading effects as random processes. One such model is the log-distance with
shadowing path loss model, which is included in GymD2D. The log-distance path
loss model extends FSPL to mimic random shadowing effects, such as caused by
buildings, with a log-normal distribution,

$PL^{LD}(f,d)[dB]=FSPL(f,d_{0})+10nlog_{10}\frac{d}{d_{0}}+\chi_{\sigma}$ (8)

where $d_{0}$ is an arbitrary close-in reference distance, typically 1–100m
and $\chi_{\sigma}$ is a zero-mean Gaussian with standard deviation $\sigma$
in dB. Empirical measurements have shown values of $n=2.7\text{ to }3.5$ to be
suitable to model urban environments [18].

###  IV-C Architecture

GymD2D consists of two main components, a network simulator and a Gym
environment. The network simulator models physical layer cellular networking.
The Gym environment provides an abstraction layer to allow researchers to
experiment with different simulation parameters and algorithms
programmatically. Users supply RRM algorithms to manage the wireless devices
under simulation. GymD2D outputs data on the state of the simulation to the
user; to allow the effectiveness of RRM algorithms to be studied, such as
through visualisation. A high level overview of the architecture of GymD2D is
depicted in Figure 2.

Figure 2: Proposed GymD2D architecture. GymD2D consists of a network
simulator, wrapped by an OpenAI Gym environment. The user creates their own
RRM algorithms to control wireless devices.

###  IV-D Network simulator

The network simulator models a single cellular cell which is populated with a
collection of randomly placed CUEs and DUE pairs. It is a configurable
component which can be customised to emulate a range of cellular offload
scenarios. This includes the number and configuration of BSs, CUEs and DUEs
and environmental parameters such as the available RBs, cell size and path
loss model.

The main components of the network simulator are: a collection of wireless
devices (BSs, CUEs, DUEs), a path loss model and a traffic model as shown in
the class diagram in Figure 3.

Figure 3: Network simulator architecture. The main components of the network
simulator are a collection of CUEs, DUEs and BS, the path loss model and the
traffic model.

Each simulation, the actions of BS and UEs within the cell can be generated
internally by the traffic model or externally from a user defined RRM
algorithm. A typical use case would be to use the internal traffic model to
control BS and CUEs and the user RRM algorithm the DUEs.

GymD2D uses a discrete-event simulation model. This method is congruent with
the Gym API in which the incoming actions are the events and the Gym step()
method calls equate to the system update intervals and model a single LTE or
NR frame. At each step, each device may transmit, receive, or take no action.
An action is tuple consisting of a transmitter, receiver, communication mode,
RB and transmission power. The simulator consolidates the actions from both
the traffic model and the RRM algorithm, then calculates the resulting
propagation and interference. After calculating propagation, metrics on the
state of the network, such as SINR and throughput, are output to the Gym
environment.

###  IV-E Gym environment

The Gym environment has been designed to be configuration driven, to
facilitate the programmatic scheduling and reproducibility of experiments.
When instantiating a new Gym environment, configuration can be provided to
specify the BSs, CUEs and DUEs that inhabit the simulation and the
environmental conditions.

The Gym environment provides RRM algorithms with an observation and action
space. These provide a mapping to configure for the expected format of inputs
and outputs. For example, when using a DRL algorithm, this would allow DRL to
configure its neural networks for the shape of incoming observations and
output actions of the correct dimension.

At each step, the Gym environment receives actions from the RRM algorithm and
converts them to a format suitable for the network simulator. Once a
simulation step is complete, the environment uses the state of the simulator
to create the observations and rewards RRM algorithms consume to make their
decisions.

##  V Evaluation

###  V-A Methodology

Figure 4: Evaluation results. To evaluate GymD2D we compared the performance
of several state-of-the-art DRL algorithms in their efficiency allocating
radio resources as D2D demand increased. Solid lines indicate the mean
algorithm performance across ten trials with the shaded area the 95%
confidence interval. The dashed red line indicates the baseline total system
capacity without D2D communication. (a) The total system capacity of the DRL
and random agent. (b) The total system capacity of just the DRL agents. (c)
The total DUE capacity of the DRL agents. (d) The mean transmit power of all
agents.

We evaluated GymD2D with several leading DRL algorithms to determine their
efficiency allocating radio resources as D2D demand increased. The objective
was to maximise the total system capacity, that is the sum data rate of all
CUE and DUE, calculated for each transmitter/receiver pair $i,j$ by

$C_{i,j}[Mbps]=\begin{cases}Blog_{2}(1+\xi_{i,j})&\xi_{i,j}\geq\rho_{j}\\\
0&\xi_{i,j}<\rho_{j}\\\ \end{cases}$ (9)

where $B=0.18$ is the RB bandwidth in MHz and $\rho_{b}=-123.4$ and
$\rho_{d}=-107.5$ is the receiver sensitivity of a BS and DUE respectively in
dBm. Our evaluation simulated a single cell under full load. The scenario
contained 25 RBs and CUEs, with each CUE allocated an individual RB. We
employed a centrally managed control mode in which DUE communicated in the
uplink frame, with the resource allocation managed by the network operator.
Each RRM algorithm was evaluated with 10, 20, 30, 40 and 50 communicating D2D
pairs. Algorithms were compared by training to convergence, then evaluating
for 100 episodes. For each algorithm–D2D link density comparison, we conducted
ten trials, retraining from scratch and evaluating, to account for variations
in performance. Each episode lasted for ten steps or equivalently ten LTE/NR
frames to simulate short bursts of traffic on a busy network. In each episode
all CUE and DUE remained geographically fixed, but at the end of each episode,
all CUE and DUE were randomly repositioned within the cell to simulate new
devices accessing the network. Wireless propagation was modelled using the
Log-Distance Shadowing model (8) with PLE $n=2.0$ and $\chi_{\sigma}=2.7$. The
simulation parameters are detailed in Table I.

TABLE I: Simulation parameters Parameter | Value
---|---
Cell radius |  $500$ m
Maximum D2D pair distance |  $30$ m
Carrier frequency |  $2.1$ GHz
RB bandwidth |  $180$ kHz
Number of RBs | $25$
Number of CUEs | $25$
Number of DUE pairs | $10,20,30,40,50$
CUE transmit power |  $23$ dBm
DUE min, max transmit power |  $0,20$ dBm
Path loss model | Log-Distance Shadowing
Path loss exponent | $2.0$
Shadowing SD $\chi_{\sigma}$ | 2.7

We evaluated three DRL algorithms, Rainbow DQN [19], Discrete Soft Actor-
Critic (SAC) [20], and Advantage Actor-Critic (A2C) [21], and a random agent
baseline. All DRL algorithms used a fully connected neural network with two
hidden layers trained using the Adam optimiser. Each hidden layer contained
128 units and used ReLU activation between layers. Policies used a reward
discounting factor of $\gamma=0.9$. Our Rainbow DQN (Table II) implementation
used distributional, dueling, double-Q and noisy networks with a prioritised
replay buffer and single step returns. The discrete action space variant of
the SAC (Table III) was used. A2C (Table IV) is the synchronous version of A3C
and used Generalised Advantage Estimator (GAE) [22] with $\lambda=1.0$.

TABLE II: Rainbow DQN hyperparameters Parameter | Value
---|---
Discounting factor $\gamma$ | $0.9$
Learning rate $\alpha$ | $5\cdot 10^{-4}$
Batch size | $32$
Online network update period |  $4$ steps
Learning start |  $1,000$ steps
Target network sync period |  $500$ steps
Distributional atoms | 51
Distributional bounds $v_{min}$, $v_{max}$ | [-1,10]
Replay buffer capacity | $50,000$
Replay buffer prioritisation exponent $\omega$ | $0.6$
Replay buffer importance sampling $\beta$ | $0.6\rightarrow 0.4$
Importance sampling annealing | 20,000 steps
TABLE III: SAC hyperparameters Parameter | Value
---|---
Discounting factor $\gamma$ | $0.9$
Learning rate $\alpha$ | $3\cdot 10^{-4}$
Batch size | $256$
Learning start |  $1,500$ steps
Target smoothing coefficient | 0.005
Replay buffer capacity | $50,000$
TABLE IV: A2C hyperparameters Parameter | Value
---|---
Discounting factor $\gamma$ | $0.9$
Learning rate $\alpha$ | $1\cdot 10^{-4}$
Rollout length | $10$
Entropy coefficient $\beta$ | $0.01$
GAE $\lambda$ | $1.0$

###  V-B Results

The results of our evaluation can be seen in Figure 4. Baseline total system
capacity measures the efficiency of the system without D2D communication. For
our scenario this was 94.75 Mbps. We found that all the DRL algorithms
achieved a similar level of performance, increasing system capacity over the
baseline by more than 11%. Our results show that the system capacity continued
to increase sublinearly as number of active D2D links demand grew. Conversely,
the performance of the random agent shows that without careful resource
allocation, the system capacity drops sharply.

We found that despite allowing DUE to communicate up to half the power of CUE
(20 vs. 23 dBm), they typically converged into operating ranges between 7 and
15 dBm. This resulted in a negligible decrease in the total CUE capacity, 1–2
Mbps or $\approx$1.84% below the baseline system capacity. This decrease was
approximately constant across D2D density.

###  V-C Discussion

Inspecting the actions the DRL algorithms selected, we found that they
converged to allocating all DUE onto one or two RBs. This is surprising as we
had anticipated the DUE to be evenly distributed amongst all available RBs.
Investigating further, we observed that over the course of a training run, the
DQN converging from an even RB distribution to the focused allocation
strategy. This behaviour developed in the later stages of training and only
contributed modest increases to the system capacity.

As expected, the optimal strategy for resource allocation was to assign DUE to
share RBs with the most geographically distant CUE. When combined with the
focused RB allocation described above, this typically resulted in the RRM
algorithm choosing to allocate DUE to share with the one or two most isolated
CUE.

Despite the random UE positioning, the DRL agents were able to learn policies
that generalised much better than we anticipated when using fully connected
neural networks. We were also surprised how quickly agents adapted during an
episode, improving their performance over the course of the ten-step episode.

##  VI Conclusion

In this research we have presented GymD2D, a network simulator and evaluation
platform for RRM in D2D underlay cellular offload. GymD2D makes it easy for
researchers to build, benchmark and share RRM algorithms and results. Our
toolkit is designed to quickly prototype physical layer resource allocation
algorithms, without the complexity of higher layer protocols. GymD2D is
configurable and extensible, allowing it to be employed to simulate a range of
D2D research needs.

We have evaluated GymD2D with several leading DRL algorithms and demonstrated
the performance gains of intelligent RRM, increasing system capacity by more
than 11%. There was no clear winner amongst the DRL algorithms which performed
similarly. The results also demonstrated that D2D cellular offload can
significantly minimise its impact on primary networks.

In the future we plan to increase the simulation complexity in GymD2D, adding
more realistic modelling. Other interesting research challenge include
investigating the impacts of CUE power control on cellular offload and
supporting D2D relay. We continue to use GymD2D in ongoing research,
developing methods for scaling up DRL based D2D RRM.

## References

  * [1] T. S. Rappaport, W. Roh, and K. Cheun, “Mobile’s millimeter-wave makeover,” _IEEE Spectrum_ , vol. 51, no. 9, pp. 34–58, 2014.
  * [2] A. Gupta and R. K. Jha, “A survey of 5g network: Architecture and emerging technologies,” _IEEE access_ , vol. 3, pp. 1206–1232, 2015.
  * [3] B. Kaufman and B. Aazhang, “Cellular networks with an overlaid device to device network,” in _2008 42nd Asilomar conference on signals, systems and computers_. IEEE, 2008, Conference Proceedings, pp. 1537–1541.
  * [4] A. Asadi, Q. Wang, and V. Mancuso, “A survey on device-to-device communication in cellular networks,” _IEEE Communications Surveys & Tutorials_, vol. 16, no. 4, pp. 1801–1819, 2014.
  * [5] P. Janis, V. Koivunen, C. Ribeiro, J. Korhonen, K. Doppler, and K. Hugl, “Interference-aware resource allocation for device-to-device radio underlaying cellular networks,” in _VTC Spring 2009-IEEE 69th Vehicular Technology Conference_. IEEE, 2009, Conference Proceedings, pp. 1–5.
  * [6] C. Chakraborty and J. J. Rodrigues, “A comprehensive review on device-to-device communication paradigm: Trends, challenges and applications,” _Wireless Personal Communications_ , pp. 1–23, 2020.
  * [7] P. Janis, C.-H. Yu, C. Ribeiro, C. Wijting, K. Hugl, O. Tirkkonen, and V. Koivunen, “Device-to-device communication underlaying cellular communications systems,” _Int’l J. of Communications, Network and System Sciences_ , vol. 2009, 2009.
  * [8] C.-H. Yu, O. Tirkkonen, K. Doppler, and C. Ribeiro, “Power optimization of device-to-device communication underlaying cellular communication,” in _2009 IEEE international conference on communications_. IEEE, 2009, Conference Proceedings, pp. 1–5.
  * [9] M. Zulhasnine, C. Huang, and A. Srinivasan, “Efficient resource allocation for device-to-device communication underlaying lte network,” in _2010 IEEE 6th International conference on wireless and mobile computing, networking and communications_. IEEE, 2010, Conference Proceedings, pp. 368–375.
  * [10] R. Zhang, X. Cheng, L. Yang, and B. Jiao, “Interference-aware graph based resource sharing for device-to-device communications underlaying cellular networks,” in _2013 IEEE wireless communications and networking conference (WCNC)_. IEEE, 2013, Conference Proceedings, pp. 140–145.
  * [11] C. Xu, L. Song, Z. Han, Q. Zhao, X. Wang, and B. Jiao, “Interference-aware resource allocation for device-to-device communications as an underlay using sequential second price auction,” in _2012 IEEE international conference on communications (ICC)_. IEEE, 2012, Conference Proceedings, pp. 445–449.
  * [12] L. Su, Y. Ji, P. Wang, and F. Liu, “Resource allocation using particle swarm optimization for d2d communication underlay of cellular networks,” in _2013 IEEE wireless communications and networking conference (WCNC)_. IEEE, 2013, Conference Proceedings, pp. 129–133.
  * [13] Y. Luo, Z. Shi, X. Zhou, Q. Liu, and Q. Yi, “Dynamic resource allocations based on q-learning for d2d communication in cellular networks,” in _2014 11th International Computer Conference on Wavelet Actiev Media Technology and Information Processing (ICCWAMTIP)_. IEEE, 2014, Conference Proceedings, pp. 385–388.
  * [14] Z. Li, C. Guo, and Y. Xuan, “A multi-agent deep reinforcement learning based spectrum allocation framework for d2d communications,” _arXiv preprint arXiv:1904.06615_ , 2019.
  * [15] J. Tan, L. Zhang, and Y.-C. Liang, “Deep reinforcement learning for channel selection and power control in d2d networks,” in _2019 IEEE Global Communications Conference (GLOBECOM)_. IEEE, 2019, Conference Proceedings, pp. 1–6.
  * [16] G. Brockman, V. Cheung, L. Pettersson, J. Schneider, J. Schulman, J. Tang, and W. Zaremba, “Openai gym,” _arXiv preprint arXiv:1606.01540_ , 2016.
  * [17] P. Gawłowicz and A. Zubow, “ns3-gym: Extending openai gym for networking research,” _arXiv preprint arXiv:1810.03943_ , 2018.
  * [18] T. S. Rappaport, _Wireless communications: principles and practice_. prentice hall PTR New Jersey, 1996, vol. 2.
  * [19] M. Hessel, J. Modayil, H. Van Hasselt, T. Schaul, G. Ostrovski, W. Dabney, D. Horgan, B. Piot, M. Azar, and D. Silver, “Rainbow: Combining improvements in deep reinforcement learning,” _arXiv preprint arXiv:1710.02298_ , 2017\.
  * [20] P. Christodoulou, “Soft actor-critic for discrete action settings,” _arXiv preprint arXiv:1910.07207_ , 2019.
  * [21] V. Mnih, A. P. Badia, M. Mirza, A. Graves, T. Lillicrap, T. Harley, D. Silver, and K. Kavukcuoglu, “Asynchronous methods for deep reinforcement learning,” in _International Conference on Machine Learning_ , 2016, Conference Proceedings, pp. 1928–1937.
  * [22] J. Schulman, P. Moritz, S. Levine, M. Jordan, and P. Abbeel, “High-dimensional continuous control using generalized advantage estimation,” _arXiv preprint arXiv:1506.02438_ , 2015.