The Finite-Difference Time-Domain (FDTD) method[1,2,3] is a state-of-the-art method for solving Maxwell's equations in complex geometries. Being a direct time and space solution, it offers the user a unique insight into all types of problems in electromagnetics and photonics. In addition, FDTD can also obtain the frequency solution by exploiting Fourier transforms, thus a full range of useful quantities can be calculated, such as the complex Poynting vector and the transmission / reflection of light.
Solver Physics
This section will introduce the basic mathematical and physics formalism behind the FDTD algorithm.
FDTD solves Maxwell's curl equations in non-magnetic materials:
\begin{array}{l}{\frac{\partial \overrightarrow{D}}{\partial t}=\nabla \times \overrightarrow{H}} \\ {\overrightarrow{D}(\omega)=\varepsilon_{0} \varepsilon_{r}(\omega) \overrightarrow{E}(\omega)} \\ {\frac{\partial \overrightarrow{H}}{\partial t}=-\frac{1}{\mu_{0}} \nabla \times \overrightarrow{E}}\end{array}
where H, E, and D are the magnetic, electric, and displacement fields, respectively, while \(\epsilon_{r}(\omega) \) is the complex relative dielectric constant (\( \epsilon_{r}(\omega)=n^2 \), where n is the refractive index).
In three dimensions, Maxwell equations have six electromagnetic field components: Ex, Ey, Ez and Hx, Hy, and Hz. If we assume that the structure is infinite in the z dimension and that the fields are independent of z, specifically that
\begin{array}{l}{\varepsilon_{r}(\omega, x, y, z)=\varepsilon_{r}(\omega, x, y)} \\ {\frac{\partial \overrightarrow{E}}{\partial z}=\frac{\partial \overrightarrow{H}}{\partial z}=0}\end{array}
then Maxwell's equations split into two independent sets of equations composed of three vector quantities each which can be solved in the x-y plane only. These are termed the TE (transverse electric), and TM (transverse magnetic) equations. We can solve both sets of equations with the following components:
TE: Ex, Ey, Hz
TM: Hx, Hy, Ez
For example, in the TM case, Maxwell's equations reduce to:
\begin{array}{l}{\frac{\partial D_{z}}{\partial t}=\frac{\partial H_{y}}{\partial x}-\frac{\partial H_{x}}{\partial y}} \\ {D_{z}(\omega)=\varepsilon_{0} \varepsilon_{r}(\omega) E_{z}(\omega)} \\ {\frac{\partial H_{x}}{\partial t}=-\frac{1}{\mu_{0}} \frac{\partial E_{z}}{\partial y}} \\ {\frac{\partial H_{y}}{\partial t}=\frac{1}{\mu_{0}} \frac{\partial E_{z}}{\partial x}}\end{array}
The FDTD method solves these equations on a discrete spatial and temporal grid. Each field component is solved at a slightly different location within the grid cell (Yee cell), as shown below. By default, data collected from the FDTD solver is automatically interpolated to the origin of each grid point, so the end user does not have to deal with this issue in their analysis.
Dispersive materials with tabulated refractive index (n,k) data as a function of wavelength can be incorporated by using the multi-coefficient material models that automatically generates a material model based on the tabulated data. Alternatively, specific models such as Plasma (Drude), Debye or Lorentz can be used. See the Materials section for details. The FDTD solver supports a range of boundary conditions, such as PML, periodic, and Bloch. See the boundary conditions section here for the complete list. The FDTD solver supports a number of different types of sources such as point dipoles, beams, plane waves, a total-field scattered-field (TFSF) source, a guided-mode source for integrated optical components, and an imported source to interface with external photonic design softwares. Detailed information about each of these sources is contained in the sources section.
You can find more information on Wikipedia.
Meshing
FDTD uses a rectangular, Cartesian style mesh, like the one shown in the following screenshot. It's important to understand that of the fundamental simulation quantities (material properties and geometrical information, electric and magnetic fields) are calculated at each mesh point. Obviously, using a smaller mesh allows for a more accurate representation of the device, but at a substantial cost. As the mesh becomes smaller, the simulation time and memory requirements will increase. The FDTD solver provides a number of features, including the conformal mesh algorithm, that allow you to obtain accurate results, even when using a relatively coarse mesh.
By default, the simulation mesh is automatically generated. To maintain accuracy, the meshing algorithm will create a smaller mesh in high index (to maintain a constant number of mesh points per wavelength) and highly absorbing (resolve penetration depths) materials. In some cases, it is also necessary to manually add additional meshing constraints. Usually, this involves forcing the mesh to be smaller near complex structures (often metal) where the fields are changing very rapidly.
Mesh refinement options, such as conformal mesh, can allow the solver to maintain high accuracy without the need for a very small mesh.
Integrating over lines, surfaces and volumes
The electric and magnetic fields are recorded on the finite-difference mesh, as shown below for a 2D monitor, where the grey dots represent the positions where the fields are recorded. The thick black outline shows the limits of the surface monitor as seen in the Layout Editor. This monitor has an x-span of 12dx and a y-span of 5dy. This monitor records a total of 13x6 data points.
A typical calculation with this monitor might be to integrate the total power flow across the surface of the monitor, or
$$Power(w) = \frac{1}{2}\int real{\{\overrightarrow{P}(w)\}} \cdot d \overrightarrow{S}$$
In order to calculate this quantity, we provide the scripting function integrate. If Pz is a variable of dimension 13x6x1, and x and y are the corresponding position vectors, then the desired quantity is:
power = 0.5*integrate(real(Pz),1:2,x,y);
The integrate script command can be used to integrate over spatial dimensions even when several frequency points have been recorded. For example, if Pz is a variable of dimension 13x6x1x10, representing 10 frequency points, then the following can be used to integrate over y (ie dimension 2) and then x (ie dimension 1)
power = 0.5*integrate(real(Pz),1:2,x,y);
power will be a matrix of dimension 10x1.
Units and normalization
Unless otherwise specified, all quantities are returned in SI units. Please see Units and normalization for more information.
References
1. Dennis M. Sullivan, Electromagnetic simulation using the FDTD method. New York: IEEE Press Series, (2000).
2. Allen Taflove, Computational Electromagnetics: The Finite-Difference Time-Domain Method. Boston: Artech House, (2005).
3. Stephen D. Gedney, Introduction to the Finite-Difference Time-Domain (FDTD) Method for Electromagnetics. Morgan & Claypool publishers, (2011).