Yee-update

The Yee scheme staggers electric and magnetic fields in space and time. Each curl is evaluated at the position of the field it updates.

A staggered grid

On a uniform Cartesian grid, one common indexing convention is:

Component	Position in grid units	Time
$E_x$	$(i+\tfrac12,j,k)$	$n\Delta t$
$E_y$	$(i,j+\tfrac12,k)$	$n\Delta t$
$E_z$	$(i,j,k+\tfrac12)$	$n\Delta t$
$H_x$	$(i,j+\tfrac12,k+\tfrac12)$	$(n+\tfrac12)\Delta t$
$H_y$	$(i+\tfrac12,j,k+\tfrac12)$	$(n+\tfrac12)\Delta t$
$H_z$	$(i+\tfrac12,j+\tfrac12,k)$	$(n+\tfrac12)\Delta t$

Electric components lie along cell edges; magnetic components lie at face centers. Translating the entire lattice gives equivalent conventions. BeamZ uses component-specific Yee supports; arrays for different components need not have identical shapes.

Staggering makes a centered derivative local. For example,

\left.\partial_y E_z\right|_{i,j+1/2,k+1/2} \approx\frac{E_z(i,j+1,k+1/2)-E_z(i,j,k+1/2)}{\Delta y}.

Both samples are half a cell from the evaluation point. In smooth media this derivative has second-order spatial accuracy. The staggered construction originates in Yee’s 1966 paper.

From derivatives to updates

Let $\nabla_h\times$ denote the discrete curl. Centering the time derivatives in the two evolution equations, first in a lossless medium, gives

\boxed{\mathbf H^{n+1/2}=\mathbf H^{n-1/2} -\frac{\Delta t}{\mu}\nabla_h\times\mathbf E^n}

\boxed{\mathbf E^{n+1}=\mathbf E^n +\frac{\Delta t}{\varepsilon} \left(\nabla_h\times\mathbf H^{n+1/2}-\mathbf J_{\mathrm s}^{n+1/2}\right)}.

Evaluate coefficients at the component being updated. Advance $\mathbf H$ and its corrections, then $\mathbf E$ and its corrections. For scalar nondispersive materials, this leapfrog requires no global linear solve.

A complete two dimensional example

For fields invariant along $z$ , the TMz system contains $E_z,H_x,H_y$ . Choose $E_z$ at $(i,j)$ , $H_x$ at $(i,j+1/2)$ , and $H_y$ at $(i+1/2,j)$ . For a source-free, lossless interior:

H_{x,i,j+1/2}^{n+1/2}=H_{x,i,j+1/2}^{n-1/2} -\frac{\Delta t}{\mu_{x,i,j+1/2}} \frac{E_{z,i,j+1}^{n}-E_{z,i,j}^{n}}{\Delta y},

H_{y,i+1/2,j}^{n+1/2}=H_{y,i+1/2,j}^{n-1/2} +\frac{\Delta t}{\mu_{y,i+1/2,j}} \frac{E_{z,i+1,j}^{n}-E_{z,i,j}^{n}}{\Delta x},

\begin{aligned} E_{z,i,j}^{n+1}=E_{z,i,j}^{n} +\frac{\Delta t}{\varepsilon_{z,i,j}}\bigg(& \frac{H_{y,i+1/2,j}^{n+1/2}-H_{y,i-1/2,j}^{n+1/2}}{\Delta x}\\ &-\frac{H_{x,i,j+1/2}^{n+1/2}-H_{x,i,j-1/2}^{n+1/2}}{\Delta y}\bigg). \end{aligned}

The subscripts on $\varepsilon_z,\mu_x,\mu_y$ identify the field supports. Even an isotropic material may have different effective coefficients on these supports at an interface. The complementary TEz system contains $H_z,E_x,E_y$ ; neither reduced system describes arbitrary three dimensional polarization coupling.

Conductivity

Center the Ohmic term at the midpoint using $(\mathbf E^{n+1}+\mathbf E^n)/2$ . Rearranging gives

\mathbf E^{n+1}=C_a\mathbf E^n+ C_b\left(\nabla_h\times\mathbf H^{n+1/2}-\mathbf J_{\mathrm s}^{n+1/2}\right),

C_a=\frac{1-\sigma\Delta t/(2\varepsilon)}{1+\sigma\Delta t/(2\varepsilon)}, \qquad C_b=\frac{\Delta t/\varepsilon}{1+\sigma\Delta t/(2\varepsilon)}.

Setting $\sigma=0$ recovers the lossless update. BeamZ precomputes these decay and curl coefficients. CPML introduces additional derivative corrections, described on the CPML page.

Stability and resolution

For the standard explicit Yee scheme with positive scalar, nondispersive materials on a uniform grid, the Courant limit is

\Delta t\leq\frac{1}{v_{\max} \sqrt{\Delta x^{-2}+\Delta y^{-2}+\Delta z^{-2}}}.

Omit invariant directions in one or two dimensions. With equal spacing $h$ in $d$ dimensions, $\Delta t\le h/(v_{\max}\sqrt d)$ . Use the fastest material wave speed in the domain, not the slower speed in a high-index core. A Courant factor below one provides margin.

Stability does not guarantee accuracy: numerical phase velocity depends on wavelength and direction. Resolve the shortest material wavelength and refine until the measured quantity converges. Interfaces can reduce convergence order.

BeamZ also supports rectilinear graded grids. There the derivatives use local edge or center distances and the solver computes a corresponding stability bound; replacing all those distances by an average spacing is incorrect.

Continue to Rasterization.