Maxwell’s equations in heterogeneous and time-varying mediaD. P. San Roman
Alerigi, D. I. Ketcheson, and B. S. Ooi

# A high-order finite volume method for Maxwell’s equations in heterogeneous
and time-varying media

Damian P. San Roman Alerigi  Photonics Laboratory, King Abdullah University of
Science and Technology _(KAUST)_ , Thuwal 21534, Saudi Arabia (). Now at EXPEC
Advanced Research Center, Aramco, Dhahran 31311, Saudi Arabia).
<EMAIL_ADDRESS>David I. Ketcheson Numerical Mathematics Group,
King Abdullah University of Science and Technology _(KAUST)_ , Thuwal 21534,
Saudi Arabia<EMAIL_ADDRESS>Boon S. Ooi Photonics
Laboratory, King Abdullah University of Science and Technology _(KAUST)_ ,
Thuwal 21534, Saudi Arabia<EMAIL_ADDRESS>
###### Abstract

We develop a finite volume method for Maxwell’s equations in materials whose
electromagnetic properties vary in space and time. We investigate both
conservative and non-conservative numerical formulations. High-order methods
are employed to accurately resolve fine structures that develop due to the
varying material properties. Numerical examples demonstrate the effectiveness
of the proposed method in handling temporal variation and its efficiency
relative to traditional 2nd-order FDTD.

###### keywords:

fvm, finite-volume, heterogeneous, spacetime, time-varying, photonics,
maxwell, electromagnetic, riemann, clawpack, emclaw

35Q61, 35Q15, 35L05, 78A25, 78A48, 78M12, 65Z05

##  1 Introduction

Maxwell’s equations in a charge- and current-free space are defined as:

(1a) $\displaystyle\partial_{t}\vec{D}-\nabla\times\vec{H}$ $\displaystyle=0,$
(1b) $\displaystyle\partial_{t}\vec{B}+\nabla\times\vec{E}$ $\displaystyle=0;$

where $\vec{E},\vec{B}:\mathbb{R}^{3}\times\mathbb{R}^{3+1}\to\mathbb{R}^{3}$
are the electric and magnetic fields, respectively. They are related to the
material fields $\vec{D}\in\mathbb{R}^{3}$ and $\vec{H}\in\mathbb{R}^{3}$ via
the constitutive relations:

(2a) $\displaystyle\vec{D}$
$\displaystyle=\vec{D}\left(\bar{\varepsilon}(\vec{r},t),\vec{E}\right),$ (2b)
$\displaystyle\vec{B}$
$\displaystyle=\vec{B}\left(\bar{\mu}(\vec{r},t),\vec{H}\right).$

where $\bar{\varepsilon},\bar{\mu}:\mathbb{R}^{3\times
3}\times\mathbb{R}^{3+1}\to\mathbb{R}^{3\times 3}$ are symmetrical and second-
rank tensors with non-zero off-diagonal entries, and describe the electric and
magnetic response of the material; they are known as permittivity and
permeability.

In this work, we are interested in the general case in which
$\bar{\varepsilon},\bar{\mu}$ vary in space and time, and $\vec{D}$, $\vec{B}$
may be nonlinear.

We are motivated by recent interest in materials whose properties can change
at a pace comparable to that of electromagnetic wave oscillation; _e.g._[5, 6,
23]. The effect of spatially heterogeneous materials on wave propagation has
received much attention, but the effect of temporal material variations is
much less studied. Novel effects arising from wave propagation in space-time-
varying media include trapping, confinement, ultra-short pulse coupling, beam
transformation, negative refraction, eigen–polarization, and the optical
Bohm–Aharonov effect [11, 20, 19, 22, 15, 16, 21, 2, 3, 1]. This has opened
the theoretical possibility to control light in new ways, _e.g._ Alcubierre’s
warp drives [25].

A numerical scheme for electromagnetic wave propagation in time-dependent
linear and heterogeneous media is developed in [8]. The algorithm is based on
the finite-difference time-domain (FDTD) method presented in [10], with an
extension to account for the temporal variation of the medium. This is
attained by applying an operator splitting method and Magnus’ series expansion
[8]. Due to the splitting, the scheme requires that certain iteration matrices
be recomputed at every time step.

We apply the high-order wave propagation method developed in [14]. This method
is based on weighted essentially non-oscillatory (WENO) reconstruction, Runge-
Kutta time integration, and wave-propagation Riemann solvers. We expect these
high-order algorithms to be more efficient than traditional second-order
finite differences (i.e., the FDTD method) for high-frequency waves in
rapidly-varying media. Furthermore, they can handle nonlinear media, do not
require operator splitting, and avoid the costly computation of iteration
matrices required in [8].

The rest of the paper is organized as follows. In Section 2, we consider two
forms of Maxwell’s equations with explicit time-dependence: first as a
homogeneous and conservative hyperbolic system in which the conserved
quantities appear implicitly, and second as a non-conservative system of
balance laws. Section 3 we review the high-order wave propagation method
introduced in [14] and present the general algorithm to solve the Riemann
problem. In Section 4, we derive the order of convergence for the algorithm
using some fundamental 1D and 2D examples. Finally, in Section 5, we discuss
future work, advantages, and opportunities.

##  2 Maxwell’s equations in nonlinear, time-varying media

To illustrate the difficulties involved in discretizing this system, let us
consider the one-dimensional case. Then $D,B,H,E$ are scalar fields, and the
system may be written in the form:

(3a) $\displaystyle\partial_{t}{D}+\partial_{x}{H}$ $\displaystyle=0,,$ (3b)
$\displaystyle\partial_{t}{B}+\partial_{x}{E}$ $\displaystyle=0.,$

In general, the constitutive relations in (2) can take any shape. For the
purpose of the following example, we follow a customary electromagnetic
approach and expand them in a power series [9]:

(4a) $\displaystyle D\left(\varepsilon(x,t),E\right)$
$\displaystyle=\varepsilon_{o}\left(E+\sum_{i}\chi_{e}^{(i)}(x,t)E^{i}\right),$
(4b) $\displaystyle B\left(\mu(x,t),H\right)$
$\displaystyle=\mu_{o}\left(H+\sum_{i}\chi_{h}^{(i)}(x,t)H^{i}\right),.$

where $\varepsilon_{o},\mu_{o}\in\mathbb{R}$ are the vacuum permittivity and
permeability, and $\chi^{(n)}\in\mathbb{R}$ are the electric (e) and magnetic
(h) susceptibility scalar functions.

Equation 4 allows us to consider two types of constitutive relations. For
linear media, we have

(5a) $\displaystyle D$ $\displaystyle=\varepsilon E,$ (5b) $\displaystyle B$
$\displaystyle=\mu H;$

whereas for non-linear meadia we can consider a system with cubic
nonlinearity, _i.e._ :

(6a) $\displaystyle D$ $\displaystyle=\varepsilon E+\chi^{3}E^{3},$ (6b)
$\displaystyle B$ $\displaystyle=\mu H+\chi^{3}H^{3}.$

Finally, we observe that equation 3 can be rewritten in conservative and non-
conservative form.

###  2.1 Conservative form

Let us define the conserved vector $\vec{q}$and flux
$\vec{f}\left(\vec{q}\right)$as:

(7) $\displaystyle\vec{q}$ $\displaystyle=\begin{pmatrix}D\\\ B\end{pmatrix},$
$\displaystyle\vec{f}\left(\vec{q}\right)$
$\displaystyle=\begin{pmatrix}H(\mu(x,t),B)\\\
E(\varepsilon(x,t),D)\end{pmatrix};$

hence we can summarize (3) as an homogeneous first-order hyperbolic system:

