RC-FDTD Simulations

Recursive convolution finite difference time domain (RC-FDTD) simulations have long been used to numerically solve Maxwell’s equations. This simulation technique discretizes the time domain and evolves the electric and magnetic fields in time using a set of update equations. Within the simulation, space is discretized into intervals of length \(\Delta z\) and time into intervals of length \(\Delta t\). A specific point in time and space is accessed via \(z=i\Delta z\) and \(t=n\Delta t\). The simulation relies on a number of assumptions:

  • All materials are linear dielectrics such that \(P(z,\omega)=\epsilon_0\chi E(z,\omega)\).
  • The electric and magnetic fields are plane waves propagating along spatial coordinate \(z\).
  • Materials are uniform along spatial coordinates \(x\) and \(y\).
  • The electric and magnetic fields are zero for all time prior to the start of the simulation (\(E(z,t)=0\) for all \(t<0\)).
  • The electric field \(E(z,t)\) is approximately constant over all time intervals of duration \(\Delta t\).
  • The magnetization of all materials is zero (\(\vec{M}=\vec{0}\)).

These assumptions allow the derivation of the discretized displacement field \(D^{i,n}\). The displacement field \(\vec{D}(\vec{r},\omega)\), with the requirement that simulated materials are linear dielectrics such that \(P(z,\omega)=\epsilon_0\chi(z,\omega) E(z,\omega)\) and the requirement that the field varies over only the spatial coordinate \(z\) we find that \(D(z,\omega)\) is

\[D(z,\omega)=\epsilon_0\left[1+\chi (z,\omega)\right]E(z,\omega)\]

The displacement field \(D(z,\omega)\) can be transformed to the time domain via

\[\begin{split}D(z,t)=&\mathcal{F}^{-1}\left\{D(z,\omega)\right\} \\ =&\mathcal{F}^{-1}\left\{\epsilon_0\left[1+\chi (\omega)\right]E(z,\omega)\right\} \\ =&\mathcal{F}^{-1}\left\{\epsilon_0\mathcal{F}\left\{1+\chi (t)\right\}\mathcal{F}\left\{E(z,t)\right\}\right\}\end{split}\]

where \(\mathcal{F}\left\{a(t)\right\}\) and \(\mathcal{F}^{-1}\left\{a(\nu)\right\}\) to denote Fourier and inverse Fourier transforms. Thus via the convolution theorem

\[\begin{split}D(z,t)&=\mathcal{F}^{-1}\left\{\epsilon_0\mathcal{F}\left\{1+\chi (t)\right\}\mathcal{F}\left\{E(z,t)\right\}\right\} \\ &=\epsilon_0\left[1+\chi (t)\right]*\left[E(z,t)\right] \\ &=\epsilon_0\left[\epsilon_\infty E(z,t)+\int_0^t\chi (\tau)E(z,t-\tau) d\tau\right]\end{split}\]

where \(*\) denotes a convolution. It is assumed that \(E(z,t)=0\) for all \(t<0\). We discretize this result by replacing the \(z\) and \(t\) coordinates via \(z=i\Delta z\) and \(t=n\Delta t\) where \(i,n\in\mathbb{R}\), yielding

\[\begin{split}D(i\Delta z,n\Delta t)=&\epsilon_0\epsilon_\infty E(i\Delta z,n\Delta t) \\ &+\epsilon_0\int_0^{n\Delta t}\chi (\tau)E(i\Delta z,n\Delta -\tau) d\tau\end{split}\]

Assuming that \(E(i\Delta z,n\Delta -\tau)\) is constant over all time intervals of duration \(\Delta t\) the integral is replaced with a sum

\[D^{i,n}=\epsilon_0\epsilon_\infty E^{i,n}+\epsilon_0\sum^{n-1}_{m=0}E^{i,n-m}\chi ^m \label{eq:disp}\]

where

\[\chi ^m=\int_{m\Delta t}^{(m+1)\Delta t}\chi (\tau) d\tau\]

