Der

In the realm of computational electromagnetics, we face a fundamental paradox: the laws of physics are described by continuous partial differential equations (PDEs), yet the machines we use to solve them are inherently discrete. Computers cannot process an infinite continuum of points; they can only manipulate finite sets of numbers. Therefore, the bridge between theoretical electromagnetics and practical simulation lies in discretization—the process of transforming continuous field equations into a solvable system of algebraic equations.

The most influential framework for this transformation is the Finite-Difference Time-Domain (FDTD) method. At its heart lies the ingenious discretization scheme proposed by Kane Yee, which allows us to track the evolution of electromagnetic waves through space and time with remarkable efficiency and accuracy.
To simulate the propagation of electromagnetic waves in an isotropic, linear, and lossless medium, we rely on the differential form of Maxwell's equations. While there are four fundamental equations, the temporal evolution of the fields is primarily driven by the two curl equations:

  1. Faraday’s Law of Induction:
    $$\nabla \times \mathbf{E} = -\mu \frac{\partial \mathbf{H}}{\partial t}$$
  2. Ampère’s Circuital Law:
    $$\nabla \times \mathbf{H} = \epsilon \frac{\partial \mathbf{E}}{\partial t} + \mathbf{J}$$

Here, $\mathbf{E}$ and $\mathbf{H}$ represent the electric and magnetic field vectors, $\epsilon$ and $\mu$ are the permittivity and permeability of the medium, and $\mathbf{J}$ is the current density. The Gauss laws ($\nabla \cdot \mathbf{E} = \rho/\epsilon$ and $\nabla \cdot \mathbf{H} = 0$) ensure the divergence-free nature of the fields but are often implicitly satisfied by the curl equations in a well-posed FDTD implementation.

The Yee Discretization Strategy

A naive approach to discretization—placing all field components on a standard Cartesian grid and using simple central differences—often leads to numerical instability or a failure to capture the essential "curl" relationship between the fields. To solve this, the Yee algorithm employs a staggered grid approach, decoupling the fields in both space and time.

Spatial Staggering: The Yee Cell

In a Yee grid, the components of the electric and magnetic fields are not co-located. Instead, they are offset such that the components required to calculate a curl at a specific point are positioned exactly where they are needed.

  • The Electric field components ($E_x, E_y, E_z$) are located on the edges of the grid cells. For instance, $E_x$ is positioned at $(i+1/2, j, k)$.
  • The Magnetic field components ($H_x, H_y, H_z$) are located on the faces of the cells, offset from the electric field. For example, $H_y$ is positioned at $(i, j+1/2, k)$.

This spatial arrangement ensures that the discrete approximation of the curl operator is naturally centered and spatially accurate.

Temporal Staggering: The Leapfrog Method

To achieve second-order accuracy in time, the FDTD method utilizes a "leapfrog" scheme. Rather than updating both fields at the same time step, we stagger them by a half-step ($\Delta t/2$):

  • The Magnetic field $\mathbf{H}$ is updated at integer time steps ($n, n+1, n+2, \dots$).
  • The Electric field $\mathbf{E}$ is updated at half-integer time steps ($n+1/2, n+3/2, \dots$).

This temporal offset allows us to use central difference approximations for the time derivatives, which significantly enhances the accuracy of the simulation.

Mathematical Derivation of the Update Equations

To understand how this works in practice, let us derive the update equation for the $y$-component of the magnetic field using Faraday's Law.

1. Discretizing the Spatial Curl

We begin with the $y$-component of the curl of $\mathbf{E}$:
$$(\nabla \times \mathbf{E})_y = \frac{\partial E_x}{\partial z} - \frac{\partial E_z}{\partial x}$$

Using a central difference approximation on the Yee grid, the spatial derivatives at the location of $H_y$ are expressed as:
$$\left. \frac{\partial E_x}{\partial z} \right|{i,j,k+1/2} \approx \frac{E_x(i, j, k+1) - E_x(i, j, k)}{\Delta z}$$
$$\left. \frac{\partial E_z}{\partial x} \right|{i,j,k+1/2} \approx \frac{E_z(i+1, j, k) - E_z(i, j, k)}{\Delta x}$$

2. Discretizing the Time Derivative

Because the curl of $\mathbf{E}$ is defined at the half-integer time step $n+1/2$, while $\mathbf{H}$ is defined at integer steps, the time derivative for $H_y$ is approximated as:
$$\left. \frac{\partial H_y}{\partial t} \right|_{n+1/2} \approx \frac{H_y^{n+1} - H_y^n}{\Delta t}$$

3. The Explicit Update Formula

By substituting these approximations back into Faraday's Law, we derive the explicit update equation for $H_y$:
$$H_y^{n+1} = H_y^n - \frac{\Delta t}{\mu} \left[ \frac{E_x^{n+1/2}(i, j, k+1) - E_x^{n+1/2}(i, j, k)}{\Delta z} - \frac{E_z^{n+1/2}(i+1, j, k) - E_z^{i, j, k)}{\Delta x} \right]$$

A symmetric derivation using Ampère's Law yields the update equation for the electric field $\mathbf{E}$, allowing the fields to "leapfrog" each other through time.

Algorithmic Execution Flow

In a computational environment, the FDTD algorithm follows a highly structured iterative loop:

  1. Initialization: Define the grid dimensions, material properties ($\epsilon, \mu$), and initial field distributions (usually zero). Set up boundary conditions.
  2. Magnetic Field Update: Calculate the $\mathbf{H}$ field at time $n+1$ using the $\mathbf{E}$ field from time $n+1/2$.
  3. Electric Field Update: Calculate the $\mathbf{E}$ field at time $n+1/2$ using the newly updated $\mathbf{H}$ field.
  4. Boundary Condition Application: Apply specialized algorithms at the edges of the simulation domain, such as Perfectly Matched Layers (PML) to simulate open space or reflective boundaries for closed cavities.
  5. Iteration: Repeat steps 2 through 4 until the desired simulation time is reached.

Numerical Stability and the CFL Condition

Because the FDTD method is an explicit scheme, it is not unconditionally stable. The time step $\Delta t$ cannot be chosen arbitrarily; it is strictly constrained by the spatial grid size and the speed of light in the medium. This constraint is known as the Courant-Friedrichs-Lewy (CFL) condition.

For a 3D Cartesian grid, the condition is:
$$\Delta t \le \frac{1}{c \sqrt{\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} + \frac{1}{\Delta z^2}}}$$

If $\Delta t$ exceeds this threshold, the numerical information travels faster than the physical wave can be resolved by the grid, leading to exponential growth of errors and a total divergence of the simulation. Physically, this ensures that the "numerical domain of dependence" encompasses the "physical domain of dependence."

Conclusion

The derivation of the FDTD update equations represents a masterclass in numerical engineering. By cleverly staggering fields in both space and time via the Yee cell, we transform complex, continuous differential equations into a series of simple, efficient, and second-order accurate arithmetic operations. This approach remains the bedrock of modern electromagnetic solvers, enabling everything from antenna design to the analysis of complex photonic structures.