(8) $\partial_{t}\vec{q}+\partial_{x}\vec{f}(\vec{q};x,t)=0,$

with

(9) $\displaystyle\partial_{q}\vec{f}$
$\displaystyle=\begin{pmatrix}0&\partial_{B}H\\\
\partial_{D}E&0\end{pmatrix},$

and wave speeds ($s$) given by:

(10) $s=\pm\sqrt{\partial_{D}E\partial_{B}H}.$

####  2.1.1 Linear media

In the linear example we substitute (5) into (7), and obtain:

(11) $\displaystyle\vec{f}(\vec{q};x,t)$
$\displaystyle=\begin{pmatrix}B/\mu(x,t)\\\ D/\varepsilon(x,t)\end{pmatrix}.$

Consequently, the wave speeds are:

(12) $s=\pm\sqrt{1/(\varepsilon\mu)}.$

####  2.1.2 Nonlinear media

In the nonlinear case, we must determine derivatives of the inverse functions
$H(\mu,B)$ and $E(\varepsilon,D)$ to compute the speed of the waves in (10)
and other relevant quantities. We can approximate the terms $\partial_{B}H$
and $\partial_{D}E$ by means of an iterative method. For example, consider the
nonlinear constitutive relation for the electric field given in (6), and solve
for $E(D)$ as:

(13) $E(D)=\frac{1}{\varepsilon}\left(D-\chi^{3}(E(D))^{3}\right).$

Substituting the latter into itself yields

(14)
$E(D)=\frac{1}{\varepsilon}\left(D-\chi^{3}D^{3}+{\mathcal{O}}(D^{5})\right).$

Thus, finally we obtain the partial derivative of $E$ with respect to $D$ as:

(15)
$\partial_{D}E\approx\frac{1}{\varepsilon}\left(1-3\chi^{3}_{e}D^{2}\right).$

Similarly, we can follow the same algorithm to calculate the derivative of $H$
with respect to $B$ and obtain:

(16) $\partial_{B}H\approx\frac{1}{\mu}\left(1-3\chi^{3}_{h}B^{2}\right).$

So the wave speeds are:

(17)
$s\approx\pm\sqrt{\frac{1}{\varepsilon\mu}\left(1-3\chi^{3}_{e}D^{2}\right)\left(1-3\chi^{3}_{h}B^{2}\right)}+{\mathcal{O}}(D^{3})$

###  2.2 Non-conservative form

As we can see from the previous treatment, the conservative form can be
intricate to solve complex nonlinear materials. An alternative is to derive a
non-conservative form to (3). In this case, we proceed by applying the chain
rule to the time derivatives in (3), obtaining:

(18a) $\displaystyle\partial_{E}D\partial_{t}E+\partial_{x}H$
$\displaystyle=-\partial_{\varepsilon}D\partial_{t}\varepsilon,$ (18b)
$\displaystyle\partial_{H}B\partial_{t}H+\partial_{x}E$
$\displaystyle=-\partial_{\mu}B\partial_{t}\mu.$

Let us re-define the conserved quantity vector and flux as:

(19) $\displaystyle\vec{q}$ $\displaystyle=\begin{pmatrix}E\\\
H\end{pmatrix},$ $\displaystyle\vec{f}\left(\vec{q}\right)$
$\displaystyle=\begin{pmatrix}H\\\ E\end{pmatrix};$

and introduce a second-rank tensor known as the capacity, $\bar{\kappa}$, and
vector term named source, $\vec{\psi}$ term:

(20) $\displaystyle\bar{\kappa}(\vec{q};x,t)$
$\displaystyle=\begin{pmatrix}\partial_{E}D\\\ \partial_{H}B\end{pmatrix},$
$\displaystyle\vec{\psi}(\vec{q};x,t)$
$\displaystyle=\begin{pmatrix}-E\partial_{t}\varepsilon\\\
-H\partial_{t}\mu\end{pmatrix}.$

Then, we can rewrite equation 18 as a balance law of the form:

(21)
$\bar{\kappa}(\vec{q};x,t)\cdot\partial_{t}\vec{q}(x,t)+\partial_{x}\vec{f}\left(\vec{q}\right)=\vec{\psi}(\vec{q};x,t),$

This balance law results in a simpler form of linear and nonlinear cases. For
linear media, substitute (5) into (18) to produce:

(22a) $\displaystyle\varepsilon(x,t)E_{t}+\partial_{x}H$
$\displaystyle=-E\partial_{t}\varepsilon,$ (22b)
$\displaystyle\mu(x,t)H_{t}+\partial_{x}E$ $\displaystyle=-H\partial_{t}\mu.$

Whereas for nonlinear media we can substitute (6) instead and obtain:

(23a)
$\displaystyle\left(\varepsilon+2\chi^{3}_{e}E^{2}\right)\partial_{t}E+\partial_{x}H$
$\displaystyle=-E\partial_{t}\varepsilon,$ (23b)
$\displaystyle\left(\mu+2\chi^{3}_{h}H^{2}\right)\partial_{t}H+\partial_{x}E$
$\displaystyle=-H\partial_{t}\mu.$

Note that in either case the speed of the waves is given

(24) $s=\pm 1$

##  3 Semi-discrete wave propagation

Numerical wave propagation methods for systems similar to (21) have been
developed previously [18, 17, 14]. However, in those methods, the capacity
function $\kappa$ was assumed to depend only on $x$; whereas here, it may
depend also on $q$ and $t$. In this Section, we review those methods and
discuss their extension to the case of (8) and (21). Specifically, we extend
the scheme presented in [14] to handle time-varying fluxes and capacity
functions that depend on $q,x,t$. Since that scheme is based on the method of
lines, time-varying coefficients can be handled in a straightforward way. We
focus first on the one-dimensional method and then briefly describe the
extension to more spatial dimensions.

###  3.1 Basic scheme

We consider a grid of cells with centers $x_{i}$ and interfaces
$x_{i\pm\frac{1}{2}}$ (see figure 1), and define the average of $q$ over cell
$i$:

(25) $Q_{i}=\int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}q(x,t)dx.$

$x_{i-1/2}$$x_{i+1/2}$$i$$i-1$$i+1$ Figure 1: Schematic representation of 1D
uniform cell grid.

From the cell averages we construct a piecewise-polynomial approximation to
$q$:

(26) $\hat{\vec{q}}(x)=\hat{\vec{q}}_{i}\left(x\right),\qquad\text{for}\quad
x\in\left(x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right),$

where

(27)
$\hat{\vec{q}}_{i}\left(x\right)=\vec{q}\left(x,t\right)+\mathcal{O}\left(\Delta
x^{r+1}\right).$

At each interface $x_{i-\frac{1}{2}}$, the approximation is discontinuous;
hence, we define the interface values as:

(28) $\displaystyle q^{L}_{i}-\frac{1}{2}$
$\displaystyle=\hat{q}_{i-1}(x_{i}-\frac{1}{2})$ $\displaystyle
q^{R}_{i}-\frac{1}{2}$ $\displaystyle=\hat{q}_{i}(x_{i}-\frac{1}{2}).$

####  3.1.1 Conservative scheme

When applied to (8), the semi-discrete scheme presented in [14] takes the form

(29) $\frac{\partial Q_{i}}{\partial t}=-\frac{1}{\Delta
x}\left(\mathcal{A}^{+}\Delta q_{i}-\frac{1}{2}+\mathcal{A}^{-}\Delta
q_{i}+\frac{1}{2}+f(q^{L}_{i}+\frac{1}{2})-f(q^{R}_{i}-\frac{1}{2})\right).$