It is not assumed \(\chi(t)\) is constant over any time interval. This result is consistent with the result derived in Luebbers et al. and Beard et al..

\[D^{i,n}=\epsilon_0\epsilon_\infty E^{i,n}+\epsilon_0\sum^{n-1}_{m=0}E^{i,n-m}\chi ^m\]

where

\[\chi ^m=\int_{m\Delta t}^{(m+1)\Delta t}\chi (\tau) d\tau \label{eq:chi_conv}\]

This result is significant for the RC-FDTD simulation framework implmented here as it means a material can be simulated as long as one can define $chi(t)$ for that material.

We proceed by deriving the update equations for the electric and magnetic fields. With the requirement that \(\vec{M}=\vec{0}\) and the requirement that the electric and magnetic fields are uniform in spatial coordinates \(x\) and \(y\), Faraday’s law of induction and Ampere’s law with Maxwell’s addition reduce to

\[\frac{\partial E}{\partial z}=-\mu_0\frac{\partial H}{\partial t} \qquad -\frac{\partial H}{\partial z}=I_f+\frac{\partial D}{\partial t}\]

where \(I_f\) is along \(\hat{z}\). Noting the definition of a derivative we find

\[\begin{split}\lim_{\Delta z\to0}\frac{E(z+\Delta z,t)-E(z,t)}{\Delta z}=-\mu_0\lim_{\Delta t\to0}\frac{H(z,t+\Delta t)-H(z,t)}{\Delta t} \\ -\lim_{\Delta z\to0}\frac{H(z+\Delta z,t)-H(z,t)}{\Delta z}=I_f+\lim_{\Delta t\to0}\frac{D(z,t+\Delta t)-D(z,t)}{\Delta t}\end{split}\]

From here the discretization process is simple. We simply remove each limit from the equations, define an appropriate value of \(\Delta z\) and \(\Delta t\), and replace the fields with their discretized forms.

\[\begin{split}\frac{E^{i+1,n}-E^{i,n}}{\Delta z}=-\mu_0\frac{H^{i,n+1}-H^{i,n}}{\Delta t} \label{eq:faraday} \\ -\frac{H^{i+1,n}-H^{i,n}}{\Delta z}=I_f+\frac{D^{i,n+1}-D^{i,n}}{\Delta t} \label{eq:ampere}\end{split}\]

If \(\Delta z\) and \(\Delta t\) aren’t small enough such that the derivative is accurate then the RC-FDTD simulation will break down.

We solve Eq.(ref{eq:faraday}) for \(H^{i,n+1}\), finding the following update equation

\[H^{i,n+1}=H^{i,n}-\frac{1}{\mu_0}\frac{\Delta t}{\Delta z}\left[E^{i+1,n}-E^{i,n}\right]\]

In order to solve Eq.(ref{eq:ampere}) we use the result of Eq.(ref{eq:disp}) to determine \(D^{i,n+1}-D^{i,n}\) in terms of \(E^{i+1,n}\) and \(E^{i,n}\)

\[\begin{split}D^{i,n+1}-D^{i,n}&=\epsilon_0\epsilon_\infty E^{i,n+1}+\epsilon_0\sum^{n}_{m=0}E^{i,n+1-m}\chi ^m-\epsilon_0\epsilon_\infty E^{i,n}-\epsilon_0\sum^{n-1}_{m=0}E^{i,n-m}\chi ^m \\ &=\epsilon_0\epsilon_\infty\left[E^{i,n+1}-E^{i,n}\right]+\epsilon_0\left[\sum^{n}_{m=0}E^{i,n+1-m}\chi ^m-\sum^{n-1}_{m=0}E^{i,n-m}\chi ^m\right]\end{split}\]

Noting that

