text:
\begin{document}

\title[Improved FDTD Method for Maxwell's Equations]
{An Improved Finite Difference Time Domain Method for Two
Dimensional Maxwell's Equations}

\author{Timothy Meagher$^{1}$, Bin Jiang$^{2}$, Peng Jiang$^{3}$}

\footnotetext[{1}]{Department of Mathematical Sciences, University
of Nevada, Las Vegas, NV 89154, USA}

\footnotetext[{2}]{Fariborz Maseeh Department of Mathematics and
Statistics, Portland State University, Portland OR 97201, USA}

\footnotetext[{3}]{Department of Chemical Engineering, University of
Florida, Gainesville, FL 32611, USA}

\bibliographystyle{unsrt}

\begin{abstract} An efficient finite-difference time-domain (FDTD) algorithm is built to solve the transverse electric 2D Maxwell's equations
with inhomogeneous dielectric media where the electric fields are discontinuous across the dielectric interface. The new algorithm is derived
based upon the integral version of the Maxwell's equations as well as the relationship between the electric fields across the interface. It is an improvement
over the contour-path effective-permittivity algorithm by including some extra terms in the formulas. The scheme is validated in solving the scattering
of a dielectric cylinder with exact solution from Mie theory and is then compared with the above contour-path method, the usual staircase and
the volume-average method. The numerical results demonstrate that the new algorithm has achieved significant improvement in accuracy over
the other methods. Furthermore, the algorithm has a simple structure and can be merged into any existing FDTD software package very easily.
\end{abstract}

\keywords{finite difference time domain; Maxwell's equations;
effective permittivity}

\maketitle

\section{Introduction}

The finite-difference time-domain (FDTD) algorithm is one of the
most popular numerical methods to solve Maxwell's equations, as
proposed by Yee \cite{Yee1966}. It is a time marching forward method
which easily makes a visual representation of the fields. It also
scales very well since the number of computation required is
proportional to the size of the model and the method requires no
large-scale linear algebra computation \cite{Taflove2005}. The FDTD
method requires a completely structured grid also known as the Yee
grid. However, Many real world problems have complicated geometries
and numerical inaccuracy happens when the modeled objects do not fit
the grid very well.

In this paper, we focus on two-dimensional transverse electric (TE)
Maxwell's equation within multiple non-magnetic media where the
electric fields are not continuous across the media interface so the
numerical schemes need to be designed to take care of such
discontinuity. The other case is the transverse magnetic equation
where the electric field is always continuous across the media
interface and the original algorithm will work as in the homogeneous
medium.

The original FDTD algorithm proposed in Yee's paper assigns each of
the field component an electric permittivity and a magnetic
permeability solely based on the material properties at its
location. This is commonly known as the staircasing or pixelation
method. In general, the staircasing method has an error that scales
with $\mathcal{O}(\Delta x^2)$ given a grid size of $\Delta x$ for
cells that are homogenous. However, in cells that contain a medium
interface, the local error becomes $\mathcal{O}(1)$. Even though the
number of cells that contain an interface is often a small fraction
of the total cells, these local errors can cause a rough global
error of $\mathcal{O}(\Delta x)$.

Subsequent research in modern optics and electromagnetism dealing
with the Maxwell's equations have focused on more and more complex
geometries with different materials so that the complicated
dielectric tensor can cause serious pixelation issues. In turn,
there have been many improvements in FDTD algorithms to overcome
these pixelation issues. One way to overcome the pixelation problem
is by fitting the mesh to the device, thereby generating the
so-called non-orthogonal FDTD \cite{Harms1992, Liu2009}. The second
method is to refine the Cartesian mesh around the interfaces and is
thus called subgridding method \cite{Zivanovic1991}. Even though
these two methods converge faster than the standard FDTD method,
their complexity are larger than original FDTD method and they also
have time stability issues.

In order to keep a completely structured grid and maintain all the
benefits of FDTD, another method redefines Maxwell's integral
equations around the interface of the media. This is known as the
effective-permittivity (EP) method. This method uses sub-pixel
smoothing techniques to change the permittivities of the field
components around the interface to produce better results by taking
account of many different factors, such as the interface conditions,
permittivities of adjacent field components and so on. The revised
effective permittivities smooth out the pixel error. The main
benefit of these methods is that there is a negligible increase of
the numerical load since all the revised permittivities have been
calculated before the main FDTD algorithm loop.

The first method of sub-pixel smoothing is the volume average
effective permittivity (V-EP) \cite{Dey1999}. This method assigns
effective permittivities for all the field components in a cell that
contains media interface. The effective permittivities are
calculated by taking a weighted average based on how much percent of
the volume each medium occupies. This method is simple for
implementation and also stabilizes the error fluctuation that occurs
in the staircasing method. However, it does not decrease the error
rapidly and its error can be worse than staircasing method sometimes
due to ignorance of the interface orientation.

The second method proposes dynamic formulas for the EP based on the
orientation of the media interface and their expressions are
accurate when the interface is perpendicular or parallel to the mesh
axes \cite{Kaneda1997, Hwang2001, Liuchang2012}. These EPs improve
the accuracy of the FDTD method, while keeping the simple structure
of the original algorithm. However, their performance will
deteriorate for a curved or flat interface not perpendicular or
parallel to the mesh axes.

The third method uses the reflection coefficients to derive
effective permittivities and shows the method can also achieve
second-order convergence with several special slanted angles between
the interface and the Yee grid rather than just orthogonal or
parallel to the Yee grid \cite{Hirono2000, Hirono2010}. However, its
derivation is quite complicated and has no ways to be extended to
arbitrary interfaces.

In a landmark paper \cite{Mohammadi2005}, Mohammadi et al. proposed
the so-called contour-path effective permittivities (CP-EP) method.
The CP-EP method first addresses the orthogonal and parallel cases
as shown above, then extends the idea to handle with more general
geometry. Unlike the previous methods, it could handle arbitrary
boundary with any orientation of the boundary. Compared with
traditional staircasing methods, CE-EP has almost no additional
runtime cost since determination of the effective permittivities has
been done before the main FDTD loop.

However, we will demonstrate that CP-EP does not incorporate one
important term from the interface conditions into the FDTD
algorithm. Then we will develop a revised EP algorithm where the
missing term will be re-considered in order to produce a more
accurate and stable algorithm. Since this method incorporates the
boundary conditions very well, it will be denoted as BC-EP.

In fact, such deficiency of CP-EP has also been identified in
\cite{Farjadpour2006, Oskooi2009} where effective dielectric
permittivity tensors are constructed to take care of both isotropic
and anisotropic materials so as to achieve second-order convergence
under the open source software MEEP. However, such proposed method
which satisfies the interface conditions for electromagnetic fields
has been shown to have late-time instabilities, and many
possibilities to average the effective dielectric tensor are
explored to avoid late time instabilities \cite{Werner2007}. Further
efforts have been reported in \cite{Bauer2011} to construct a new
second-order scheme by taking the average of eight triplets.
Although the new method is highly accurate, its effective dielectric
permittivity tensor can still be asymmetric thus unstable for
certain conditions. Therefore, the last improvement of the effective
dielectric tensor has been given in \cite{Werner2013} to make it
symmetric and stable. The numerical results show that this scheme
gives the best result in general and the error in practice still
lies in between first and second-order in most cases.

In the following, we will introduce the new BC-EP method by adding a
few terms to the established CP-EP method for the cells around the
interface so as to keep the numerical load increase as small as
possible and improve the accuracy significantly. Then the numerical
tests verify that BC-EP has a much better performance than CP-EP,
V-EP and staircasing methods while still maintaining numerical
stability. BC-EP has a very simple structure and can be merged into
any existing FDTD packages easily, compared with other established
FDTD software \cite{Oskooi2009, Werner2013} which can be run only
within their framework. Conclusion and some future research
consideration will be provided in the end.

\section{New Algorithm Design}

In this section, we will propose a new effective permittivity method
to solve Maxwell's equations in a domain composed of two different
non-magnetic media with different electric permittivities. To simply
the algorithm description, we will focus on the most important
transverse electric (TE) scenario. That is,
$\mathbf{E}(x,y)=E_x\mathbf{i}+E_y\mathbf{j}$ is located in the
incident xy-plane while $\mathbf{H}(x,y)=H_z\mathbf{k}$ is along
z-axis. Moreover, we restrict the discussion to dielectric and
non-magnetic media such that $\mathbf{D}(x,y)=\epsilon
\mathbf{E}(x,y)$ and $\mathbf{B}(x,y)= \mathbf{H} (x,y)$ where
$\epsilon$ stays as a constant in different media and
$\epsilon_0=\mu_0=1$. Based on the integral version of the Maxwell's
equations, we will derive the new algorithm for any orientation of
the interface. This is a significant improvement over all current
methods which are only accurate under specific angles between the
interface and Cartesian coordinate system, such as parallel or
orthogonal to Yee-axis.

As a preliminary step, let us set up the relation of
$\mathbf{D}(x,y)$ and $\mathbf{E}(x,y)$ across the interface of
these two media. It is well known that $\mathbf{H}(x,y)$ is
continuous in the whole domain. Meanwhile, across the media
interface, the tangential component $\mathbf{E}_\tau(x,y)$ of the
electric field and the normal component $\mathbf{D}_n (x,y)$ of the
electric flux are continuous as well, as shown in
Figure~\ref{interface} where $\mathbf{n}=(n_x, n_y)$ is the outward
unit norm vector from region 1.
\begin{figure}[h]
    \begin{center}
            \includegraphics[scale=0.2]{interface}
                \caption{Relation of $\mathbf{D}$ and $\mathbf{E}$ across interface $\Gamma$}
                \label{interface}
    \end{center}
\end{figure}
Choose $A \in \Gamma$. $\mathbf{D}_1 (A)=(D_{x1}(A), D_{y1}(A))$ and
$\mathbf{D}_2 (A)=(D_{x2}(A), D_{y2}(A))$ are the electric fluxes at
A from different media, and $\mathbf{E}_1(A)=(E_{x1}(A), E_{y1}(A))$
and $\mathbf{E}_2(A)=(E_{x2}(A), E_{y2}(A))$ are the corresponding
electric fields.

Suppose $\mathbf{n}=(n_x, n_y)$ is the unit norm vector at $A$ from
region 1, then the corresponding tangent unit vector is
$\mathbf{\tau}=(n_y, -n_x).$ Therefore, we can calculate $E_\tau$
and $D_n$ at point $A$ in different regions as follows.
\begin{align}
   E_{\tau_1}= \mathbf{E}_1 \cdot  \mathbf{\tau}=E_{x_1}n_y-E_{y_1}n_x
\end{align}
\begin{align}
   E_{\tau_2}= \mathbf{E}_2 \cdot  \mathbf{\tau}=E_{x_2}n_y-E_{y_2}n_x
\end{align}
Meanwhile,
\begin{align}
   D_{n_1}= \epsilon_1 \mathbf{E}_{n_1}=\epsilon_1 E_{x_1}n_x+ \epsilon_1E_{y_1}n_y
   \label{Dinterface}
\end{align}
\begin{align}
   D_{n_2}= \epsilon_2 \mathbf{E}_{n_2}=\epsilon_2 E_{x_2}n_x+ \epsilon_2 E_{y_2}n_y
\end{align}
Based on the fact that $E_{\tau_1}=E_{\tau_2}$ and
$D_{n_1}=D_{n_2}$, we have
\begin{align*}
    \begin{cases}
    E_{x_1}n_y-E_{y_1}n_x & =E_{x_2}n_y-E_{y_2}n_x\\
    \epsilon_1  E_{x_1}n_x+\epsilon_1E_{y_1}n_y &= \epsilon_2 E_{x_2}n_x+\epsilon_2 E_{y_2}n_y
    \end{cases}
\end{align*}
Therefore, by using $n_x^2+n_y^2=1$, we get
\begin{align*}
    \begin{bmatrix}
        E_{x_1}   \\
    E_{y_1}
     \end{bmatrix}=
      \begin{bmatrix}
        n_y &  - n_x    \\
        \epsilon_1 n_x   &  \epsilon_1  n_y
     \end{bmatrix} ^{-1}
    \begin{bmatrix}
        n_y &  - n_x    \\
        \epsilon_2 n_x   &  \epsilon_2  n_y
     \end{bmatrix}
     \begin{bmatrix}
        E_{x_2}   \\
    E_{y_2}
     \end{bmatrix}
\end{align*}
Or
\begin{align}
    \begin{cases}
    E_{x_1} & = (\frac{ \epsilon_2}{ \epsilon_1}n_x^2+n_y^2) E_{x_2} + (\frac{ \epsilon_2}{ \epsilon_1} -1)n_xn_yE_{y_2}\\
    E_{y_1} & = (\frac{ \epsilon_2}{ \epsilon_1} -1)n_xn_y E_{x_2}   + (n_x^2+\frac{\epsilon_2}{\epsilon_1}n_y^2)E_{y_2}
     \end{cases}
     \label{eintro}
\end{align}
Furthermore, by using $\mathbf{D}_i=\epsilon_i\mathbf{E}_i$ at $A$,
we obtain
 \begin{align}
    \begin{cases}
    D_{x_1} & = (\frac{\epsilon_1}{\epsilon_2}n_y^2+n_x^2) D_{x_2} + (1-\frac{\epsilon_1}{\epsilon_2})n_xn_y D_{y_2}\\
    D_{y_1} & = (1-\frac{\epsilon_1}{\epsilon_2})n_xn_y D_{x_2} + (\frac{\epsilon_1}{\epsilon_2}n_x^2+n_y^2) D_{y_2}
     \end{cases}
     \label{dintro}
\end{align}
If the interface $\Gamma$ is parallel or orthogonal to $x$-axis,
either $n_y=0$ or $n_x=0$ so $n_x \cdot n_y=0$. Then the relation
\eqref{eintro} involving $\mathbf{E}_1$ and $\mathbf{E}_2$ will be
simplified such that $E_{x_1}$ and $E_{y_1}$ will only depend on
$E_{x_2}$ or $E_{y_2}$, respectively. Same holds for $\mathbf{D}_1$
and $\mathbf{D}_2$ in relation \eqref{dintro}.

By using \eqref{dintro}, we can express $D_{x_2}$ in region 2 based
on its neighbor $D_{x_1}$ in region 1 and its corresponding
$D_{y_2}$ in the same region 2. Similar formula can also be derived
for $D_{y_2}$ based on $D_{y_1}$ and $D_{x_2}$.
\begin{align}
    D_{x_2} & = \frac{\epsilon_2}{\epsilon_2 n_x^2+\epsilon_1 n_y^2}D_{x_1}
                 + \frac{(\epsilon_1-\epsilon_2)n_xn_y}{\epsilon_2 n_x^2+\epsilon_1 n_y^2}D_{y_2} \label{dintroDx} \\
    D_{y_2} & = \frac{(\epsilon_1-\epsilon_2)n_xn_y}{\epsilon_1 n_x^2+\epsilon_2 n_y^2}D_{x_2}
                 + \frac{\epsilon_2}{\epsilon_1 n_x^2+\epsilon_2 n_y^2}D_{y_1} \label{dintroDy}
\end{align}
\eqref{eintro} will be used to derive the new scheme for Faraday's
law while \eqref{dintroDx} and \eqref{dintroDy} will be used to
derive the new scheme for Ampere's law.

Now let us express Ampere's law in integral form
\begin{align}
   \frac{\partial }{\partial t} \iint_S \mathbf{D} \cdot  \mathbf{n} dS= \oint_{\partial S}  \mathbf{H} \cdot d\mathbf{l}
   \label{Ampere's}
\end{align}
and Faraday's law integral form
\begin{align}
   \frac{\partial }{\partial t} \iint_S \mathbf{B} \cdot \mathbf{n} dS=- \oint_{\partial S}  \mathbf{E} \cdot d\mathbf{l}
   \label{Faraday's}
\end{align}

Firstly, let us discretize Ampere's law \eqref{Ampere's} around the
interface $\Gamma$ to update $E^n_x$ and $E^n_y$ as below.
\begin{figure}[ht!]
    \begin{center}
            \includegraphics[scale=0.2]{amex}
                \caption{Discretization of Ampere's law across $\Gamma$}
                \label{amex}
    \end{center}
\end{figure}
To update $E_x$, we choose $S$ in \eqref{Ampere's} to be the
rectangle in $yz$-plane whose projection in $xy$-plane is the line
segment $OC$. Suppose the interface $\Gamma$ intersect $OC$ at point
$A$ and $OA=f$ as in Figure~\ref{amex}. Applying time discretization
to Ampere's law along $S$ at time $t_{n+\frac{1}{2}}$ yields
\begin{align}
    \frac{\int_O^C D^{n+1}_{x_{(i,y)}}dy-\int_O^C D^{n}_{x_{(i,y)}}dy}{\Delta t}=H^{n+\frac{1}{2}}_{z_{(i,j)}}-H^{n+\frac{1}{2}}_{z_{(i,j-1)}}
    \label{newdisAmpere's}
\end{align}
Based on formula \eqref{dintroDx}, we can compute the integrals of
left hand side as follows. From now on, all the higher-order terms
of $dx$ and $dy$ are ignored for the clarification of the formulas.
\begin{align}
    \int_O^C D^{n+1}_{x_{(i,y)}}dy &= \int_O^A D^{n+1}_{x_{(i,y)}}dy+\int_A^C D^{n+1}_{x_{(i,y)}}dy \nonumber \\
    &= f D^{n+1}_{x_{(i,j-\frac{1}{2})}} + (\Delta y - f) D^{n+1}_{x_{(A,\Omega_2)}} \nonumber \\
    &= f D^{n+1}_{x_{(i,j-\frac{1}{2})}} + (\Delta y - f) \left[\frac{\epsilon_2}{\epsilon_2 n_x^2+\epsilon_1 n_y^2}D^{n+1}_{x_{(A,\Omega_1)}}
       + \frac{(\epsilon_1-\epsilon_2)n_xn_y}{\epsilon_2 n_x^2+\epsilon_1 n_y^2}D^{n+1}_{y_{(A,\Omega_2)}}\right] \nonumber \\
    &= \frac{f(\epsilon_2 n_x^2+\epsilon_1 n_y^2) + (\Delta y- f)\epsilon_2}{\epsilon_2 n_x^2+\epsilon_1 n_y^2} D^{n+1}_{x_{(i,j-\frac{1}{2})}}
       + \frac{(\epsilon_1-\epsilon_2)n_xn_y}{\epsilon_2 n_x^2+\epsilon_1 n_y^2}(\Delta y- f)D^{n+1}_{y_{(C)}}
    \label{newdisAmpere'sdn1}
\end{align}
Similarly,
\begin{align}
    \int_O^C D^{n}_{x(i,y)}dy & = \frac{f(\epsilon_2 n_x^2+\epsilon_1 n_y^2) + (\Delta y- f)\epsilon_2}{\epsilon_2 n_x^2+\epsilon_1 n_y^2}
    D^n_{x_{(i,j-\frac{1}{2})}} \nonumber \\
    & \quad + \frac{(\epsilon_1-\epsilon_2)n_xn_y}{\epsilon_2 n_x^2+\epsilon_1 n_y^2}(\Delta y- f)D^n_{y_{(C)}}
    \label{newdisAmpere'sdn}
\end{align}
Furthermore,
\begin{align}
    D^{n+1}_{y_{(C)}} & = D^n_{y_{(C)}}+\frac{\partial D_y}{\partial t}\rvert_C^{t_n} \cdot \Delta t \nonumber \\
                      & = D^n_{y_{(C)}}-\frac{\partial H_z}{\partial x}\rvert_C^{t_{n+\frac{1}{2}}} \cdot \Delta t \nonumber \\
     &= D^n_{y_{(C)}}-\left(\frac{H^{n+\frac{1}{2}}_{z_{(i,j)}}-H^{n+\frac{1}{2}}_{z_{(i-1,j)}}}{\Delta x}\right) \Delta t
     \label{newdisAmpere'sen1}
\end{align}
By putting \eqref{newdisAmpere'sdn1}, \eqref{newdisAmpere'sdn} and
\eqref{newdisAmpere'sen1} into \eqref{newdisAmpere's}, we obtain
\begin{align}
    \frac{f(\epsilon_2 n_x^2+\epsilon_1 n_y^2) + (\Delta y- f)\epsilon_2}{\epsilon_2 n_x^2+\epsilon_1 n_y^2}
    \left(D^{n+1}_{x_{(i,j-\frac{1}{2})}}-D^{n}_{x_{(i,j-\frac{1}{2})}}\right)
    =\Delta t\left(H^{n+\frac{1}{2}}_{z_{(i,j)}}-H^{n+\frac{1}{2}}_{z_{(i,j-1)}}\right) \nonumber \\
    + \frac{\Delta t}{\Delta x} \frac{(\epsilon_1-\epsilon_2)n_xn_y}{\epsilon_2 n_x^2+\epsilon_1 n_y^2}(\Delta y-f)
    \left(H^{n+\frac{1}{2}}_{z_{(i,j)}}-H^{n+\frac{1}{2}}_{z_{(i-1,j)}}\right)
\end{align}
which, based on $D_x=\epsilon_1 E_x$ in $\Omega_1$, becomes
\begin{align}
    E^{n+1}_{x_{(i,j-\frac{1}{2})}} &= E^{n}_{x_{(i,j-\frac{1}{2})}}
    + \Delta t \frac{\frac{\epsilon_2}{\epsilon_1} n_x^2+ n_y^2}{f(\epsilon_2 n_x^2+\epsilon_1 n_y^2) + (\Delta y- f)\epsilon_2}
                                        \left(H^{n+\frac{1}{2}}_{z_{(i,j)}}-H^{n+\frac{1}{2}}_{z_{(i,j-1)}}\right) \nonumber \\
    &+ \frac{\Delta t}{\Delta x} \frac{(\Delta y-f)(1-\frac{\epsilon_2}{\epsilon_1})n_xn_y}{f(\epsilon_2 n_x^2+\epsilon_1 n_y^2) + (\Delta y- f)\epsilon_2}
    \left(H^{n+\frac{1}{2}}_{z_{(i,j)}}-H^{n+\frac{1}{2}}_{z_{(i-1,j)}}\right)
\label{newdisAmpere'slaw1}
\end{align}
which shows improved discretization of Ampere's law for updating
$E_x.$

It should also be noted that the last term on the right hand side of
\eqref{newdisAmpere'slaw1} have been ignored by CP-EP for algorithm
simplification but it indeed provides necessary corrections for the
interface conditions across multiple media. Furthermore, the
numerical tests have also verified that in order to make the
algorithm for Ampere's Law more accurate, these two $H_z$ values in
the last term should stay close to the interface. That is, if the
normal direction of the interface $n$ satisfies $n_x\cdot n_y
> 0$ as shown in Figure~\ref{amex}, we use
$(H^{n+\frac{1}{2}}_{z_{(i,j)}}-H^{n+\frac{1}{2}}_{z_{(i-1,j)}})$
since these two $H_z$ points are closer to the interface, while for
the other $n_x\cdot n_y\leq 0$ case, we use
$(H^{n+\frac{1}{2}}_{z_{(i+1,j)}}-H^{n+\frac{1}{2}}_{z_{(i,j)}})$ in
\eqref{newdisAmpere'slaw1} instead since these two $H_z$ nodes are
closer to the interface than the other pair. This improvement is
feasible since $H_z$ in continuous across the interface.
\begin{figure}[ht!]
    \begin{center}
            \includegraphics[scale=0.2]{amey}
                \caption{Discretization of the integral along $OC$}
                \label{amey}
    \end{center}
\end{figure}

Similarly, update of $E_y$ at point $B=(i-\frac{1}{2},j)$ by
Ampere's law \eqref{Ampere's} can be taken care analogously as shown
in Figure~\ref{amey}, where the interface $\Gamma$ intersects the
line segment $OC$, projection of $S$ at $xy$-plane, at $A$. Let
$OA=d$. By repeating the same procedure as above while based on
formula \eqref{dintroDy}, we can also obtain the discretization of
Ampere's law for $E_y$ as below:
\begin{align}
    E^{n+1}_{y_{(i,j-\frac{1}{2})}} &= E^{n}_{y_{(i,j-\frac{1}{2})}}
    + \Delta t \frac{\frac{\epsilon_2}{\epsilon_1} n_y^2+ n_x^2}{d(\epsilon_1 n_x^2+\epsilon_2 n_y^2) + (\Delta x- d)\epsilon_2}
                               \left(H^{n+\frac{1}{2}}_{z_{(i-1,j)}}-H^{n+\frac{1}{2}}_{z_{(i,j)}}\right) \nonumber \\
    &+ \frac{\Delta t}{\Delta y} \frac{(\Delta x-d)(1-\frac{\epsilon_2}{\epsilon_1})n_xn_y}{d(\epsilon_1 n_x^2+\epsilon_2 n_y^2) + (\Delta x- d)\epsilon_2}
                               \left(H^{n+\frac{1}{2}}_{z_{(i,j-1)}}-H^{n+\frac{1}{2}}_{z_{(i,j)}}\right)
\label{newdisAmpere'slaw2}
\end{align}
where
$(H^{n+\frac{1}{2}}_{z_{(i,j-1)}}-H^{n+\frac{1}{2}}_{z_{(i,j)}})$ in
the last term will be replaced by
$(H^{n+\frac{1}{2}}_{z_{(i,j)}}-H^{n+\frac{1}{2}}_{z_{(i,j+1)}})$
when\textbf{} $n_x\cdot n_y\leq 0$ holds.

Secondly, let's discretize Faraday's law \eqref{Faraday's} around
the interface $\Gamma$.
\begin{figure}[ht!]
    \begin{center}
            \includegraphics[scale=0.2]{fahz}
                \caption{Discretization of Faraday's law across $\Gamma$}
                \label{fahz}
    \end{center}
\end{figure}
We choose $S$ to be the rectangle centered at $(i,j)$. If the medium
interface $\Gamma$ does not cut through interior of $S$, the
standard Yee scheme should apply without modification. So we only
consider the case where $\Gamma$ intersects with $S$. It's obvious
that $\Gamma$ will only intersect with at most two sides, e.g. the
bottom and left side as in Figure~\ref{fahz}. All other cases can be
handled in the same way. Similarly as for Ampere's law, all the
higher-order terms of dx and dy are ignored for the clarification of
the formulas. Discretization of the left hand side term of
\eqref{Faraday's} at the time $t_n$ via finite difference scheme
yields
\begin{align*}
   \frac{ \iint_S \mathbf{H}^{n+\frac{1}{2}} \cdot \mathbf{n} dS - \iint_S \mathbf{H}^{n-\frac{1}{2}} \cdot \mathbf{n} dS  }{\Delta t}
   = -\oint_{\partial S}  \mathbf{E}^n \cdot d\mathbf{l}
\end{align*}
or
\begin{align}
   \iint_S \mathbf{H}^{n+\frac{1}{2}} \cdot \mathbf{n} dS= \iint_S \mathbf{H}^{n-\frac{1}{2}} \cdot \mathbf{n} dS
   -\Delta t \oint_{\partial S}  \mathbf{E}^n \cdot d\mathbf{l}
   \label{newhupdate}
\end{align}
since $\mathbf{H}$ is always perpendicular to the incident
$xy$-plane and so is always continuous across the interface, then by
using midpoint rule,
\begin{align}
   \iint_S \mathbf{H}^{n+\frac{1}{2}} \cdot \mathbf{n} dS= H^{n+\frac{1}{2}}_{z_{(i,j)}} \Delta x \Delta y
   \label{hmid1}
\end{align}
\begin{align}
   \iint_S \mathbf{H}^{n-\frac{1}{2}} \cdot \mathbf{n} dS= H^{n-\frac{1}{2}}_{z_{(i,j)}} \Delta x \Delta y
   \label{hmid2}
\end{align}
However, the situation for $\mathbf{E}$ is much more complicated. In
Figure~\ref{fahz}, the line integral of $\mathbf{E}^n$ along the
right and top sides can be handled as Standard Yee scheme since they
are completely inside $\Omega_2$. But the line integrals along the
bottom and left sides need to be handled separately.
 \begin{align}
   \oint_{\partial S}  \mathbf{E}^n \cdot d\mathbf{l} &= \int_{bottom}  \mathbf{E}^n \cdot d\mathbf{l} +\int_{left}  \mathbf{E}^n \cdot d\mathbf{l}
   \nonumber \\   & \quad+\int^{(j+\frac{1}{2})\Delta y}_{(j-\frac{1}{2})\Delta y}E^n_{(i+\frac{1}{2},y)}dy-\int^{(i+\frac{1}{2})\Delta x}_{(i-\frac{1}{2})\Delta x}E^n_{(x,j+\frac{1}{2})}dx \nonumber \\
&=\int^C_O E^n_{(x,j-\frac{1}{2})}dx-\int^D_O E^n_{(i-\frac{1}{2},y)}dy +E^n_{y_{(i+\frac{1}{2},j)}}\Delta y -E^n_{x_{(i,j+\frac{1}{2})}} \Delta x
\label{cire}
\end{align}
Suppose the portion OA of the bottom side in $\Omega_1$ has length
$d$, and the portion OB of the left side in $\Omega_1$ has length
$f$ as in Figure~\ref{fahz}, then the first two terms on the right
hand side of \eqref{cire} can be handled as below. For the first
term on the right hand side of \eqref{cire}, we have
\begin{align}
    \int^C_O E^n_{(x,j-\frac{1}{2})}dx &= \int^A_O E^n_{(x,j-\frac{1}{2})}dx+\int^C_A E^n_{(x,j-\frac{1}{2})}dx \nonumber \\
    &= E^n_{x_{(i,j-\frac{1}{2})}} \cdot d+ E^n_x \rvert_{A, \Omega_2} \cdot (\Delta x -d)
    \label{rcire}
\end{align}
By using formula \eqref{eintro} with $\Omega_1$ and $\Omega_2$
interchanged, we have
\begin{align}
    E^n_x \rvert_{A, \Omega_2} = (\textstyle{\frac{ \epsilon_1}{ \epsilon_2}}n_x^2+n_y^2)E^n_x \rvert_{A, \Omega_1}+( \frac{ \epsilon_1}{ \epsilon_2} -1)n_xn_yE^n_y \rvert_{A, \Omega_1}
\end{align}
Putting it into \eqref{rcire} yields
\begin{align}
    \int^C_O E^n_{(x,j-\frac{1}{2})}dx &= \textstyle E^n_{x_{(i,j-\frac{1}{2})}} \cdot d+  (\frac{ \epsilon_1}{ \epsilon_2}n_x^2+n_y^2)(\Delta x -d)E^n_x \rvert_{A, \Omega_1} \nonumber \\
    & \textstyle + (\frac{\epsilon_1}{\epsilon_2} -1)n_xn_y(\Delta x -d)E^n_y \rvert_{A, \Omega_1} \nonumber \\
    & \textstyle = [d+(\Delta x-d)(\frac{ \epsilon_1}{\epsilon_2}n_x^2+n_y^2)]E^n_{x_{(i,j-\frac{1}{2})}} \nonumber \\
    & \textstyle + (\frac{\epsilon_1}{\epsilon_2} -1)n_xn_y(\Delta x -d)E^n_{y_{(i,j-\frac{1}{2})}}
    \label{ebottom}
\end{align}
where $E^n_{y_{(i,j-\frac{1}{2})}}$ is the value of $E^n_y$
evaluated at point $(i,j-\frac{1}{2})$ in $\Omega_1$. However, only
$E_x$ instead of $E_y$ is calculated at the mesh point
$(i,j-\frac{1}{2})$. Therefore, we will take the average value of
$E^n_y$ at the neighboring pair of mesh points $(i-\frac{1}{2},j-1)$
and $(i+\frac{1}{2},j)$ if they are both located within $\Omega_1$,
or take the average value of $E^n_y$ at the other neighboring pair
mesh points $(i-\frac{1}{2},j)$ and $(i+\frac{1}{2},j-1)$ if they
are both located within $\Omega_1$. If none of the pairs are located
within $\Omega_1$ , this term will be set to zero as in CP-EP.
Therefore, BC-EP will try to make corrections for the missing terms
in CP-EP if possible.

Similarly, we can calculate the second term on the right hand side
of \eqref{cire}. Since the $E_y$ node $(i-\frac{1}{2},j)$ is located
within $\Omega_2$ instead of $\Omega_1$, so the new $\mathbf{n}$
should be  negative of the old $\mathbf{n}$ used above. Therefore,
applying integration midpoint rule, Taylor expansion and formula
\eqref{eintro} where $(n_x,n_y)$ replaced by $(-n_x,-n_y)$ and $
\epsilon_1, \epsilon_2$ interchanged, we have
\begin{align}
    \int^D_O E^n_{(i-\frac{1}{2},y)}dy &= [f(n_x^2+\textstyle{\frac{ \epsilon_2}{ \epsilon_1}n_y^2)}+(\Delta  y-f)]E^n_{y_{(i-\frac{1}{2},j)}} \nonumber \\
    & \quad+( \frac{ \epsilon_2}{ \epsilon_1} -1)n_xn_yf E^n_{x_{(i-\frac{1}{2},j)}}
    \label{eleft}
\end{align}
when $E^n_{x_{(i-\frac{1}{2},j)}}$ will be taken care similarly as
$E^n_{y_{(i,j-\frac{1}{2})}}$ from the first term. By putting
\eqref{ebottom}, \eqref{eleft} into \eqref{cire}, we get
\begin{align}
       \oint_{\partial S} \mathbf{E}^n \cdot d\mathbf{l} =& \textstyle \ [d+(\Delta x-d)(\frac{ \epsilon_1}{ \epsilon_2}n_x^2+n_y^2)]E^n_{x_{(i,j-\frac{1}{2})}}
       +E^n_{y_{(i+\frac{1}{2},j)}}\Delta y -E^n_{x_{(i,j+\frac{1}{2})}} \Delta x \nonumber \\
       & \textstyle -[f(n_x^2+\frac{\epsilon_2}{\epsilon_1}n_y^2)+(\Delta  y-f)]E^n_{y_{(i-\frac{1}{2},j)}}
       + (\frac{\epsilon_1}{\epsilon_2} -1)n_xn_y(\Delta x -d) E^n_{y_{(i,j-\frac{1}{2})}} \nonumber \\
       & \textstyle - (\frac{\epsilon_2}{\epsilon_1} -1)n_xn_yf E^n_{x_{(i-\frac{1}{2},j)}}
        \label{newcire}
\end{align}
Finally, by putting \eqref{hmid1}, \eqref{hmid2} and \eqref{newcire}
into \eqref{newhupdate} and dividing both sides by $\Delta x \Delta
y$, we obtain the new discretization of Faraday's law:
\begin{align}
    H^{n+\frac{1}{2}}_{z_{(i,j)}} &=  \textstyle H^{n-\frac{1}{2}}_{z_{(i,j)}}  - \frac{\Delta t}{\Delta x \Delta y}\{[d+(\Delta x-d)
       (\frac{\epsilon_1}{\epsilon_2}n_x^2+n_y^2)]E^n_{x_{(i,j-\frac{1}{2})}}+E^n_{y_{(i+\frac{1}{2},j)}}\Delta y  \nonumber \\
       & \textstyle  \quad -E^n_{x_{(i,j+\frac{1}{2})}} \Delta x - [f(n_x^2+\frac{ \epsilon_2}{ \epsilon_1}n_y^2)+(\Delta  y-f)]E^n_{y_{(i-\frac{1}{2},j)}} \nonumber \\
       & \textstyle \quad + (\frac{\epsilon_1}{\epsilon_2} -1)n_xn_y(\Delta x -d)E^n_{y_{(i,j-\frac{1}{2})}}
                          - (\frac{\epsilon_2}{\epsilon_1} -1)n_xn_yf E^n_{x_{(i-\frac{1}{2},j)}} \}
\label{newdisFaraday'slaw}
\end{align}
where the last two terms are expressed by
\begin{align*}
E^n_{y_{(i,j-\frac{1}{2})}} & =
\begin{cases}
    \frac{1}{2}(E^n_{y_{(i-\frac{1}{2},j-1)}} + E^n_{y_{(i+\frac{1}{2},j)}}) & \mbox{ if } (i-\frac{1}{2},j-1)\in \Omega_1, \quad (i+\frac{1}{2},j)\in \Omega_1 \\
    \frac{1}{2}(E^n_{y_{(i-\frac{1}{2},j)}} + E^n_{y_{(i+\frac{1}{2},j-1)}}) & \mbox{ if } (i-\frac{1}{2},j)\in \Omega_1, \quad (i+\frac{1}{2},j-1)\in \Omega_1 \\
    0                                                                        & \mbox{ otherwise }
\end{cases}
\intertext{ and }
E^n_{x_{(i-\frac{1}{2},j)}} & =
\begin{cases}
    \frac{1}{2}(E^n_{x_{(i-1,j-\frac{1}{2})}} + E^n_{x_{(i,j+\frac{1}{2})}}) & \mbox{ if } (i-1,j-\frac{1}{2})\in \Omega_2, \quad (i,j+\frac{1}{2})\in \Omega_2 \\
    \frac{1}{2}(E^n_{x_{(i-1,j+\frac{1}{2})}} + E^n_{x_{(i,j-\frac{1}{2})}}) & \mbox{ if } (i-1,j+\frac{1}{2})\in \Omega_2, \quad (i,j-\frac{1}{2})\in \Omega_2 \\
    0                                                                        & \mbox{ otherwise }
\end{cases}
\end{align*}
If we set $\epsilon_1= \epsilon_2$ in the new algorithm
\eqref{newdisAmpere'slaw1}, \eqref{newdisAmpere'slaw2} for Amepre's
law and \eqref{newdisFaraday'slaw} for Faraday's law, we can
retrieve the standard FDTD scheme in homogeneous medium. Therefore,
BC-EP is a simple extension of the original schemes in order to take
care of the media interface. Meanwhile, the extra terms added in the
formulas will greatly improve the accuracy compared with CP-EP.
Furthermore, BC-EP has a very simple structure and can be merged
into any FDTD software by just revising the effective permittivities
and adding some extra terms when necessary.

It should be mentioned that due to complicated expressions of those
ignored higher-order terms in the newly derived formulas, it is
difficult to demonstrate the convergence order of BC-EP
theoretically. However, subsequent numerical results will show that
BC-EP can achieve the highest convergence (close to second order)
compared with other methods. Another point is that it is still an
open challenge to accurately handle objects with sharp corners,
where the resulting field singularities are known to degrade the
accuracy of all numerical schemes \cite{Oskooi2009}. Therefore, just
as other methods, BC-EP performs well for straight or curved
interfaces instead of interfaces with corners.

\section{Numerical Results}

In this section, we will demonstrate convergence order and stability
of BC-EP numerically. To accomplish this task, we implement the
algorithm to solve the Maxwell's equations on a dielectric cylinder,
together with the staircasing, V-EP, and CP-EP methods. All the
implementations are set up for the 2D FDTD TE case. The total
scattering cross sections (SCS) are calculated by all four methods
and then compared with the well-known analytic solution by Mie
Theory \cite{Bohren1983} so as to measure the accuracy of each
method.

\begin{figure}[ht!]
    \begin{center}
    \includegraphics[scale=0.5]{numericalsetup}
    \end{center}
    \caption{The numerical test setup using a dielectric cylinder, giving enough space for the total-field/scattered-field and integration line before the PML}
    \label{numericalsetup}
\end{figure}

The dielectric cylinder is simulated in an area that is $10$ times
the radius of the cylinder plus some extra space for PML condition.
The cylinder is centered at $((\frac{N_x}{2}+0.5)\Delta
x,(\frac{N_y}{2}+0.5)\Delta y)$ where $N_x$ and $N_y$ are the number
of the cells along the $x$ and $y$ directions, $\Delta x$ and
$\Delta y$ are the horizontal and vertical mesh sizes of each cell
as seen in Figure~\ref{numericalsetup}. We run the test for two
cases: one case for a larger cylinder with radius $r=400\textrm{nm}$
and another case for a smaller cylinder with $r=150\textrm{nm}$. The
total-field/scattered-field line is three times the length of the
cylinder radius $r$ away from the center, and the integration line
for calculating SCS is four times $r$ away from the center, while
the PML starts at five times $r$ away from the center.

The source wave is a planar wave in the $E_y$ direction with a
Gaussian envelope $e^{\frac{-(t-t_0)^2}{2}\cdot
\frac{1}{2\pi(c/\lambda_0-c/\lambda_1)^2}}\cdot \cos(2\pi ct)$,
where $t_0$ is six times $\Delta t$, $c$ is speed of light in a
vacuum, and $[\lambda_0, \lambda_1]$ is the testing range of
wavelengths. In our simulation, we calculate SCS for $601$
wavelengths ranging from $\lambda_0=400\textrm{nm}$ to
$\lambda_1=1000\textrm{nm}$ of visual light spectrum with equal
distance. Fast Fourier transform (FFT) is then applied to transform
the electric and magnetic fields from spatial domain into light
frequency domain so as to calculate SCS over all related
wavelengths.

We use square unit cells with $\Delta x= \Delta y$ to simplify the
simulation. The time step is set to $\Delta t= \frac{S \Delta x}{c}$
where $S=\frac{0.98}{ \sqrt{3}}$ to ensure stability\
\cite{Taflove2005}. Also for stability concerns, $\Delta x$ were
chosen to divide the smallest tested wave length
($\lambda_0=400\textrm{nm}$) by at least twenty times. The total SCS
is calculated by using the Poynting Vector along with the
integration line as seen in Figure~\ref{numericalsetup}. Meanwhile,
Mie theory is used to calculate the true value of total SCS for each
wavelength. The number of iterations for the main FDTD loop is set
large enough in order to give enough time for the electric and
magnetic waves to leave the simulated area for stable SCS
calculation.

Firstly, we investigate the accuracy and the convergence order of
BC-EP and compare it with other similar methods. We test all the
mthods for a circular cylinder with radius $r=400\textrm{nm}$ and
permittivity $\epsilon = 3$. The background media ia always set as
air with $\epsilon=1$. The mesh size is originally set to $\Delta
x=\Delta y = 10\textrm{nm}$. Each method has been run to calculate
the total SCS with different mesh sizes. As seen in
Figure~\ref{simulation3}, the relative errors between the numerical
values from all methods and the true solution are quite small. But
we can still see that BC-EP has the smallest error while V-EP has
the largest error. Meanwhile, for each mesh size  $\Delta x$, the
average relative error of SCS is calculated among all the
wavelengths for each method. In order to observe the error
convergence order clearly, the average relative errors versus the
mesh sizes are converted into a log scale where $N_{\lambda}$
denotes the log of the number of mesh points for the radius of the
cylinder, as seen in Figure~\ref{researchtest3}.

\begin{figure}[ht!]
  \centering
  \subfloat[Overall SCS]{\includegraphics[width=0.485\textwidth]{scse3}}
  \hfill
  \subfloat[Relative error]{\includegraphics[width=0.485\textwidth]{relerror3}}
  \caption{Total SCS calculation of a cylinder, $\epsilon = 3$ and radius $r = 400\textrm{nm}$. The grid size is $\Delta x=\Delta y = 2.8\textrm{nm}.$}
   \label{simulation3}
\end{figure}
\begin{figure}[ht!]
    \begin{center}
            \includegraphics[scale=0.5]{overallerror3}
    \end{center}
    \caption{The average relative error of a cylinder in log scale, $\epsilon = 3$ and radius $r = 400\textrm{nm}$. $N_{\lambda} =\log(\frac{400\textrm{nm}}{\Delta x}).$}
    \label{researchtest3}
\end{figure}

This whole process is then repeated for other larger permittivity
values $\epsilon=6, 10$. The numerical results are reported in
Figure~\ref{simulation6} and Figure~\ref{researchtest6} for for
$\epsilon=6$, and in Figure~\ref{simulation10} and
Figure~\ref{researchtest10} for $\epsilon=10$. other $\epsilon$
values have similar outcomes.

\begin{figure}[ht!]
  \centering
  \subfloat[Overall SCS]{\includegraphics[width=0.485\textwidth]{scse6}}
  \hfill
  \subfloat[Relative error]{\includegraphics[width=0.485\textwidth]{relerror6}}
  \caption{Total SCS calculation of a cylinder, $\epsilon = 6$ and radius $r = 400\textrm{nm}$. The grid size is $\Delta x=\Delta y = 2.8\textrm{nm}.$}
   \label{simulation6}
\end{figure}
\begin{figure}[ht!]
    \begin{center}
            \includegraphics[scale=0.5]{overallerror6}
    \end{center}
    \caption{The average relative error of a cylinder in log scale, $\epsilon = 6$ and radius $r = 400\textrm{nm}$. $N_{\lambda} =\log(\frac{400\textrm{nm}}{ \Delta x}).$}
    \label{researchtest6}
\end{figure}

\begin{figure}[ht!]
  \centering
  \subfloat[Overall SCS]{\includegraphics[width=0.485\textwidth]{scse10}}
  \hfill
  \subfloat[Relative error]{\includegraphics[width=0.485\textwidth]{relerror10}}
  \caption{Total SCS calculation of a cylinder, $\epsilon = 10$ and radius $r = 400\textrm{nm}$. The grid size is $\Delta x=\Delta y = 2.8\textrm{nm}.$}
   \label{simulation10}
\end{figure}
\begin{figure}[ht!]
    \begin{center}
            \includegraphics[scale=0.5]{overallerror10}
    \end{center}
    \caption{The average relative error of a cylinder in log scale, $\epsilon = 10$ and radius $r = 400\textrm{nm}$. $N_{\lambda} =\log(\frac{400\textrm{nm}}{ \Delta x}).$}
    \label{researchtest10}
\end{figure}

It can be seen that BC-EP gives the least amount of error for most
cases, followed by CP-EP. Staircasing results are quite erratic and
its error doesn't go down in a uniform manner due to random nature
of some mesh sizes conforming with the cylinder better than others.
V-EP has the worst error but gives a uniform error decrease with
smaller mesh sizes, and thus gives more consistent results over
staircasing.

To further investigate the order of convergence of all the methods,
we estimate the convergence order by computing the slopes from the
log figures for all the given permittivity. Table~\ref{convergence}
shows BC-EP converges significantly faster than all the other
methods. It approaches second-order in most cases. CE-EP is the
second in terms of convergence order and outperforms both
staircasing and V-EP as demonstrated in \cite{Mohammadi2005}
already.

\begin{table}[ht!]
    \centering
    \begin{tabular}{lllll} \hline
        \multicolumn{1}{c}{Permittivity} & \multicolumn{1}{c}{Staircasing} &
        \multicolumn{1}{c}{V-EP} & \multicolumn{1}{c}{CP-EP}  & \multicolumn{1}{c}{BC-EP}\\ \hline \hline
        $\epsilon = 3$ & 1.6621 & 1.0860 & 1.6977 & 2.4386 \\

       $\epsilon = 4$ & 1.0700 & 1.0850 & 1.4560 & 1.9251 \\

        $\epsilon = 5$ & 1.2160 & 1.1294 & 1.5267 & 1.6974 \\

        $\epsilon = 6$ & 1.3723 & 1.1476 & 1.3607 & 1.5247 \\

         $\epsilon = 7$ & 1.4322 & 1.1834 & 1.3800 & 1.6014 \\

        $\epsilon = 8$ & 1.4283 & 1.1860 & 1.3564 & 1.6423 \\

        $\epsilon = 9$ & 1.4867 & 1.2262 & 1.4333 & 1.6837 \\

        $\epsilon = 10$ & 1.6172 & 1.3183 & 1.5323 & 1.8006 \\ \hline
    \end{tabular}
    \caption{Order of convergence for each FDTD algorithm and for each given permittivity in a cylinder of $r=400\textrm{nm}$.}
    \label{convergence}
\end{table}

To complete the final evaluation of the performance for all the
methods, the whole test suite is then repeated for another circular
cylinder with smaller radius $r=150\textrm{nm}$. All the results are
very similar and for brevity, we only report the case of
$\epsilon=6$, as seen in Figure~\ref{2simulation6} and
Figure~\ref{2researchtest6}. It can be seen that BC-EP has a higher
convergence order than others. To be more specific, BC-EP, CP-EP,
Staircasing and V-EP have the order of convergence $1.4005$,
$0.92869$, $0.92802$ and $1.0059$, respectively, based on
Figure~\ref{2researchtest6}.

\begin{figure}[ht!]
  \centering
  \subfloat[Overall SCS]{\includegraphics[width=0.485\textwidth]{e6r1501}}
  \hfill
  \subfloat[Relative error]{\includegraphics[width=0.485\textwidth]{e6r1502}}
  \caption{Total SCS calculation of a cylinder, $\epsilon = 6$ and radius $r = 150\textrm{nm}$. The grid size is $\Delta x=\Delta y = 2.8\textrm{nm}.$}
   \label{2simulation6}
\end{figure}

\begin{figure}[ht!]
    \begin{center}
            \includegraphics[scale=0.6]{e6r1503}
    \end{center}
    \caption{The average relative error of a cylinder in log scale, $\epsilon = 6$ and radius $r = 150\textrm{nm}$. $N_{\lambda} =\log(\frac{150\textrm{nm}}{\Delta x}).$}
    \label{2researchtest6}
\end{figure}

Therefore, BC-EP is more accurate and converges faster than the
other three popular FDTD methods.

Secondly, we investigate the stability of BC-EP since the stability
of other methods have been demonstrated in the literature. The
algorithm has been run for $200,000$ iterations under FDTD main loop
which is long enough for the electric and magnetic waves completely
leave the region as seen in Figure~\ref{stable}. This process has
been repeated for $\epsilon = 3, 10, 30$ and for mesh size $\Delta
x=10\textrm{nm}, 7.3\textrm{nm}, 4.6\textrm{nm},$ and
$2.8\textrm{nm}$. At every $2,000$ iterations, a calculation of the
SCS for all $601$ wavelengths is conducted. Then the average
relative error between the calculated SCS at that iteration and the
exact SCS over those wavelengths is calculated. We observe that the
result are very similar for different $\epsilon$ and different
$\Delta x$ in the sense that the relative errors drop down to zero
very quickly. This shows that there is no electromagnetic power left
in the region and the fields return back to zero as excepted,
meaning no late time instability even at very high contrast.

\begin{figure}[ht!]
  \centering
  \subfloat[$\epsilon = 3$]{\includegraphics[width=0.3\textwidth]{stabletest3}}
  \hfill
  \subfloat[$\epsilon = 10$]{\includegraphics[width=0.3\textwidth]{stabletest10}}
  \hfill
  \subfloat[$\epsilon = 30$]{\includegraphics[width=0.3\textwidth]{stabletest30}}
  \caption{Relative error in log format for the SCS calculation at a given iteration with $r = 400\textrm{nm}$. $\Delta x=\Delta y=4.6\textrm{nm}.$}
  \label{stable}
\end{figure}

\section{Conclusions and Future Research}

We have built a new BC-EP method for the challenging two-dimensional
transverse electric Maxwell's equations involving multiple
dielectric media. Based on a relation for the electric fields across
the dielectric media interface and the integral version of the
Maxwell's equations, we derived this method which represents
Maxwell's equations accurately while the previous numerical methods
skipped some terms from the boundary conditions for the purpose of
algorithm simplification. We then apply BC-EP together with three
other popular FDTD methods, BC-EP, staircase and V-EP, to solve a
numerical example of a dielectric cylinder with given analytic
solution. The numerical results clearly demonstrated that it always
achieves the highest convergence (close to second order) compared
with other methods. The stability of BC-EP has also been verified
numerically. Due to its simplicity, BC-EP can be merged into other
FDTD software packages easily.

Our next research consideration is to extend BC-EP to solve
two-dimensional Maxwell's equations involving multiple dispersive
media. Numerical performance comparison of this method with other
reported dispersive methods (\cite{Mohammadi2008, Liu2012,
Nguyen2016}) will be provided.


text:
\begin{document}

\title[Zeros of sparse polynomials]{Zeros of sparse polynomials
	over local fields of characteristic~$p$}

\author{Bjorn Poonen}

\address{Department of Mathematics, University of California, Berkeley, CA 94720-3840, USA}
\email{poonen@math.berkeley.edu}
\date{March 2, 1998}

\maketitle

\section{Statement of results}
\label{introduction}

Let $K$ be a field of characteristic~$p>0$
equipped with a valuation $v: K^\ast \rightarrow G$
taking values in an ordered abelian group $G$.
Let $\OO_K=\{\alpha \in K: v(\alpha) \ge 0\}$
and $\mm_K=\{\alpha \in K: v(\alpha) > 0\}$
be the valuation ring and maximal ideal, respectively,
and suppose that the residue field $\OO_K/\mm_K$
is finite, with $q$ elements.

\begin{theorem}
\label{maintheorem}
If $f(x) = a_0 x^{n_0} + a_1 x^{n_1} + \cdots + a_k x^{n_k}$
is a polynomial with $k+1$ nonzero coefficients $a_i \in K^\ast$,
then $f$ has at most $q^k$ distinct zeros in $K$.
\end{theorem}

This upper bound is sharp:
if $K$ is $\F_q((T))$ with the usual
discrete valuation $v:K^\ast \rightarrow \Z$,
if $V \subset K$ is an $\F_q$-subspace of dimension $k$,
and if $c \in K$ is nonzero,
then the polynomial $f(x) := c \prod_{\alpha \in V} (x-\alpha)$
has the form $a_0 x + a_1 x^q + \cdots + a_k x^{q^k}$
for some $a_0, a_1, \ldots, a_k \in K^\ast$.

Theorem~\ref{maintheorem} is the case $d=1$
of the following generalization, which bounds the number
of distinct zeros of {\em bounded degree}.
Let $\mu(n)$ be the M\"obius $\mu$-function.

\begin{theorem}
\label{degreetheorem}
Fix $d \ge 1$.
If $f(x) = a_0 x^{n_0} + a_1 x^{n_1} + \cdots + a_k x^{n_k}$
is a polynomial with $k+1$ nonzero coefficients $a_i \in K^\ast$,
then the number of distinct zeros of $f$ in $\Kbar$ of degree
at most $d$ over $K$ is at most $\sum_{j=1}^d \sum_{i|j} q^{ik} \mu(j/i)$.
\end{theorem}

This upper bound is sharp as well, for every $q$, $k$, and $d$.
Let $K=\F_q((T))$ and $v$ be as before.
Let $\F$ be a finite field containing $\F_{q^i}$ for $i \le d$.
Let $V \subset \F((T))$ be a $k$-dimensional $\F$-vector space
that is $\Gal(\F/\F_q)$-stable (or equivalently, has an $\F$-basis
of elements of $K$).
Then equality is attained in Theorem~\ref{degreetheorem}
for $f(x) := c \prod_{\alpha \in V} (x-\alpha)$
for any $c \in K^\ast$.
(The inner sum in Theorem~\ref{degreetheorem} performs
the inclusion-exclusion to count zeros of {\em exact} degree $j$.)

We make no claim that these are the only polynomials that attain
equality; in fact there are many others.
For example, if $K$, $V$, and $f$ are as in the previous paragraph, and if
the $\F$-basis of $V$ consists of elements of $K$ of distinct valuation,
with all these valuations divisible by a single integer $e \ge 1$,
then $f(x^e)$ also attains equality,
as a short argument involving Hensel's lemma shows.
Other examples can be constructed using the observation
that if $f(x) \in K[x]$ has $N$ zeros in a given field extension $L$ of $K$,
one of which is 0,
then the same holds for $x^m f(1/x)$ when $m > \deg f$.

\begin{rem}
H.~W.~Lenstra, Jr.~\cite{lenstra}
proves related facts for finite extensions $L$ of $\Q_p$,
using very different methods.
One of his results is that for any such $L$ and any positive integer $k$,
there exists a positive integer $B=B(k,L)$ with the following property:
if $f \in L[x]$ is a nonzero polynomial with at most $k+1$ nonzero
terms and $f(0) \not= 0$, then $f$ has at most $B$ zeros in $L$,
{\em counted with multiplicities}.
His bound $B(k,L)$ is explicit, but almost certainly not sharp.
Finding a sharp bound seems difficult in general, although Lenstra
does this for the case $k=2$ and $L=\Q_2$ (the bound then is 6).
He also applies his local result to bound uniformly
the number of factors of given degree over number fields.
In~\cite{lenstra2} he shows that if $f$ is represented sparsely,
then these factors can be found in polynomial time.
\end{rem}

\begin{rem}
We cannot count multiplicities in either of {\em our} theorems and hope to
obtain a bound depending only on $k$ and $K$
(and $d$, for Theorem~\ref{degreetheorem}),
because of examples like $f(x)=(1+x)^{q^m}$ with $m \rightarrow \infty$.
Requiring that $f$ not be a $p$-th power would not eliminate the problem,
because one could also take $f(x)=(1+x)^{q^m+1}$.
\end{rem}

\section{Proof of Theorem~\ref{maintheorem}}
\label{proof1}

By a {\em disk} in a valued field $K$,
we mean either
an ``open disk'' $D(x_0,g):=\{x \in K: v(x-x_0)>g\}$,
or a ``closed disk'' $\Dbar(x_0,g):=\{x \in K: v(x-x_0) \ge g\}$
where $x_0 \in K$ and $g \in G$.

Let $\sigma_1$, $\sigma_2$, \dots, $\sigma_t$ be the
non-vertical segments of the Newton polygon of $f$.
Let $-g_j \in G \tensor \Q$ be the slope of $\sigma_j$.
If $e_1,e_2,\ldots,e_r$ are the exponents
of the monomials in $f$ corresponding to points on a given $\sigma_j$,
define $N_j$ as the largest integer for which
the images of $(1+x)^{e_1}$, $(1+x)^{e_2}$, \dots, $(1+x)^{e_r}$
in $\F_p[x]/(x^{N_j})$ are linearly {\em dependent} over $\F_p$.
We say that the $\sigma_j$ are in a {\em proper order}
if $N_1 \ge N_2 \ge \cdots \ge N_t$.
This particular ordering is crucial to the proof,
but it is hard to motivate its definition.
It was discovered by analyzing proofs
of many special cases of Theorem~\ref{maintheorem}.
For instance, if the Newton polygon of $f$ has $k$ non-vertical segments
(each associated with exactly two exponents),
then the segments are being ordered according to the $p$-adic absolute values
of their horizontal lengths.

\begin{lemma}
\label{distances}
Let $L$ be a field of characteristic~$p>0$
with a valuation $v: L^\ast \rightarrow G$.
Suppose $f(x) = a_0 x^{n_0} + a_1 x^{n_1} + \cdots + a_k x^{n_k} \in L[x]$
with each $a_i$ nonzero.
List the segments of the Newton polygon of $f$ in a proper order as above.
Fix $u$ and let $-g_u \in G \tensor \Q$ be the slope of
the $u$-th segment $\sigma_u$.
Suppose $r \in L$ is {\em not} a zero of $f$, and $v(r)=g_u$.
Let $S$ be the set of zeros of $f$ in $L$ lying inside $D(r,g_u)$.
Then $\# \{v(\alpha-r): \alpha \in S \} \le k+1-u$.
\end{lemma}

\begin{proof}
Replacing $f(x)$ by $cf(rx)$ for suitable $c \in L^\ast$,
we may reduce to the case in which $r=1$, $g_u=0$, $f \in \OO_L[x]$,
and $f \bmod \mm_L$ is nonzero.
Write $f(1+x)=\sum_{j=0}^{n_k} b_j x^j$,
and let $M$ be the smallest integer for which $b_M$ is nonzero modulo $\mm_L$.
By definition of $N_u$, $f(1+x) \not\equiv 0 \bmod (\mm_L,x^{N_u+1})$.
Hence $M \le N_u$.

For each $i \le u$, we have $N_u \le N_i$, so $M \le N_i$,
and there is some $\F_p$-linear relation in $\F_p[x]/(x^M)$
between the $(1+x)^e$ for the exponents $e$ associated to $\sigma_i$.
The subspace of $\F_p[x]/(x^M)$ spanned by the $(1+x)^e$,
where $e$ ranges over {\em all} the exponents in $f$,
then has dimension at most $(k+1)-u$,
since the $u$ relations above are independent,
the largest $e$ involved in each relation being distinct from
the others.
It follows that the $\F_p$-subspace of $K$
spanned by $b_0$, $b_1$, \dots, $b_{M-1}$
is at most $(k+1-u)$-dimensional.
Then $\#\{v(b_i): 0 \le i < M \text{ and } b_i \not=0 \} \le k+1-u$,
because nonzero elements of distinct valuations
are automatically $\F_p$-independent.
The left endpoints of the
negative slope segments of the Newton polygon of $f(1+x)$
correspond to $b_i$ of distinct valuations for $i<M$,
so there are at most $k+1-u$ such segments.
Hence at most $k+1-u$ positive elements of $G$ can
be valuations of zeros of $f(1+x)$,
which is what we needed to prove.
\end{proof}

\begin{rem}
Note that there is no assumption on the residue field
in Lemma~\ref{distances}; $L$ could even be algebraically closed.
\end{rem}

Let $S$ be any finite subset of a field $L$ with valuation $v$.
We associate a tree $\TT$ to $S$ as follows.
(See~\cite{stanley} for arboreal terminology.)
Let $\TT$ be the Hasse diagram of the finite poset
(ordered by inclusion) of nonempty sets
of the form $S \cap D$ where $D$ is a disk.
Clearly $\TT$ is a tree, whose leaves are the singleton subsets of $S$.
We would obtain the same tree if we required the disks $D$ to be
open (resp.\ closed), since $S$ is finite.

Suppose $r$ and $S$ are as in Lemma~\ref{distances}.
Let $T_0 > T_1 > \cdots > T_\ell$ be the longest chain in $\TT$.
Then $T_\ell$ is a leaf, and $\# T_\ell=1$.
Choose $r_0 \in D(r,g_u) \setminus S$ closer to the element of $T_\ell$
than to any other element of $S$.
For various $g>g_u$, the set $S \cap D(r_0,g)$ can equal
$T_0$, $T_1$, \dots, $T_\ell$, or $\emptyset$.
Hence
	$$\# \{v(\alpha-r_0): \alpha \in S\} \ge \ell+1.$$
On the other hand, Lemma~\ref{distances} applied to $r_0$ yields
	$$\# \{v(\alpha-r_0): \alpha \in S\} \le k+1-u.$$
Combining these, we have that {\em the length $\ell=\ell(\TT)$
of the tree satsifies $\ell \le k-u$.}

Suppose $S_0 \in \TT$ is not a leaf (i.e. $\#S_0>1$),
and let $g=\min\{v(s-t):s,t \in S_0\}$,
so that for any $s \in S_0$,
$\Dbar(s,g)$ is the smallest disk containing $S_0$.
Then the children of $S_0$ in the tree
are nonempty sets of the form $S \cap D(x_0,g)$
for some $x_0 \in \Dbar(s,g)$.
In particular the number of children is at most the
size of the residue field of $L$.

\begin{proof}[Proof of Theorem~1]
Let notation be as in Lemma~\ref{distances}, but take $L=K$.
By the theory of Newton polygons, each nonzero zero of $f$
has valuation equal to $g_u$ for some $u$.
Let us now fix $u$ and let $Z_u$
be the number of zeros in $K$ of valuation $g_u$.
We may assume $g_u \in G$, since otherwise $Z_u=0$.
Then $\{x \in K: v(x)=g_u\}$ is the union of $q-1$ open disks $D_j$
of the form $D(x_j,g_u)$.
As above, the tree corresponding to the set of zeros in $D_j$
has length at most $k-u$, and each vertex has at most $q$ children.
Hence the tree has at most $q^{k-u}$ leaves,
and $Z_u \le (q-1)q^{k-u}$.
Allowing for the possibility that $0$ also is a zero of $f$,
we find that the number of zeros of $f$ in $K$ is at most
	$$1+\sum_{u=1}^t Z_u \le 1 + \sum_{u=1}^t (q-1)q^{k-u} \le 1 + \sum_{u=1}^k (q-1)q^{k-u} = q^k.$$
\end{proof}

\section{Valuation theory}

Before proving Theorem~\ref{degreetheorem},
we will need to recall some facts from valuation theory.
We write $(K,v)$ for a field $K$ with a valuation $v$.
We say that $(L,w)$ is an {\em extension} of $(K,v)$
if $K \subseteq L$ and $w|_K=v$.
In this case,
when we say that $L$ has the same value group (resp.\ residue field)
as $K$,
we mean that the inclusion of value groups (resp.\ residue fields)
induced from the inclusion of $(K,v)$ in $(L,w)$
is an isomorphism.
Recall that any valuation on a field $K$ admits at least
one extension to any field containing $K$.
An abelian group $G$ is {\em divisible} if for all $g \in G$ and $n \ge 1$,
the equation $nx=g$ has a solution $x$ in $G$.

\begin{prop}
\label{makedivisible}
Any valued field can be embedded in another valued field
having the same residue field, but divisible value group.
\end{prop}

\begin{proof}
Let $v: K^\ast \rightarrow G$ be the original valuation.
If $G$ is not already divisible, then there exists $g \in G$
and a prime number $n$ such that $n x=g$ has no solution in $G$.
Pick $\alpha \in K^\ast$ with $v(\alpha)=g$, and extend
$v$ to a valuation on $L=K(\alpha^{1/n})$.
Let $e$ and $f$ denote the ramification index and residue class degree
for $L/K$.
Then $e=n$, and the inequality $ef \le n$

(Lemma~18 in Chapter~1 of~\cite{schilling})
forces $f=1$.
An easy Zorn's lemma argument now shows that $v$ extends
to a valuation $v:M^\ast \rightarrow G \tensor \Q$ where
$M$ is an extension with the same residue field as $K$,
but with divisible value group.
\end{proof}
Recall that $(L,w)$ is called an {\em immediate extension} of $(K,v)$ if
\begin{enumerate}
	\item $(L,w)$ is an extension of $(K,v)$;
	\item $(L,w)$ has the same value group as $(K,v)$; and
	\item $(L,w)$ has the same residue field as $(K,v)$.
\end{enumerate}
Also recall that $(K,v)$ is called {\em maximally complete}
if it has no nontrivial immediate extensions.

\begin{prop}
\label{krullresult}
Every valued field has a maximally complete immediate extension.
\end{prop}

\begin{proof}
This is an old result of Krull:
see Theorem~5 of Chapter~2 in~\cite{schilling}.
\end{proof}

\begin{prop}
\label{embedFq}
Suppose that $(K,v)$ is maximally complete of characteristic $p>0$,
and that $\F_q$ is contained in the residue field.
Then $\F_q$ can be embedded in $K$.
\end{prop}

\begin{proof}
Apply a suitable version of Hensel's lemma
(combine Theorems~6 and~7 of Chapter~2 of~\cite{schilling})
to the factorization of $x^q-x$ over $\F_q$.
\end{proof}

\begin{prop}
\label{finiteextensions}
Suppose that $(K,v)$ is maximally complete of characteristic $p>0$,
with divisible value group $G$ and with residue field $\F_q$.
If $L \subset \Kbar$ is a finite extension of $K$ of degree $n$, then
$L$ is the compositum of $\F_{q^n}$ and $K$ in $\Kbar$.
\end{prop}

\begin{proof}
Extend $v$ to $L$.
Theorem~11 in Chapter~2 of~\cite{schilling} shows that
$L$ is maximally complete,
and that $ef=n$ holds for $L/K$.
Since $G$ is divisible, there are no ordered abelian groups $G'$
with $1<(G':G)<\infty$.
Hence $e=1$, $f=n$, and the residue field of $L$ is $\F_{q^n}$.
Proposition~\ref{embedFq} implies that
the subfield $\F_{q^n}$ of $\Kbar$ is contained in $L$.
But the compositum of the linearly disjoint
fields $\F_{q^n}$ and $K$ in $\Kbar$ is already
$n$-dimensional over $K$, so the compositum must equal $L$.
\end{proof}

\begin{rem}
	Lenstra notes that if one is interested in proving
	Theorem~\ref{degreetheorem}
	only for polynomials over $K_0=\F_q((T))$, then one can circumvent
	the theory of maximally complete fields by choosing
	$\sigma \in \Gal(\Kbar_0/K_0)$ that acts as $x \mapsto x^q$
	on $\Fbar_q$, and by taking $K$ to be the fixed field of $\sigma$.
	This $K$ contains $K_0$, still has residue field $\F_q$, and satisfies
	the conclusion of Proposition~\ref{finiteextensions}.
\end{rem}

\section{Proof of Theorem~\ref{degreetheorem}}

In proving Theorem~\ref{degreetheorem},
we may first apply Propositions~\ref{makedivisible} and~\ref{krullresult}
to assume that the value group $G$ is divisible and that
$(K,v)$ is maximally complete (still with residue field $\F_q$).
Let $\F=\F_{q^{d!}} \subset \Kbar$.
Proposition~\ref{finiteextensions}
shows that all elements of $\Kbar$ of degree at most $d$ over $K$
lie inside the compositum $L:=\F \cdot K$ of fields in $\Kbar$.
Extend $v$ to $L$.

For each $g \in G$, choose $\beta_g \in K$ with $v(\beta_g)=g$.
Now suppose $\Dbar:=\Dbar(x_0,g)$ is a closed ball in $L$.
Let $I$ be the subgroup of $\Gal(L/K) \isom \Gal(\F/\F_q)$
that maps $\Dbar$ into $\Dbar$.
Division by $\beta_g$ induces an isomorphism of $I$-modules
$\Dbar(0,g)/D(0,g) \isom \F$,
so the cohomology group $H^1(I,\Dbar(0,g)/D(0,g))$ is trivial.
The long exact sequence associated with the
exact sequence
	$$0 \rightarrow \frac{\Dbar(0,g)}{D(0,g)}
		\rightarrow \frac{L}{D(0,g)}
		\rightarrow \frac{L}{\Dbar(0,g)}
		\rightarrow 0$$
of $I$-modules shows that $\Dbar$
contains an open disk $D(x_1,g)$ mapped to itself by $I$.
We then have a bijection of $I$-sets
$\phi_{\Dbar}: \Dbar/D(0,g) \rightarrow \F$
that maps the coset $y+D(0,g)$ to the residue class of $(y-x_1)/\beta_g$.
We assume that the elements $\beta_g$ and the maps $\phi_{\Dbar}$
are fixed once and for all.

Now let $g_u$, $r$, and $S$ be as in Lemma~\ref{distances},
and let $\TT$ be the tree associated to $S$ as in Section~\ref{proof1},
so that $\ell(\TT) \le k-u$.
We now describe a labelling of the vertices of $\TT$ by elements of $\F$.
Recall that if $S_0 \in \TT$ is not a leaf,
and if $\Dbar=\Dbar(s,g)$ is the smallest disk containing $S_0$,
then the children of $S_0$ are nonempty sets of
the form $S \cap D(x_0,g)$ for some $x_0 \in \Dbar$.
Label each child by $\phi_{\Dbar}(D(x_0,g))$.
Note that the children of $S_0$ are labelled with {\em distinct}
elements of $\F$.
Finally, label the root of $\TT$
with the residue of $r/\beta_{g_u}$ in $\F^\ast$.

Let $R$ be the set of all roots of $f$ in $L$,
and let $\F[X]_{<k}$ denote the set of polynomials
of the form $a_0 + a_1 X + \cdots + a_{k-1} X^{k-1}$ with $a_i \in \F$.
We now define a map $\Phi: R \rightarrow \F[X]_{<k}$.
First, if $0 \in R$, define $\Phi(0)=0 \in \F[X]_{<k}$.
If $z \in R$ is nonzero, then $v(z)=g_u$ for some $u$.
Let $T_0 > T_1 > \cdots > T_n$ be the maximal chain ending at $T_n=\{z\}$
in the tree $\TT$ associated to $S:=R \cap D(z,g_u)$.
Define $\Phi(z)=X^{u-1} \sum_{i=0}^n \text{label}(T_i) X^i.$
Since $n \le \ell(\TT) \le k-u$, we have $\Phi(z) \in \F[X]_{<k}$.

\begin{lemma}
\label{Philemma}
$\nichts$
\begin{enumerate}
\item The map $\Phi: R \rightarrow \F[X]_{<k}$ is injective.
\item If $z \in R$ is of degree $j$ over $K$, then $\Phi(z) \in \F_{q^j}[X]$.
\end{enumerate}
\end{lemma}

\begin{proof}
To prove injectivity, we describe how to reconstruct $z$ from $\Phi(z)$.
If $\Phi(z)=0$, then $z$ must be $0$.
Otherwise its lowest degree monomial involves $X^{u-1}$ where $v(z)=g_u$.
Hence, assuming from now on that $z \not=0$,
we can reconstruct $v(z)$ from $\Phi(z)$.
Next, the {\em coefficient} of $X^{u-1}$ determines which
(nontrivial) coset of $D(0,g_u)$ in $\Dbar(0,g_u)$ $z$ belongs to.
The other coefficients uniquely determine a path ending at
the leaf $\{z\}$ in the tree associated to this coset.
Thus $\Phi(z)$ determines $z$.

For the second part, it suffices to show that
if $H$ is the subgroup of $\Gal(L/K)$ fixing $z \in R$,
then $H$ (or equivalently the isomorphic subgroup of $\Gal(\F/\F_q)$)
fixes the coefficients of $\Phi(z)$ also.
We may assume $z \not=0$.
Let $g_u=v(z)$,
and let $T_0 > T_1> \cdots > T_n=\{z\}$ be the maximal chain
in the tree $\TT$ associated to the coset $z+D(0,g_u)$ in which $z$ lies.
Since $H$ preserves the coset $z+D(0,g_u)$,
$H$ fixes the label of $T_0$.
Now suppose $1 \le i \le n$.
The smallest disk containing $T_{i-1}$ is of the form $\Dbar:=\Dbar(z,g)$
for some $g > g_u$, so $H$ is contained
in the subgroup $I \subseteq \Gal(L/K)$ preserving this disk.
The label of the child $T_i$ is $\phi_{\Dbar}(z+D(0,g))$,
and $\phi_{\Dbar}$ respects the action of $H \subseteq I$,
so $H$ fixes this label.
This holds for all $i$, so $H$ fixes all coefficients of $\Phi(z)$.
\end{proof}

Lemma~\ref{Philemma} shows that the number of zeros of $f$ in $\Kbar$
of degree at most $d$ is less than or equal to the number of polynomials
in $\F[X]_{<k}$ that are defined over $\F_{q^j}$ for some $j \le k$.
The number of such polynomials defined over $\F_{q^j}$
but no subfield is  $\sum_{i|j} q^{ik} \mu(j/i)$, by M\"obius inversion.
Theorem~2 follows upon summing over $j$.

\section*{Acknowledgements}
I thank Hendrik Lenstra for bringing the problems to my attention,
and for suggesting that $q^k$ might be the correct bound
for Theorem~\ref{maintheorem}.

config: arxiv


TITLE: Solving a general ODE
QUESTION [1 upvotes]: $ \newcommand{\der}[1]{\frac{d^#1 u}{dx^#1}} $
I'm trying to solve
$$\der n =\sum^{n-1}_{k=0}\alpha_k\der k +e^{-x}.$$
It can be transformed to an set of equations:
$$\begin{cases}
u_0'=u_1\\
u_1'=u_2\\
\vdots\\
u_{n-1}'= \sum^{n-1}_{k=0}\alpha_ku_k +e^{-x}
\end{cases}$$
And in matrix form:
$$\begin{pmatrix}u_0\\u_1\\u_2\\\vdots\\u_{n-1}\end{pmatrix}'=\begin{pmatrix}
0&1&0&\dots&0\\
0&0&1&\dots & 0\\
0&0&0&\ddots&0\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
\alpha_0&\alpha_1&\alpha_1&\dots&\alpha_{n-1}\\
\end{pmatrix}\begin{pmatrix}u_0\\u_1\\u_2\\\vdots\\u_{n-1}\end{pmatrix}+\begin{pmatrix}0\\0\\0\\\vdots\\e^{-x}\end{pmatrix}$$
Now, I can start by solving the homogenous form:
$$\begin{pmatrix}u_0\\u_1\\u_2\\\vdots\\u_{n-1}\end{pmatrix}'=\begin{pmatrix}
0&1&0&\dots&0\\
0&0&1&\dots & 0\\
0&0&0&\ddots&0\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
\alpha_0&\alpha_1&\alpha_1&\dots&\alpha_{n-1}\\
\end{pmatrix}\begin{pmatrix}u_0\\u_1\\u_2\\\vdots\\u_{n-1}\end{pmatrix}$$
I've got stuck here and haven't managed to diagnalize the matrix to exponentiate it.

REPLY [1 votes]: The matrix in the homogeneous differential equation given in your question is also referred to as companion matrix.
Some interesting properties (eigenvalues, eigenvectors, etc.) are summarized in 'The Companion Matrix and Its Properties. Brand, L., The American Mathematical Monthly, 6, 71, p. 629--634, 1964'.
The eigenvalues of this matrix correspond to zeros of the polynomial
$$
f(\lambda)=\lambda^n - \alpha_{n-1} \lambda^{n-1} - \ldots - \alpha_1\lambda-\alpha_0.
$$
It is not fully clear from your question if you try to find an analytic solution or if you want to compute a numerical solution. I would not expect a analytic solution to be available in general. Considering numerical solutions, hump effects can occur for the exponential of the companion matrix which can result in numerical difficulties. E.g. for $f(\lambda)=(\lambda+w)^n$ with $w\gg 0$ the companion matrix with the respective choice of coefficients $\alpha_0,\ldots,\alpha_{n-1}$ does have an eigenvalue $-w\ll 0$ with multiplicity $n$ but the $\mu$-norm of this matrix scales with $w^n$ (hump effects for short time steps).


TITLE: Boiling potatoes
QUESTION [1 upvotes]: If the temperature of boiling water is about 100 ºC (depending on air pressure), why are boiled potatoes on a high stove flame cooked faster?
From the English Wikipedia:

Acrylamide was discovered in foods in April 2002 by Eritrean scientist
  Eden Tareke in Sweden when she found the chemical in starchy foods,
  such as potato chips (potato crisps), French fries (chips), and bread
  that had been heated higher than 120 °C (248 °F) (production of
  acrylamide in the heating process was shown to be
  temperature-dependent). It was not found in food that had been boiled
  or in foods that were not heated.

If I boil potatoes on the highest stove flame, will the potatoes not reach the temperature of 120 ºC?

REPLY [1 votes]: Whether your stove is on 'high' or 'low' doesn't change the temperature of the water, i.e. the boiling point remains $100\ \mathrm{C}$. So the potatoes can NEVER reach $120\ \mathrm{C}$.
But objects being heated by hot water follow Newton's Law of Cooling (and Heating), which describes the heat transfer per unit of time, $\frac{dQ}{dt}$:
$$\frac{dQ}{dt}=hA[T_{BP}-T_{pot}(t)]$$
We can derive from this, for potatoes plunged into boiling hot water at $t=0$, with $T_0$ their ambient (and initial) temperature and $T_e$ their end temperature, their temperature evolution in time $t$:
$$T_e=T_{BP}-(T_{BP}-T_0)e^{-\frac{hA}{mc_p}t}$$
And the heating time:
$$t=\frac{mc_p}{hA}\times \ln\Big({\frac{T_{BP}-T_0}{T_{BP}-T_e}}\Big)$$
($m$ is mass of potatoes, $c_p$ their specific heat capacity and $A$ their surface area)
The factor $h$ is the heat transfer coefficient (in $\mathrm{Wm^{-2}K}$). It is obvious that when $h$ goes up, $t$ goes down.
In fluids, like water, $h$ is known to be somewhat dependent on the Reynolds Number, which is higher in turbulent conditions. More turbulent conditions favour higher heat transfer rates.
For that reason there may be a slight increase in heat transfer (and thus a slight reduction in cooking time) when comparing a vigorous boil (high heat) to a gentle simmer (low heat).


TITLE: Calculating Basis Functions for DFTs (64 Samples)
QUESTION [1 upvotes]: I am attempting to graph some 64 sample'd basis functions in MatLab, and getting inconsistent results -- which is to say, I'm getting results that are still sinusoidal, but don't have the frequency they ought.
Here's a graph of what is supposed to be my c8 basis function:

Unfortunately, it only has 7 peaks, which indicates that I seem to have botched the frequency somehow. I'm assuming my problem lies somewhere within how I'm trying to graph in matlab, and not an error in the function itself.
Here's my code:
n = linspace(0, 2*pi*8, 64)
x = cos(2*pi*8*n/64)
plot(n,x)

I'm inclined to believe x has the correct formula, but I'm at a loss as to how else to formulate an 'n' to graph it with. Why am I getting a result with the incorrect frequency?

REPLY [1 votes]: You're plotting the function $\cos n\pi/4$, which has period $8$, and thus $8$ full periods from $0$ to $64$, but you're only substituting values from $0$ to $16\pi$. Since $16\pi\approx50$, you're missing a bit less than two of the periods. From what it seems you're trying to do, you should be plotting the function from $0$ to $64$, i. e. replace 2*pi*8 by 64 in the first line.