The flux difference arises from integrating the flux term over the cell. The
terms $\mathcal{A}^{-}\Delta q,\mathcal{A}^{+}\Delta q$ are referred to as
fluctuations and incorporate the effects of waves emanating from the cell
interfaces, where we must solve a Riemann problem.

At each cell interface, we approximate the solution of a Riemann problem with
initial states given by the reconstructed values at the interface:

(30) $q(x)=\left\\{\begin{matrix}q^{L}_{i}-\frac{1}{2}&x<x_{i-\frac{1}{2}},\\\
q^{R}_{i}-\frac{1}{2}&x>x_{i-\frac{1}{2}}.\end{matrix}\right.$

The solution is approximated by a sequence of jump discontinuities; this
amounts to decomposing the jump as:

(31)
$q^{R}_{i}-\frac{1}{2}-q^{L}_{i}-\frac{1}{2}=\sum_{p=1}^{m}\alpha^{p}_{i}-\frac{1}{2}r^{p}_{i}-\frac{1}{2}=\sum_{p}\mathcal{W}^{p}_{i}-\frac{1}{2}.$

Here the vectors $r^{p}$ are the eigenvalues of some approximate flux Jacobian
at the interface. The vectors $\mathcal{W}^{p}$ are known as waves, and the
collective effect of the right- and left-going waves yields the fluctuations:

(32a) $\displaystyle\mathcal{A}^{-}\Delta q_{i}-\frac{1}{2}$
$\displaystyle\equiv\sum_{p}\left(s^{p}_{i}-\frac{1}{2}\right)^{-}\mathcal{W}^{p}_{i}-\frac{1}{2},$
(32b) $\displaystyle\mathcal{A}^{+}\Delta q_{i}-\frac{1}{2}$
$\displaystyle\equiv\sum_{p}\left(s^{p}_{i}-\frac{1}{2}\right)^{+}\mathcal{W}^{p}_{i}-\frac{1}{2},$

Here $s^{p}$ is the speed of the wave $\mathcal{W}^{p}$ (often given by an
eigenvalue of an approximate flux Jacobian), and

$\left(s^{p}\right)^{-}=\min\left(s^{p},0\right)\qquad\left(s^{p}\right)^{+}=\max\left(s^{p},0\right).$

####  3.1.2 Non-conservative scheme

Following a similar derivation to its conservative counterpart, the non-
conservative form (21) becomes

(33) $\frac{\partial Q_{i}}{\partial t}=-\frac{1}{\bar{K}_{i}\Delta
x}\left(\mathcal{A}^{+}\Delta q_{i-\frac{1}{2}}+\mathcal{A}^{-}\Delta
q_{i+\frac{1}{2}}+f(q^{L}_{i}+\frac{1}{2})-f(q^{R}_{i}-\frac{1}{2})\right)+\Psi_{i}.$

Here $K_{i}$ and $\Psi_{i}$ are averages of $\kappa$ and $\psi$ over cell $i$.
For $\Psi$ we use the trapezoidal rule:

(34) $\displaystyle\Psi^{1}_{i}$
$\displaystyle=-\frac{\dot{\epsilon}^{R}_{i-\frac{1}{2}}E^{R}_{i}-\frac{1}{2}+\dot{\epsilon}^{L}_{i}+\frac{1}{2}E^{L}_{i}+\frac{1}{2}}{2}$
(35) $\displaystyle\Psi^{2}_{i}$
$\displaystyle=-\frac{\dot{\mu}^{R}_{i-\frac{1}{2}}H^{R}_{i}-\frac{1}{2}+\dot{\mu}^{L}_{i}+\frac{1}{2}H^{L}_{i}+\frac{1}{2}}{2}$

where for readability $\dot{\eta}=\eta_{t}$. For $\kappa$, we also use a
trapezoidal average. Let $\tilde{\eta}_{i}$ be a suitable average of $\eta$
across the $i$th cell at some time $t^{n}$. Then, at the interface
$x_{i-\frac{1}{2}}$ we can approximate the average value $K_{i-\frac{1}{2}}$
using the trapezoidal rule,

(36)
$K_{k,i-\frac{1}{2}}^{n}\approx\frac{\tilde{\eta}^{R}_{i-\frac{1}{2}}+\tilde{\eta}^{L}_{i-\frac{1}{2}}}{2}+\frac{p}{2^{p}}\>\chi^{(p)}\left(q^{R}_{k,i-\frac{1}{2}}+q^{L}_{k,i-\frac{1}{2}}\right)^{p-1},\quad
p\geq 2,$

###  3.2 Riemann solver

Following (LeVeque, 2002), we can design a Riemann solver for this system as
follows. We assume the values of $\epsilon,\mu$ are constant in each cell so
that we can write $\epsilon_{i},\mu_{i}$. The eigenvector matrix for $f_{q}$
is

(37) $\displaystyle R_{i}-\frac{1}{2}$
$\displaystyle=\begin{pmatrix}-Z_{i-1}&Z_{i}\\\ 1&1\end{pmatrix},$

where $Z_{i}=\sqrt{H_{B}/E_{D}}$. We can use an $f$-wave solver by solving the
system

(38) $R\beta=f_{i}(Q_{i})-f_{i-1}(Q_{i-1}).$

###  3.3 Reconstruction and time integration

It remains to specify the method of time integration and spatial
reconstruction. These ingredients determine the order of accuracy of the
scheme.

We use the ten-stage fourth-order strong-stability-preserving Runge-Kutta
scheme described in [12] for time integration. Although it uses many stages,
this method has a large region of absolute stability, allowing the use of
large CFL numbers.

We use fifth-order component-wise WENO reconstruction in space as described in
[24]. Traditionally, methods based on finite differences use several points
per wavelength to correctly model the high-frequency waves that occur in some
physical problems, _e.g._ electromagnetism. We use WENO or any high-degree
polynomial because it allows us to keep a high level of accuracy with
relatively coarser grids.

All of this is implemented in Clawpack’s PyClaw module [13, 4]. The Riemann
solvers for Maxwell’s equations are implemented in a separate package called
EMClaw, available at [7].

###  3.4 Algorithm

At every Runge-Kutta stage, the numerical implementation of (33) follows the
steps:

  1. 1.

For the non-conservative system only: update cell averages of the time-
dependent capacity $\bar{K}_{i}^{n}$, and source $\Psi_{i}^{n}$ functions;

  2. 2.

Reconstruct the interface values $q_{i-\frac{1}{2}}^{R},q_{i+\frac{1}{2}}^{L}$
using fifth-order WENO reconstruction;

  3. 3.

Solve the Riemann problem with initial states
$(q_{i-\frac{1}{2}}^{L},q_{i-\frac{1}{2}}^{R})$ to compute the fluctuations,
$\mathcal{A}^{\pm}\Delta q_{i-\frac{1}{2}}$;

  4. 4.

Calculate the flux difference
$f(q^{L}_{i}+\frac{1}{2})-f(q^{R}_{i}-\frac{1}{2})$;

  5. 5.

Compute $\partial Q_{i}/\partial t$ using the semi-discrete scheme (33).

###  3.5 Extension to two dimensions

Using a dimension-by-dimension approach, we extend the numerical wave
propagation method to two dimensions where the main equation is (on the
Cartesian grid):

(39)
$\bar{\kappa}\cdot\vec{q}_{t}+\vec{f}(\vec{q})_{x}+\vec{g}(\vec{q})_{y}=\vec{\psi}(\vec{q},x,y,t).$

In two dimensions, the semi-discrete scheme (33) takes the form:

(40) $\begin{split}\frac{\partial Q_{ij}}{\partial
t}=-\frac{1}{\bar{K}_{ij}\Delta x\Delta y}\Bigg{(}&\mathcal{A}^{+}\Delta
q_{i-\frac{1}{2},j}+\mathcal{A}^{-}\Delta
q_{i+\frac{1}{2},j}+f(q^{L}_{i+\frac{1}{2},j})-f(q^{R}_{i-\frac{1}{2},j})\\\
&+\mathcal{B}^{+}\Delta q_{i,j-\frac{1}{2}}+\mathcal{B}^{-}\Delta
q_{i,j+\frac{1}{2}}+g(q^{L}_{i,j+\frac{1}{2}})-g(q^{R}_{i,j-\frac{1}{2}})\Bigg{)}\\\
&+\Psi_{ij}.\end{split}$

The fluctuation terms are determined by solving the Riemann problem in the
corresponding direction and initial data; for example, $\mathcal{B}^{+}\Delta
q_{i,j-\frac{1}{2}}$ is calculated by solving the Riemann problem in the
$y$-direction with initial states
$(q^{L}_{i,j-\frac{1}{2}},q^{R}_{i,j-\frac{1}{2}})$. As noted in [14], for the
method to be high-order accurate, the fluctuations and flux differences must
be computed based on a multidimensional reconstruction of $q$. In the present
work, we instead use a much cheaper dimension-by-dimension reconstruction,
described in [14]. Although this leads to a formally 2nd-order accurate
scheme, it is still much less dissipative than traditional second-order
schemes and allows high-frequency waves to be resolved on much coarser grids.

##  4 Numerical applications and convergence

We test the above-described techniques numerically by studying the propagation
of electromagnetic waves in media with spacetime variations. We choose
examples where the exact solution can be calculated with great accuracy.

To obtain the results in this section, we developed EMClaw [7], a multi-
dimensional numerical solver for Maxwell’s equations in spacetime-varying and
nonlinear media, based on the scheme (33). EMClaw was designed as an extension
of the Clawpack package [4].

### Miscellaneous definitions

We use the numerical examples to measure the $L_{1}$ errors, $\vec{E}$, and
rate of convergence, $p$, of the numerical schemes developed in previous
sections. Specifically, we compute

(41) $E_{L_{1}}=\Delta x\sum_{i}|Q_{i}-\hat{Q_{i}}|,$

where $\hat{Q_{i}}$ is a highly accurate solution cell average computed by
characteristics (exact solution) or by using grid refinement.

For the subsequent numerical studies, let $E_{e}(h)$ denote the error of the
numerical solution with respect to the exact solution for some grid spacing
$h$, and similarly, let us denote by $E_{s}(h)$ the error of the numerical
solution on a grid with spacing $h$ with respect to the numerical solution
obtained on the refined grid, $h/2$.

Let $I$ denote the 2-norm of $\vec{q}$, the magnitude of the field,

(42) $I(\vec{r},t)=\sqrt{\sum_{i=1}^{m}q_{i}^{2}(\vec{r},t)}.$

Similarly, let $n$ denote the magnitude of the material profile

(43) $n=\sqrt{\eta_{e}\cdot\eta_{h}}.$

Finally let $F=I,n$, then we can define the maximum
$F_{max}(t)=\max_{\vec{r}}{F(\vec{r},t)}$, and the path described by it
$\vec{r}_{F}(t)=\\{\vec{r}:F(\vec{r},t)=F_{max}\\}$.

###  4.1 1D Electromagnetic

Here we apply scheme (33) to one-dimensional Maxwell’s equations in the case
of spacetime-varying media with linear and nonlinear components. We assume the
materials to be isotropic and piecewise homogeneous, leading to the 1D system
(18).

For the examples below, we set the initial condition to be the right-moving
pulse

(44) $E(x,0)=H(x,0)=\exp{\left(\frac{x-x_{o}}{\sigma}\right)^{2}}.$

Reflecting boundary conditions are used in all 1D examples.

####  4.1.1 Time-varying medium

We consider 1D electromagnetic waves in a linear, spatially homogeneous, and
time-varying medium. Namely, we solve (18) in the interval $x\in[0,100]$, with

$\varepsilon=\mu=1+\delta\eta\sin\left(\frac{12\pi}{100}t\right),$

where $\delta\eta\in\mathbb{R}$ is the oscillation amplitude.

The initial condition is given by (44) with $x_{o}=10$ and $\sigma=2$.

Figure 2: Time-staggered plot of $I$ for 1D wave propagation in a vibrating
medium with $n_{\omega}=12$.

Table 1 shows the error and convergence rates for propagation in this problem
at different grid sizes $h$ defined by the number of cell grids $m_{x}$.

|  | Conservative | Non-Conservative
---|---|---|---
$m_{x}$ | $h$ | $E_{e}(h)$ | $p_{e}$ | $E_{e}(h)$ | $p_{e}$
128 | 7.812e-01 | 2.814e-01 | 3.657 | 2.821e-01 | 3.597
256 | 3.906e-01 | 2.230e-02 | 3.904 | 2.331e-02 | 5.519
512 | 1.953e-01 | 1.490e-03 | 2.457 | 5.085e-04 | 5.713
1024 | 9.766e-02 | 2.714e-04 | 5.819 | 9.694e-06 | 5.265
2048 | 4.883e-02 | 4.809e-06 | 3.820 | 2.521e-07 | 5.078
4096 | 2.441e-02 | 3.405e-07 | 6.160 | 7.464e-09 | 5.028
Table 1: Errors and convergence rate for 1D electromagnetic wave propagation
in a linear, spatially homogeneous, and time-varying medium. Comparison
between the Conservative and Non-Conservative Riemann solvers, with $cfl=2.4$.

(a) Conservative (b) Non-Conservative

Figure 3: Convergence of the finite volume solution to the analytic solution,
$E_{e}(h)$, for 1D electromagnetic wave propagation in a linear, spatially
homogeneous and time-varying medium (blue-dashed); comparison between the
conservative and non-conservative Riemann solvers with $cfl=2.4$. Linear fit
with slope (red) $p_{e}\approx-5.173$ (non-conservative) $p_{e}\approx-3.922$
(conservative).

Figure 3 shows $E_{e}(h)$ for this case as a function of $m_{x}=L/h$, where
$L=100$ is the length of the simulation space. Using a least-square method, we
can approximate the slope of the line in the asymptotic convergence region of
$E_{e}(h)$ and find $p_{e}\approx-3.922$ for the conservative, and
$p_{e}\approx 5.173$ for the non-conservative, Riemann solvers.

####  4.1.2 Spacetime-varying medium

We now consider 1D electromagnetic waves in a linear and spacetime-varying
medium. Namely, we solve 1D equation (3) in the interval $x\in[0,300]$, with

(45)
$\eta=\eta_{o}+\delta\eta\exp{\left(\frac{x-x_{o}-v\>t}{\sigma}\right)^{2}},$

where where $\eta_{o}\in\mathbb{R}$ is the background material parameter,
$\delta\eta\in\mathbb{R}$ is the spacetime-variation’s amplitude,
$x_{o}\in\mathbb{R}$ is the offset, $v\in\mathbb{R}$ is the velocity, and
$\sigma\in\mathbb{R}$ is the full-width half maximum.

To be specific, we set $\eta_{o}=1.5$, $\delta\eta=0.15$, $x_{o}=25$,
$v=0.59$, and $\sigma=5$; for the initial conditions set $q_{o}$ with
$x_{o}=10$ and $\sigma=2$.

Figure 4: Time-staggered plots of $I$ (black) and $n$ (green) for 1D wave
propagation in a Gaussian-like spacetime-varying medium.

Table 2 shows the errors and convergence rates for propagation in this
problem.

|  | Conservative | Non-Conservative
---|---|---|---
$m_{x}$ | $h$ | $E_{e}(h)$ | $p_{e}$ | $E_{e}(h)$ | $p_{e}$
128 | 7.812e-01 | 3.648e+00 | 1.369 | 3.775e+00 | 1.365
256 | 3.906e-01 | 1.413e+00 | 1.411 | 1.465e+00 | 1.407
512 | 1.953e-01 | 5.315e-01 | 1.575 | 5.524e-01 | 1.544
1024 | 9.766e-02 | 1.784e-01 | 2.677 | 1.895e-01 | 2.658
2048 | 4.883e-02 | 2.790e-02 | 2.559 | 3.002e-02 | 4.185
4096 | 2.441e-02 | 4.735e-03 | 1.941 | 1.650e-03 | 5.217
Table 2: Errors and convergence rate for 1D wave propagation in a Gaussian-
like spacetime-varying medium. Comparison between the Conservative and Non-
Conservative Riemann solvers, with $cfl=2.4$

(a) Conservative (b) Non-Conservative

Figure 5: Convergence of the finite volume solution to the analytic solution,
$E_{e}(h)$, for 1D electromagnetic wave propagation in a linear, spatially
homogeneous and time-varying medium (blue-dashed); comparison between the
conservative and non-conservative Riemann solvers with $cfl=2.4$. Linear fit
with slope (red) $p_{e}\approx-2.782$ (non-conservative) $p_{e}\approx-2.311$
(conservative).

Figure 5 shows $E_{e}(h)$ for 1D wave propagation in a Gaussian-like
spacetime-varying medium as a function of $m_{x}=L/h$, where $L=300$ is the
length of the simulation space. Using a least-square method, we can
approximate the slope of the line in the asymptotic convergence region of
$E_{e}(h)$ and find $p_{e}\approx 3.481$.

To explain the difference in the rate of convergence between the 1D wave
propagation problems, recall that in the former, the material is only time-
dependent, which together with the results would suggest that the discrepancy
is due to the second-order approximation to $\Psi_{i}$ and $K_{i}$ within each
cell $i$ for spacetime-varying media. Note that because we use the trapezoidal
rule, we would expect the convergence rate to be second-order accurate; the
higher degree of convergence ensues from using high-order WENO interpolation
and Runge-Kutta methods.

####  4.1.3 Nonlinear flowing medium

In this section, we study a hybrid material that incorporates the spacetime-
varying medium of the previous section and a background nonlinear material,
$\chi^{(3)}\neq 0$. This problem is difficult to study in the paraxial
approximation, and a numerical scheme is usually required.

The material parameters and domain are the same in this example as in the
previous section. The nonlinearity is introduced by setting
$\bar{\chi}^{(3)}=0.1$. The initial condition $q_{o}$ is again (44) with
$x_{o}=10$ and $\sigma=2$.

Figure 6 (a) shows the effect of the nonlinearity in the absence of the moving
perturbation (45). Note that as the pulse evolves, the nonlinearity in the
medium compresses the pulse where $q_{x}\geq 0$ and expands it where
$q_{x}\leq 0$; that is, we observe the formation of a shock and rarefaction,
respectively.

(a) Nonlinear and stationary medium (b) Nonlinear and spacetime-varying medium

Figure 6: Time staggered plot of $I$ for 1D wave propagation in a nonlinear
and homogeneous medium (left), and a nonlinear spacetime-varying medium
(right).

Allow us now to introduce the flowing material of equation 45 and observe the
evolution of the right-moving pulse (44) plotted in Figure 6 (b). The results
suggest that the moving perturbation balances the dispersion introduced by the
nonlinearity; in other words, the former causes the front of the pulse to
expand while the moving perturbation compresses it.

To better appreciate this dynamic, let us look at the time rate of change of
$I_{max}$ plotted in Figure 7, notice that as the pulse interacts with the
perturbation, it begins to compress and increase in amplitude, while at the
same time, the nonlinearity causes rarefaction; ultimately, this leads to
equilibrium, as exemplified by the fact that $\frac{dI_{max}}{dt}\approx 0$.

(a) $\chi^{(3)}=0$ (b) $\chi^{(3)}=0.1$

Figure 7: Time rate of change of $I_{max}$ for 1D wave propagation in a
nonlinear Gaussian-like spacetime-varying medium.

###  4.2 2D Electromagnetic

Here we apply scheme (40) to the two-dimensional electromagnetic problem in
the case of linear spacetime-varying media. We assume the materials to be
isotropic and piecewise homogeneous, leading to the 2D system:

(46a)
$\displaystyle\kappa_{1}(q_{1},x,t)\left(q_{1}\right)_{t}-\left(q_{3}\right)_{y}$
$\displaystyle=\psi_{e}(q_{1},x,t),$ (46b)
$\displaystyle\kappa_{2}(q_{2},x,t)\left(q_{2}\right)_{t}+\left(q_{3}\right)_{x}$
$\displaystyle=\psi_{e}(q_{2},x,t),$ (46c)
$\displaystyle\kappa_{3}(q_{3},x,t)\left(q_{3}\right)_{t}+\left(q_{2}\right)_{x}-\left(q_{1}\right)_{y}$
$\displaystyle=\psi_{h}(q_{3},x,t);$

with initial condition to be the right-moving pulse

(47a) $\displaystyle q_{0}(x,0)$ $\displaystyle=0.0,$ (47b) $\displaystyle
q_{1}(x,0)$ $\displaystyle=g(y,y_{o})\>q_{o}(x,x_{o},\sigma),$ (47c)
$\displaystyle q_{2}(x,0)$ $\displaystyle=g(y,y_{o})\>q_{o}(x,x_{o},\sigma),$

where $\vec{q}_{o}$ is defined following (44), and $g:\mathbb{R}\to\mathbb{R}$
defines the transversal profile of the pulse. To be specific, we choose

$g(y)=\cos{\frac{\left(y-y_{o}\right)\pi}{L_{y}}}$

where $L_{y}$ is the simulation length in the $y$ direction.

We solve (46) in the interval $\\{x,y\\}\in[0,180]\times[0,180]$ with varying
grid space and set the all the boundary conditions to be reflecting or wall.

For the sake of simplicity, we again consider propagation in a system where
$\eta_{l}^{0}=1$ and the material has unitary impedance, _i.e._
$\eta_{k}=\eta$. Namely, the material profile is

(48)
$\varepsilon=\mu=\eta_{o}+\delta\eta\exp{\left(\left(\frac{x-x_{o}-v\>t}{\sigma_{x}}\right)^{2}+\left(\frac{y-y_{o}}{\sigma_{y}}\right)^{2}\right)}$

where where $\eta_{o}\in\mathbb{R}$ is the background material parameter,
$\delta\eta\in\mathbb{R}$ is the spacetime-variation’s amplitude,
$\square_{o}\in\mathbb{R}$ is the offset in the $x$ or $y$ direction; and
$\sigma_{i}\in\mathbb{R}$ is the full-width half maximum in the corresponding
direction. Specifically we set $\eta_{o}=1.5$, $\delta\eta=0.15$, $v=0.59$,
$x_{o}=25$, $y_{o}=\frac{L_{y}}{2}$, $\sigma_{x}=5$ and $\sigma_{y}=25$.

Table 3 summarizes the errors and convergence rate for subsequent grid
refinement at different grid sizes $h$ defined by the number of cell grids
$m_{x}$.

$mx/my$ | $h/mx$ | $E_{s}(h)$ | $p_{s}$
---|---|---|---
128 | 1.978 | 1.981e-01 | 2.231
256 | 4.944e-01 | 4.222e-02 | 2.980
512 | 1.236e-01 | 5.349e-03 | 2.210
1024 | 3.090e-02 | 1.156e-03 | 2.008
2048 | 7.725e-03 | 2.875e-04 | 2.004
4096 | 1.931e-03 | 7.168e-05 | —
Table 3: Errors and rate of convergence for 2D electromagnetic wave
propagation in a linear spacetime-varying medium.

Figure 8 shows $E_{s}(h)$. Again, using a least square linear fit, we find the
convergence rate in the asymptotic region to be $p_{s}\approx 2.524$.

Figure 8: Convergence of the finite volume solution to subsequent refined
finite volume solution, $E_{s}(h)$, 2D electromagnetic wave propagation in a
linear spacetime-varying medium (blue-dashed); and liner fit with slope
$p_{s}\approx-2.524$ (red).

##  5 Summary

We have shown an approach to extend the highly accurate wave propagation
algorithm of SharpClaw to model wave propagation in nonlinear and spacetime-
varying media in one and two spatial dimensions. The scheme is second-order
accurate, as demonstrated by the test results.

Using WENO and strong-stability-preserving time integration, high-order
accurate results are obtained in one dimension even when the average value of
the coefficients in each cell is second-order-accurate.

A drawback of our implementation is that to achieve the high-order convergence
in two dimensions, as observed in SharpClaw, the spacetime-varying
coefficients need to be resolved with high accuracy, for example, using high-
degree polynomial interpolation and adequate quadrature.

## References

  * [1] F. Belgiorno, S. Cacciatori, G. Ortenzi, L. Rizzi, V. Gorini, and D. Faccio, Dielectric black holes induced by a refractive index perturbation and the Hawking effect, Physical Review D, 83 (2011), http://dx.doi.org/10.1103/PhysRevD.83.024015, http://prd.aps.org/abstract/PRD/v83/i2/e024015.
  * [2] F. Belgiorno, S. L. Cacciatori, G. Ortenzi, V. G. Sala, and D. Faccio, Quantum Radiation from Superluminal Refractive-Index Perturbations, Physical Review Letters, 104 (2010), p. 140403, http://dx.doi.org/10.1103/PhysRevLett.104.140403, http://link.aps.org/doi/10.1103/PhysRevLett.104.140403http://prl.aps.org/abstract/PRL/v104/i14/e140403, arXiv:arXiv:0910.3508v3.
  * [3] S. L. Cacciatori, F. Belgiorno, V. Gorini, G. Ortenzi, L. Rizzi, V. G. Sala, and D. Faccio, Spacetime geometries and light trapping in travelling refractive index perturbations, New Journal of Physics, 12 (2010), p. 095021, http://dx.doi.org/10.1088/1367-2630/12/9/095021, http://stacks.iop.org/1367-2630/12/i=9/a=095021%****␣maxwell-fv.bbl␣Line␣25␣****http://stacks.iop.org/1367-2630/12/i=9/a=095021?key=crossref.eb283ba4c42fb3331796b446bb8fcf27.
  * [4] Clawpack Development Team, Clawpack software, 2016, http://www.clawpack.org.
  * [5] B. J. Eggleton, Chalcogenide photonics: fabrication, devices and applications Introduction, Optics Express, 18 (2010), p. 26632, http://dx.doi.org/10.1364/OE.18.026632, http://www.opticsexpress.org/abstract.cfm?URI=oe-18-25-26632.
  * [6] B. J. Eggleton, B. Luther-Davies, and K. Richardson, Chalcogenide photonics, Nature Photonics, 5 (2011), pp. 141–148.
  * [7] EMClaw Development Team, EMClaw, 2016, https://github.com/MaxwellGEMS/emclaw.
  * [8] I. Faragó, Á. Havasi, and R. Horváth, Numerical solution of the Maxwell equations in time-varying media using Magnus expansion, Central European Journal of Mathematics, 10 (2011), pp. 137–149, http://dx.doi.org/10.2478/s11533-011-0074-3, http://www.springerlink.com/index/10.2478/s11533-011-0074-3.
  * [9] D. J. Griffiths, Introduction to Electrodynamics, Benjamin Cummings, 3rd ed., 1999, http://www.amazon.com/Introduction-Electrodynamics-3rd-David-Griffiths/dp/013805326X.
  * [10] Kane Yee, Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media, IEEE Transactions on Antennas and Propagation, 14 (1966), pp. 302–307, http://dx.doi.org/10.1109/TAP.1966.1138693, http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=1138693.
  * [11] A. Kaplan, Light-induced nonreciprocity, field invariants, and nonlinear eigenpolarizations, Optics letters, 8 (1983), pp. 560–562, http://www.opticsinfobase.org/abstract.cfm?id=8131.
  * [12] D. I. Ketcheson, Highly Efficient Strong Stability-Preserving Runge–Kutta Methods with Low-Storage Implementations, SIAM Journal on Scientific Computing, 30 (2008), pp. 2113–2136, http://dx.doi.org/10.1137/07070485X, http://epubs.siam.org/doi/abs/10.1137/07070485X.
  * [13] D. I. Ketcheson, K. Mandli, A. J. Ahmadia, A. Alghamdi, M. Q. de Luna, M. Parsani, M. G. Knepley, and M. Emmett, PyClaw: Accessible, Extensible, Scalable Tools for Wave Propagation Problems, SIAM Journal on Scientific Computing, 34 (2012), pp. C210–C231, http://dx.doi.org/10.1137/110856976, http://epubs.siam.org/doi/abs/10.1137/110856976.
  * [14] D. I. Ketcheson, M. Parsani, and R. J. LeVeque, High-Order Wave Propagation Algorithms for Hyperbolic Systems, SIAM Journal on Scientific Computing, 35 (2013), pp. A351–A377, http://epubs.siam.org/doi/abs/10.1137/110830320.
  * [15] U. Leonhardt, Optical conformal mapping., Science (New York, N.Y.), 312 (2006), pp. 1777–80, http://dx.doi.org/10.1126/science.1126493, http://www.ncbi.nlm.nih.gov/pubmed/16728596.
  * [16] U. Leonhardt and T. G. Philbin, General relativity in electrical engineering, New Journal of Physics, 8 (2006), pp. 247–247, http://dx.doi.org/10.1088/1367-2630/8/10/247, http://stacks.iop.org/1367-2630/8/i=10/a=247?key=crossref.4e433b9506a5936770423817211714f3.
  * [17] R. J. LeVeque, Finite volume methods for hyperbolic problems, Cambridge University Press, 1st editio ed., 2002, http://books.google.com/books?hl=en{&}lr={&}id=QazcnD7GUoUC{&}oi=fnd{&}pg=PR17{&}dq=Finite-volume+methods+for+hyperbolic+problems{&}ots=WYBwvhSL1m{&}sig=3nSTrIb-KR2rQKrRZqaUr1X33sYhttp://books.google.com/books?hl=en{&}lr={&}id=QazcnD7GUoUC{&}oi=fnd{&}pg=PR17{&}dq=Finite-volume+.
  * [18] R. J. LeVeque, Finite-volume methods for non-linear elasticity in heterogeneous media, International Journal for Numerical Methods in Fluids, 40 (2002), pp. 93–104, http://dx.doi.org/10.1002/fld.309, http://onlinelibrary.wiley.com/doi/10.1002/fld.309/abstracthttp://doi.wiley.com/10.1002/fld.309.
  * [19] P. Piwnicki and U. Leonhardt, Optics of moving media, Applied Physics B: Lasers and Optics, 72 (2001), pp. 51–59, http://www.springerlink.com/index/M329N4AN59WWU599.pdf.
  * [20] C. Ridgely, Applying covariant versus contravariant electromagnetic tensors to rotating media, American Journal of Physics, 67 (1999), p. 414, http://link.aip.org/link/?AJPIAS/67/414/1.
  * [21] G. Rousseaux, On the electrodynamics of Minkowski at low velocities, EPL (Europhysics Letters), 84 (2008), p. 20002, http://dx.doi.org/10.1209/0295-5075/84/20002, http://stacks.iop.org/0295-5075/84/i=2/a=20002?key=crossref.3f68616dc01be811c3c319258f3ff434.
  * [22] N. Rozanov and G. Sochilin, Geometrical optics of moving media, Optics and spectroscopy, 98 (2005), pp. 441–446, http://www.springerlink.com/index/WG1J7274U1Q472M5.pdf.
  * [23] D. P. San-Roman-Alerigi, D. H. Anjum, Y. Zhang, X. Yang, A. B. Slimane, T. K. Ng, M. N. Hedhili, M. Alsunaidi, and B. S. Ooi, Electron irradiation induced reduction of the permittivity in chalcogenide glass (As2S3) thin film, Journal of Applied Physics, 113 (2013), p. 044116, http://dx.doi.org/10.1063/1.4789602, http://link.aip.org/link/JAPIAU/v113/i4/p044116/s1{&}Agg=doi.
  * [24] C.-W. Shu, High Order Weighted Essentially Nonoscillatory Schemes for Convection Dominated Problems, SIAM Review, 51 (2009), pp. 82–126, http://dx.doi.org/10.1137/070679065, http://epubs.siam.org/doi/abs/10.1137/070679065.
  * [25] I. Smolyaninov, Metamaterial-based model of the Alcubierre warp drive, Physical Review B, 84 (2011), pp. 2–5, http://dx.doi.org/10.1103/PhysRevB.84.113103, http://link.aps.org/doi/10.1103/PhysRevB.84.113103.


$\displaystyle=\int_{0}^{t}\sum_{\begin{subarray}{c}k\in\mathbb{Z}\\\
k_{1}+k_{2}+k_{3}=k\end{subarray}}e^{ikx}e^{i(t-t^{\prime})k^{4}}e^{it^{\prime}(k_{1}^{4}-k_{2}^{4}+k_{3}^{4})}\hat{\phi}(k_{1})\overline{\hat{\phi}(-k_{2})}\hat{\phi}(k_{3})n_{3}(k_{1},k_{2},k_{3})dt^{\prime}$
$\displaystyle=\sum_{\begin{subarray}{c}k\in\mathbb{Z}\\\
k_{1}+k_{2}+k_{3}=k\end{subarray}}e^{ikx+itk^{4}}\hat{\phi}(k_{1})\overline{\hat{\phi}(-k_{2})}\hat{\phi}(k_{3})n_{3}(k_{1},k_{2},k_{3})\int_{0}^{t}e^{-it^{\prime}(k^{4}-k_{1}^{4}+k_{2}^{4}-k_{3}^{4})}dt^{\prime}.$

Here $n_{3}(k_{1},k_{2},k_{3})=(k_{1}+k_{2})^{2}+\frac{3}{2}(k_{1}+k_{3})^{2}$
is the symbol corresponding to the terms in $N_{3}$. We may now choose
$\hat{\phi}(k)=k^{-s}(\delta_{k,N}+\delta_{k,N_{0}})$, where $N_{0}\ll N$. The
choice of the $N_{0}$ parameter is not important as long as it is, say, fixed.
For all further calculations the reader may assume $N_{0}=1$. We then observe
$\|\phi\|_{\hat{H}^{s}_{r}}\sim 1$ independent of the two parameters.

Inserting this into the above expression we note that it suffices to look at
the terms that produce a resulting frequency of $k=N$. There are three such
choices for the tuple $(k_{1},k_{2},k_{3})$, namely $(N,-N,N)$,
$(N,-N_{0},N_{0})$ and $(N_{0},-N_{0},N)$. Note that for each of these three
choices the resonance relation $k^{4}-k_{1}^{4}+k_{2}^{4}-k_{3}^{4}$ cancels
and the integral in the formula above is equal to $t$ and the symbol of our
non-linearity has size on the order of $N^{2}$.

These frequency choices thus produce Fourier coefficients (at frequency $N$)
on the order of $tN^{2-3s}$ (for the first one) and $tN^{2-s}$ (for the second
and third). The remaining five frequency constellations cannot cancel these
contributions as they are of lower order in $N$.

This leaves us with the following lower bound for the Sobolev norm of the
operator that is the derivative of the flow:

$\left\|\left.\frac{\partial^{3}u}{\partial\delta^{3}}\right|_{\delta=0}\right\|_{\hat{H}^{s}_{r}}^{r^{\prime}}\gtrsim
N^{sr^{\prime}}\cdot
t^{r^{\prime}}N^{(2-s)r^{\prime}}(1+N^{-2sr^{\prime}})\geq
t^{r^{\prime}}N^{2r^{\prime}}$

for $1<r\leq\infty$. If $r=1$ we still have a lower bound of $tN^{2}$ though
with a simpler argument. Letting $N\to\infty$ we can now see, that the flow
cannot be $C^{3}$ for any $s\in\mathbb{R}$. ∎

The previous proposition shows that an approach with (just) a fixed-point
theorem to prove well-posedness must fail at any regularity in
$\hat{H}^{s}_{r}(\mathbb{T})$. As is stated in Theorem 3.8 the situation is
much more dire at lower regularities. Its proof lies in the following
proposition.

######  Proposition 6.4.

Let $j\in\mathbb{N}$, $1\leq r\leq\infty$ and $s<j-1$. The flow
$S:\hat{H}^{s}_{r}(\mathbb{T})\times(-T,T)\to\hat{H}^{s}_{r}(\mathbb{T})$ of
the Cauchy problem

$i\partial_{t}u+(-1)^{j+1}\partial_{x}^{2j}u=|u|^{2}\partial_{x}^{2j-2}u\quad\text{with}\quad
u(t=0)=u_{0}\in\hat{H}^{s}_{r}(\mathbb{T}),$ (6.5)

cannot be uniformly continuous on bounded sets.

We want to point out, that equation (6.5) is in fact a higher-order NLS-like
equation according to (2.10). More so it even fits the structure of an NLS
hierarchy equation (2.8), though it is unlikely to be one because of its
simple non-linearity.

###### Proof of Proposition 6.4.

We follow a similar argument to the one used in, for example,
[OhTzvetkov2017]*Appendix A.2.

The reader may verify that our equation (6.5) has the two-parameter family of
solutions

$u_{N,a}(x,t)=N^{-s}a\exp(i(Nx-N^{2j}t+N^{2j-2-2s}|a|^{2}t)).$

We fix $a\in\mathbb{R}$ at two different values and will only deal with the
two solution families $u_{n}(x,t)=u_{N_{n},1}(x,t)$ and
$\tilde{u}_{n}(x,t)=u_{N_{n},1+\frac{1}{n}}(x,t)$ depending on
$n\in\mathbb{N}$. $N_{n}$ will be chosen later. We find that

$\|u_{n}(\cdot,0)\|_{\hat{H}^{s}_{r}},\|\tilde{u}_{n}(\cdot,0)\|_{\hat{H}^{s}_{r}}\lesssim
1\qquad\text{and}\qquad\|u_{n}(\cdot,0)-\tilde{u}_{n}(\cdot,0)\|_{\hat{H}^{s}_{r}}\sim\frac{1}{n},$

where the implicit constant is independent of $n\in\mathbb{N}$. Now choosing

$t_{n}=\frac{\pi N_{n}^{2s+2-2j}}{(1+\frac{1}{n})^{2}-1}$

and $N_{n}$ large enough, such that $t_{n}\leq\frac{1}{n}$, we may then
observe that

$\|u_{n}(\cdot,t_{n})-\tilde{u}_{n}(\cdot,t_{n})\|_{\hat{H}^{s}_{r}}=\left|\exp(iN_{n}^{2j-2-2s}(1-(1+\frac{1}{n})^{2})t_{n})-(1+\frac{1}{n})\right|=2+\frac{1}{n}.$

Letting $n\to\infty$ this shows that the flow is not uniformly continuous.
Such a choice is possible, if $2s+2-2j<0$ or equivalently $s<j-1$ as stated. ∎

##  Appendix A The first few NLS hierarchy equations

For the reader’s convenience and future reference we will list the first few
conserved quantities $I_{k}$ derived from (2.5) and their associated non-
linear evolution equations (2.7) in terms of the potentials $q$ and $r$. In
this form both the focusing and defocusing variants of the (NLS) hierarchy can
be derived by the identifications $r=+\overline{q}$ or $r=-\overline{q}$
respectively.

Though we will not just give the even numbered equations, corresponding to the
NLS hierarchy, but also those corresponding to the mKdV hierarchy. Using the
identification $r=q$ one arrives at the real mKdV hierarchy discussed in
[AGTowers]. Deriving _a_ complex mKdV hierarchy (of which again there is a
defocusing and focusing variant) is also possible (again using the
identifications $r=\pm\overline{q}$). But there are two problems:

  1. (1)

Identifying $r=\pm\overline{q}$ for the equation induced by $I_{4}$, see
(A.1), does not lead to the well known form of the complex mKdV equation given
in (1.2). Rather the non-linearity is replaced by $\pm 6|u|^{2}\partial_{x}u$,
up to a choice of $\alpha_{3}$. For our local well-posedness theory this does
not make a difference, as we are able to estimate both non-linearities equally
well. Though for a treatment relying more on the structure of the equation
(e.g. for cancellation properties) this may be a relevant difference.

When looking at the real mKdV hierarchy, i.e. using $r=q$, this problem does
not present itself.

  2. (2)

If one wishes to use the identification $r=-\overline{q}$ the compatibility
condition for the coefficients $\alpha_{2j+1}$ reads
$\alpha_{2j+1}=-\overline{(\alpha_{2j+1})}$, as in (2.6). Meaning
$\alpha_{2j+1}$ is imaginary121212It is non-zero, as otherwise this would lead
to a trivial equation. and this introducing a factor $i$ that is usually not
present in complex mKdV-like equations.

Again, looking at the real mKdV hierarchy this is a non-issue, see also
[Alberty1982-i]*Section 3.2.2.

Not choosing an identification $r=\pm\overline{q}$ or $r=q$ also has the
advantage, that we may derive the equations in the KdV hierarchy by setting
$r=-1$, see [Alberty1982-i]*Section 3.2.1

Finally we note that our conserved quantities may differ from those given
elsewhere in the literature, as these are only determined up to (repeated)
partial integration and simplification. The equations though only differ up to
a choice of $\alpha_{k}$.

A similar listing is given in [Koch2018]*Appendix C and
[KochKlaus2023]*Appendix C.

  1. (1)

$n=1,2$. Phase shifts & Group of translations

$\displaystyle I_{1}=-\frac{1}{2i}\int qr\operatorname{d\\!}{x}$ and
$\displaystyle I_{2}=-\left(\frac{1}{2i}\right)^{2}\int
qr_{x}\operatorname{d\\!}{x}$ $\displaystyle q_{t}=2\alpha_{0}q$ and
$\displaystyle q_{t}=i\alpha_{1}q_{x}$

  2. (2)

$n=3$. cubic non-linear Schrödinger equation

$\displaystyle I_{3}=\left(\frac{1}{2i}\right)^{3}\int
q_{x}r_{x}+q^{2}r^{2}\operatorname{d\\!}{x}$ $\displaystyle
q_{t}=\frac{\alpha_{2}}{2}(-q_{xx}+2q^{2}r)$

  3. (3)

$n=4$. modified Korteweg-de-Vries equation

$\displaystyle I_{4}=\left(\frac{1}{2i}\right)^{4}\int
q_{x}r_{xx}+qq_{x}r^{2}+4q^{2}rr_{x}\operatorname{d\\!}{x}$ $\displaystyle
q_{t}=\frac{-\alpha_{3}}{4}(q_{xxx}-6qq_{x}r)$ (A.1)

  4. (4)

$n=5$. fourth order NLS hierarchy equation

$\displaystyle I_{5}=\left(\frac{1}{2i}\right)^{5}\int-
q_{xx}r_{xx}+q_{xx}r^{2}+6qq_{x}rr_{x}+5q^{2}r_{x}^{2}+6q^{2}rr_{xx}-2q^{3}r^{3}\operatorname{d\\!}{x}$
$\displaystyle
q_{t}=\frac{-\alpha_{4}}{8}(-q_{xxxx}+8qq_{xx}r+2q^{2}r_{xx}+4qq_{x}r_{x}+6q_{x}^{2}r-6q^{3}r^{2})$

  5. (5)

$n=6$. fifth order mKdV hierarchy equation

$\displaystyle I_{6}=\left(\frac{1}{2i}\right)^{6}\int$ $\displaystyle-
qr_{xxxxx}+qq_{xxx}r^{2}+8qq_{xx}rr_{x}+11qq_{x}r_{x}^{2}+12qq_{x}rr_{xx}$
$\displaystyle+18q^{2}r_{x}r_{xx}+8q^{2}rr_{xxx}-6q^{2}q_{x}r^{3}-16q^{3}r^{2}r_{x}\operatorname{d\\!}{x}$
$\displaystyle q_{t}=\frac{i\alpha_{5}}{2^{4}}(q_{xxxxx}$
$\displaystyle-10qq_{xxx}r-10qq_{xx}r_{x}-10qq_{x}r_{xx}-20q_{x}q_{xx}r-10q_{x}^{2}r_{x}$
$\displaystyle+30q^{2}q_{x}r^{2})$

  6. (6)

$n=7$. sixth order NLS hierarchy equation

$\displaystyle I_{7}=\left(\frac{1}{2i}\right)^{7}\int$ $\displaystyle-
qr_{xxxxxx}+qq_{xxxxx}r^{2}+10qq_{xxx}rr_{x}+19qq_{xx}r_{x}^{2}$
$\displaystyle+52qq_{x}r_{x}r_{xx}+20qq_{xx}rr_{xx}+20qq_{x}rr_{xxx}+19q^{2}r_{xx}^{2}$
$\displaystyle+28q^{2}r_{x}r_{xxx}+10q^{2}rr_{xxxx}+5q^{4}r^{4}$
$\displaystyle-6qq_{x}^{2}r^{3}-8q^{2}q_{xx}r^{3}-64q^{2}q_{x}r^{2}r_{x}-50q^{3}rr_{x}^{2}-30q^{3}r^{2}r_{xx}\operatorname{d\\!}{x}$
$\displaystyle q_{t}=\frac{\alpha_{6}}{2^{5}}(-q_{xxxxxx}$
$\displaystyle+12qq_{xxxx}r+2q^{2}r_{xxxx}+18qq_{xxx}r_{x}+22qq_{xx}r_{xx}+8qq_{x}r_{xxx}$
$\displaystyle+30q_{x}q_{xxx}r+20q_{x}^{2}r_{xx}+20q_{xx}^{2}r+50q_{x}q_{xx}r_{x}+20q^{4}r^{3}$
$\displaystyle-20q^{3}rr_{xx}-50q^{2}q_{xx}r^{2}-10q^{3}r_{x}^{2}-60q^{2}q_{x}rr_{x}-70qq_{x}^{2}r^{2})$

## References