\[\begin{split}\sum^{n}_{m=0}E^{i,n+1-m}\chi ^m-\sum^{n-1}_{m=0}E^{i,n-m}\chi ^m&=E^{i,n+1}\chi ^0+\sum^{n}_{m=1}E^{i,n+1-m}\chi ^m-\sum^{n-1}_{m=0}E^{i,n-m}\chi ^m \\ &=E^{i,n+1}\chi ^0+\sum^{n-1}_{m=0}E^{i,n+1-(m+1)}\chi ^{m+1}-\sum^{n-1}_{m=0}E^{i,n-m}\chi ^m \\ &=E^{i,n+1}\chi ^0+\sum^{n-1}_{m=0}E^{i,n-m}\left[\chi ^{m+1}-\chi ^m\right] \\\end{split}\]

and letting

\[\begin{split}\Delta\chi^m&=\chi^m-\chi^{m+1} \\ \psi^n&=\sum^{n-1}_{m=0}E^{i,n-m}\Delta\chi^m\end{split}\]

we find

\[\begin{split}D^{i,n+1}-D^{i,n}&=\epsilon_0\epsilon_\infty\left[E^{i,n+1}-E^{i,n}\right]+\epsilon_0\left[E^{i,n+1}\chi^0-\psi^n\right] \\ &=\epsilon_0\left[\epsilon_\infty+\chi^0\right]E^{i,n+1}-\epsilon_0\epsilon_\infty E^{i,n}-\epsilon_0\psi^n\end{split}\]

Substituting this result into Eq.(ref{eq:ampere}) and solving for \(E^{i,n+1}\) we find

\[\begin{split}E^{i,n+1}=&\frac{\epsilon_\infty}{\epsilon_\infty+\chi^0}E^{i,n}+\frac{1}{\epsilon_\infty+\chi^0}\psi^n-\frac{\Delta tI_f}{\epsilon_0\left[\epsilon_\infty+\chi^0\right]} \\ &-\frac{1}{\epsilon_0\left[\epsilon_\infty+\chi^0\right]}\frac{\Delta t}{\Delta z}\left[H^{i+1,n}-H^{i,n}\right]\end{split}\]

We then implement the Yee cell in the simulation by offsetting the electric and magnetic field cells by half a spatial and temporal incrementcite{beard}, producing

\[\begin{split}H^{i+1/2,n+1/2}=&H^{i+1/2,n-1/2}-\frac{1}{\mu_0}\frac{\Delta t}{\Delta z}\left[E^{i+1,n}-E^{i,n}\right] \label{eq:hup} \\ E^{i,n+1}=&\frac{\epsilon_\infty}{\epsilon_\infty+\chi^0}E^{i,n}+\frac{1}{\epsilon_\infty+\chi^0}\psi^n-\frac{\Delta tI_f}{\epsilon_0\left[\epsilon_\infty+\chi^0\right]} \\ &-\frac{1}{\epsilon_0\left[\epsilon_\infty+\chi^0\right]}\frac{\Delta t}{\Delta z}\left[H^{i+1/2,n+1/2}-H^{i-1/2,n+1/2}\right]\end{split}\]

where

\[\begin{split}\Delta\chi^m&=\chi^m-\chi^{m+1} \nonumber \\ \psi^n&=\sum^{n-1}_{m=0}E^{i,n-m}\Delta\chi^m \label{eq:psi}\end{split}\]

The accuracy of the derivative approximation inherent to these update equations relies on choosing some \(\Delta z\) and \(\Delta t\) small enough such that the electric and magnetic fields are approximately linear over spatial intervals \(\Delta z\) and time intervals \(\Delta t\). If this condition is not met then the accuracy of the derivative approximation breaks down. The update equations derived here are significant as they reveal that any linear dielectric can be accurately simulated via the RC-FDTD method as long as the electric susceptibility of the material \(\chi(t)\) is well defined. We turn our attention to modeling the electric susceptibility of materials in section ref{sec:susceptibility}.

rcfdtdpy provides a framework in which the user need only provide the electric susceptibility \(\chi(t)\) to run a simulation.