## 1 FDTD discretization

In finite-difference time-domain (FDTD), the discretization of the spatial domain is driven by the mesh definition. The mesh or grid is comprised of grid cells, with the size of these cells establishing the simulation\('\)s spatial resolution. The conventional grid cell dimensions in FDTD varies from \(\lambda/10\) to \(\lambda/20\), where \(\lambda\) is the smallest wavelength in the simulation source. It is possible for a coarse mesh using the Yee grid to introduce a nonphysical numerical dispersion, which results in the simulated wave travelling with a slower phase velocity. While it is important to note that structures often have fine features or curved surfaces that may require a finer mesh, this article will focus on the simulation resolution relative to the source wavelength.

### 1.2 Numerical dispersion

Solving Maxwell\('\)s equations using the Yee grid [1] results in a discrepancy between the analytical and simulated solutions. In this regard, the phase velocity of the numerical wave deviates slightly from its true value depending on the frequency, grid resolution, and the direction of the propagating wave. To compare the analytic and simulated phase velocities, one must first compute the wavevector (\(????\)), which is related to the phase velocity as \[\begin{equation} v_{p}=\frac{\omega}{k} \end{equation}\] where \(\omega\) is the angular frequency.

#### 1.2.1 Analytical Dispersion Equation

For the comparison of the analytical and simulation results this article will consider a propagating monochromatic planewave in freespace. In the 2D domain, this planewave (\(E_{y}=E_{y_ {0}}e^{j\left(\omega t-k_{x}x-k_{z}z\right)}\)) along with Maxwell\('\)s equations yields the analytic dispersion equation[1-2] \[\begin{equation} \label{analytic}k_{x}^{2}+k_{z}^{2}=\left(\frac{\omega}{c}\right)^{2} \end{equation}\] where \(c\) represent the speed of light and the analytical wavevector components along x and z are defined as \(k_{x}\) and \(k_{z}\), respectively.

#### 1.2.2 Numerical Dispersion Equation

For the simulated case, Maxwell\('\)s equations need to be expressed in their discretized form. The source free expressions for the TE polarization (i.e. Hx, Hz, Ey) are written as [1-2] \[\begin{eqnarray} \label{continuous1}\frac{\partial E_{y}}{\partial t}&=&\frac{1}{\varepsilon_{0}}\left (\frac{\partial H_{x}}{\partial z}-\frac{\partial H_{z}}{\partial x}\right)\\ \label{continuous2}\frac{\partial H_{x}}{\partial t}&=&\frac{1}{\mu_{0}}\left(\frac {\partial E_{y}}{\partial z}\right)\\ \label{continuous3}\frac{\partial H_{z}}{\partial t}&=&\frac{1}{\mu_{0}}\left(\frac {-\partial E_{y}}{\partial x}\right) \end{eqnarray}\] Discretizing the partial derivatives in the above equations is achieved through the use of the Yee cell [1-2] and yields \[\begin{eqnarray} \label{discretized1} E_{y}^{n+1}(i,k) & = & E_{y}^{n}(i,k)+\frac{\Delta t}{\varepsilon_{0}\,\Delta z}\left[H_{x}^{n+0.5}(i,k)-H_{x}^{n+0.5}(i,k-1)\right]\nonumber \\ & & -\frac{\Delta t}{\varepsilon_{0}\,\Delta x}\left[H_{z}^{n+0.5}(i,k)-H_{z}^ {n+0.5}(i-1,k)\right]\\ \label{discretized2}H_{x}^{n+0.5}(i,k)&=&H_{x}^{n-0.5}(i,k)+\frac{\Delta t}{\mu_{0}\,\Delta z}\left[E_ {y}^{n}(i,k+1)-E_{y}^{n}(i,k)\right]\\ \label{discretized3}H_{z}^{n+0.5}(i,k)&=&H_{z}^{n-0.5}(i,k)-\frac{\Delta t}{\mu_{0}\,\Delta x}\left[E_ {y}^{n}(i+1,k)-E_{y}^{n}(i,k)\right] \end{eqnarray}\] where the discretized field components are written as [2] \[\begin{eqnarray} \label{component1}E_{y}\mid_{I,K}^{n}&=&E_{y_{0}}e^{j\left(\omega n\Delta t-\tilde{k}_{x} I\Delta x-\tilde{k}_{z}K\Delta z\right)}\\ \label{component2}H_{x}\mid_{I,K}^{n}&=&H_{x_{0}}e^{j\left(\omega n\Delta t-\tilde{k}_{x} I\Delta x-\tilde{k}_{z}K\Delta z\right)}\\ \label{component3}H_{z}\mid_{I,K}^{n}&=&H_{z_{0}}e^{j\left(\omega n\Delta t-\tilde{k}_ {x}I\Delta x-\tilde{k}_{z}K\Delta z\right)} \end{eqnarray}\] where \(\tilde{k}_{x}\) and \(\tilde{k}_{z}\) are numerical wavevectors in x and z direction. The indices I and K denote the spatial indices in the x and z directions, respectively. Inserting the discretized fields into the discretized Maxwell\('\)s equations leads to \[\begin{equation} \left[\left(\frac{1}{c\Delta t}\right)\sin\left(\frac{\omega\Delta t}{2}\right)\right]^ {2}=\left[\frac{1}{\Delta x}\sin\left(\frac{\tilde{k_{x}}\Delta x}{2}\right)\right]^{2} +\left[\frac{1}{\Delta z}\sin\left(\frac{\tilde{k_{z}}\Delta z}{2}\right)\right]^{2} \end{equation}\] It can be shown that the discretized dispersion equation (12) simplifies to the analytical expression (2) for a square grid (\(\Delta x = \Delta z = \Delta\)) under the limit that \(\Delta->0\) and using the Courant stability criteria (\(\Delta t = \frac{\Delta}{\sqrt{2}c}\)).

#### 1.2.3 Numerical Dispersion Relative to Frequency and Resolution

To explore the deviation between analytical and discretized dispersion equations relative to frequency and resolution it is easier to reduce the problem to 1D (\(k_z=0\)) . This action modifies the analytical (1) and discretized (12) equations to be \[\begin{eqnarray} \frac{\omega\Delta x}{c} & = & |k_x\Delta x|\qquad\text{analytical}\\ \frac{\omega\Delta x}{c} &=& \left|\frac{2}{S}sin^{-1}\left(S\cdot sin\left(\frac {k_x\Delta x}{2}\right)\right)\right|\qquad\text{discretized} \end{eqnarray}\] where \(S=\frac{c\Delta t}{\Delta x}\). Note that the analytical expression has been normalized by \(\Delta x\) for the purpose of comparison. Plotting both of these expressions demonstrates a number of trends in how the discretized dispersion equation deviates from the analytical one, see figure 1.

The observations that can be made are that the deviation:

- proportional to the frequency.
- proportional to the size of the grid cells (inversely proportional to the number of points per wavelength).
- inversely proportional to the ratio of \(\Delta t\) to \(\Delta_x\), i.e. \(S\). Note that while S = 1 would match the analytical and discretized results this is at the edge of numerical stability and as such is not advisable.

#### 1.2.4 Numerical Dispersion Relative to Propagation Direction

To explore numerical dispersion relative to propagation direction the
analytical and discretized expressions can be rewritten relative to a
general wavevector \(k\) such that
\[\begin{eqnarray}
k_x &=& kcos\phi\\
k_z &=& ksin\phi
\end{eqnarray}\] With this identity the analytical expression (2)
becomes \[\begin{eqnarray}
\left(\frac{\omega\Delta_x}{c}\right)^2 & = & k^2(cos^2\phi +
sin^2\phi) \\
\left(\frac{\omega\Delta_x}{c}\right)^2 & = & k^2
\end{eqnarray}\] and similarly the discretized expression (12),
on a square grid, becomes \[\begin{eqnarray}
\left(\frac{\Delta}{c\Delta t}\right)^2sin^2\left(\frac{\omega\Delta
t}{2}\right) & = &
sin^2\left(\frac{\tilde{k}
\Delta cos\phi}{2}\right)+sin^2\left(\frac{\tilde{k}
\Delta sin\phi}{2}\right)
\end{eqnarray}\] Solving (19) with respect to \(\tilde{k}\) requires the use of Newton’s
method on the following expression \[\begin{equation}
\label{dispersion}f\left(\tilde{k}\right)=\sin^{2}\left(\frac{\tilde{k}\Delta\cos\varphi}{2}\right)+\sin^{2}\left(\frac{\tilde{k}\Delta\sin\varphi}{2}\right)-\left(\frac{\Delta}{c\Delta
t}\right)^{2}\sin^{2}\left(\frac{\pi c\Delta t}{\lambda_{0}}\right)=0
\end{equation}\] The algorithm [2] is given as \[\begin{equation}
\label{Newton}\tilde{k}_{n+1}=\tilde{k}_{n}-\frac{f\left(\tilde{k}\right)}{f^{'}\left(\tilde{k}\right)}=\tilde{k}_{n}-\frac{\sin^{2}\left(C_{1}\tilde{k}_{n}\right)+\sin^{2}\left(C_{2}\tilde{k}_{n}\right)-C_{3}}{C_{1}\sin\left(2C_{1}\tilde{k}_{n}\right)+C_{2}\sin\left(2C_{2}\tilde{k}_{n}\right)}
\end{equation}\] where \(f^{'}\) represents the first derivative
of the function \(f\) and

\[\begin{eqnarray}
C_{1}&=&\frac{\Delta\cos\varphi}{2}\\
C_{2}&=&\frac{\Delta\sin\varphi}{2}\\
C_{3}&=&\left(\frac{\Delta}{c\Delta
t}\right)^{2}\sin^{2}\left(\frac{\pi c\Delta t}
{\lambda_{0}}\right)
\end{eqnarray}\] The numerical phase velocity \({\tilde{v}_{p}}=\frac{\omega}{\tilde{k}}\)
can now be calculated from (1). Figure 2 shows the normalized numerical
phase velocity versus the angle of propagation, \(\varphi\), for three cases i.e. when the
number of grid cells per wavelength are 4, 8, and 16. In all plots the
time step is set as \(S=\frac{c\Delta
t}{\Delta}=\frac{1}{\sqrt{2}}\). As can be seen from the results,
the numerical phase velocity is lower than its true value, c (shown by
the dashed line), and that the deviation between the analytical and
discretized results can be minimized by increasing the number of grid
cells.