Efficient Computation of High-Order Electromagnetic Field Derivatives for Multiple Design Parameters in FDTD

TL;DR

This paper presents a highly efficient and accurate computational framework for deriving electromagnetic field derivatives of any order with respect to multiple design parameters using FDTD, requiring minimal simulations.

Contribution

The authors develop a novel FDTD-based method that computes high-order derivatives with respect to multiple parameters using only a few simulations, ensuring guaranteed accuracy.

Findings

Requires only one or two FDTD simulations for derivatives of any order.

Provides subtractive cancellation error-free derivatives.

Enables precise sensitivity analysis and uncertainty quantification.

Abstract

This paper introduces a new computational framework to derive electromagnetic field derivatives with respect to multiple design parameters up to any order with the Finite-Difference Time-Domain (FDTD) technique. Specifically, only one FDTD simulation is needed to compute the first-order field derivatives with respect to N parameters, while two FDTD simulations are needed to compute the field derivatives with respect to one parameter up to any order. The field derivatives with respect to N parameters up to any order are computed with (N+1) FDTD runs. In addition to its efficiency, this framework is based on a subtractive cancellation error-free approach, providing guaranteed accuracy toward the computation of field derivatives up to any order. With high-order field derivatives available, sensitivity analysis, parametric modelling and uncertainty quantification can be accurately performed.

Figures

Tables

Table 1. TABLE I: Operation Count for The Computation of Field Derivatives in FDTD simulations

	${\nabla }_{𝝃}{E}_{p_a}^z$	$\frac{{\partial }^M{E}_{p_a}^z}{\partial {\xi }^M}$	$\frac{{\partial }^M{E}_{p_a}^z}{\partial {𝝃}^M}$	$ℍ\left(E_{p_a}^z\right)$
This paper	$1$	$2$	$N+1$	$N+1$
CFD-FDTD	$2N$	$\left(M+1\right)$	$N\left(M+1\right)$	$2N^2+N$

Equations

$$\frac{\partial E^{x,y,z}}{\partial \xi }({\xi }_0)=\frac{E^{x,y,z}({\xi }_0+h)-E^{x,y,z}({\xi }_0-h)}{2h}+\mathcal{O}(h^2).$$

$$E^{x,y,z}(\xi )=E^{x,y,z}({\xi }_0)+\sum _{n=1}^{\infty }\frac{1}{n!}\frac{{\partial }^n{E}^{x,y,z}({\xi }_0)}{\partial {\xi }^n}(\xi -{\xi }_0{)}^n.$$

$$\text{E}=\left[\begin{array}{c}{\text{E}}^x\\ {\text{E}}^y\\ {\text{E}}^z\end{array}\right],\text{H}=\left[\begin{array}{c}{\text{H}}^x\\ {\text{H}}^y\\ {\text{H}}^z\end{array}\right],$$

$${\text{E}}^x=\left[\begin{array}{c}{E}_{1,1,1}^x\\ E_{2,1,1}^x\\ ⋮\\ E_{N_x,N_y,N_z}^x\end{array}\right]={\left[E^x\right]}_p\in {\mathbb{R}}^{P\times 1}$$

$$\begin{gathered} \mathbb{P}(i,j,k)=(k-1)N_x{N}_y+(j-1)N_x+i, \\ \forall i=1,2,..,N_x,j=1,2,..,N_y,k=1,2,...,N_z \\ =\{p\mid p=1,...,P\}, \end{gathered}$$

$${\text{d}}^x=\left[\begin{array}{c}1/\Delta x_1\\ 1/\Delta x_2\\ ⋮\\ 1/\Delta x_P\end{array}\right],{\text{d}}^y=\left[\begin{array}{c}1/\Delta y_1\\ 1/\Delta y_2\\ ⋮\\ 1/\Delta y_P\end{array}\right],{\text{d}}^z=\left[\begin{array}{c}1/\Delta z_1\\ 1/\Delta z_2\\ ⋮\\ 1/\Delta z_P\end{array}\right]\in {\mathbb{R}}^{P\times 1}$$

$$\widetilde{\Delta y_p}({\xi }_1)=\Delta y_p+\delta {\xi }_1=\Delta y_p+({\xi }_1-\overline{{\xi }_1}).$$

$${\text{d}}^y={\left[d^y\right]}_p=\{\begin{array}{l}1/\widetilde{\Delta y_p}=1/(\Delta y_p+\delta {\xi }_1),p\in P_{{\xi }_1}\\ 1/\Delta y_p,\text{otherwise.}\end{array}$$

$$\{\begin{array}{ll}{\text{E}}^{n+1}& ={\text{E}}^n+{\text{A}}_e{\text{H}}^{n+\frac{1}{2}}\\ {\text{H}}^{n+\frac{1}{2}}& ={\text{H}}^{n-\frac{1}{2}}+{\text{A}}_h{\text{E}}^n\end{array},$$

$${\text{A}}_e$$

$${\text{A}}_h$$

$$\text{D}\equiv \left[\begin{array}{ccc}\text{0}& {\text{D}}^z& {\text{D}}^y\\ {\text{D}}^z& \text{0}& {\text{D}}^x\\ {\text{D}}^y& {\text{D}}^x& \text{0}\end{array}\right]$$

$${\text{D}}^y=\text{diag}({\text{d}}^y)={\left[D^y\right]}_{ij}=\{\begin{array}{ll}1/\widetilde{\Delta y_i},& \text{if i=j∈Pξ1}\\ 1/\Delta y_i& \text{if i=j∈/Pξ1}\\ 0,& \text{otherwise}.\end{array}$$

$$\{\begin{array}{l}\frac{\partial {\text{E}}^{n+1}}{\partial {\xi }_1}=\frac{\partial {\text{E}}^n}{\partial {\xi }_1}+{\text{A}}_e\frac{\partial {\text{H}}^{n+\frac{1}{2}}}{\partial {\xi }_1}+\mathcal{J}\\ \frac{\partial {\text{H}}^{n+\frac{1}{2}}}{\partial {\xi }_1}=\frac{\partial {\text{H}}^{n-\frac{1}{2}}}{\partial {\xi }_1}+{\text{A}}_h\frac{\partial {\text{E}}^n}{\partial {\xi }_1}+\mathcal{K}\end{array}$$

$$\mathcal{J}=\frac{\partial {\text{A}}_e}{\partial {\xi }_1}{\text{H}}^{n+\frac{1}{2}},\mathcal{K}=\frac{\partial {\text{A}}_h}{\partial {\xi }_1}{\text{E}}^n$$

$$\{\begin{array}{l}E(x,y,z,t;{\xi }_1)=\sum _{i,j,k,n}{E}_{i,j,k}^n({\xi }_1)h_i(x)h_j(y)h_k(z)h_n(t)\\ H(x,y,z,t;{\xi }_1)=\sum _{i,j,k,n}{H}_{i,j,k}^n({\xi }_1)h_i(x)h_j(y)h_k(z)h_n(t),\end{array}$$

$$\{\begin{array}{l}\frac{\partial E(x,y,z,t;{\xi }_1)}{\partial {\xi }_1}=\sum _{i,j,k,n}\frac{\partial E_{i,j,k}^n({\xi }_1)}{\partial {\xi }_1}{h}_i(x)h_j(y)h_k(z)h_n(t)\\ \frac{\partial H(x,y,z,t;{\xi }_1)}{\partial {\xi }_1}=\sum _{i,j,k,n}\frac{\partial H_{i,j,k}^n({\xi }_1)}{\partial {\xi }_1}{h}_i(x)h_j(y)h_k(z)h_n(t),\end{array}$$

$$h_i(x)=h(\frac{x}{\Delta x_i}-i),\text{and}$$

$$h(\tau )=\{\begin{array}{l}1,\frac{-1}{2}\le \tau \le \frac{1}{2}\\ 0,\text{otherwise}.\end{array}$$

$$\frac{\partial {\text{A}}_{e,h}}{\partial {\xi }_1}=\Delta t{\text{M}}_{\epsilon ,\mu }^{-1}({\text{C}}_{e,h}\circ \frac{\partial \text{D}}{\partial {\xi }_1}).$$

$$\begin{array}{rl}\frac{\partial \text{D}}{\partial {\xi }_1}& =\left[\begin{array}{ccc}\text{0}& \partial {\text{D}}^z/\partial {\xi }_1& \partial {\text{D}}^y/\partial {\xi }_1\\ \partial {\text{D}}^z/\partial {\xi }_1& \text{0}& \partial {\text{D}}^x/\partial {\xi }_1\\ \partial {\text{D}}^y/\partial {\xi }_1& \partial {\text{D}}^x/\partial {\xi }_1& \text{0}\end{array}\right]\\ & =\left[\begin{array}{ccc}\text{0}& \text{0}& \partial {\text{D}}^y/\partial {\xi }_1\\ \text{0}& \text{0}& \text{0}\\ \partial {\text{D}}^y/\partial {\xi }_1& \text{0}& \text{0}\end{array}\right].\end{array}$$

$$\begin{gathered} \frac{\partial \left(\frac{1}{\widetilde{\Delta y_p}}\right)}{\partial {\xi }_1}=\frac{\partial \left(\frac{1}{\Delta y_p+{\xi }_1-\overline{{\xi }_1}}\right)}{\partial {\xi }_1} \\ =\frac{-1}{{\left(\Delta y_p+{\xi }_1-\overline{{\xi }_1}\right)}^2}\approx \frac{-1}{{\Delta y_p}^2}. \end{gathered}$$

$$\frac{\partial {\text{D}}^y}{\partial {\xi }_1}={\left[\frac{\partial D^y}{\partial {\xi }_1}\right]}_{ij}=\{\begin{array}{ll}1/\Delta y_i^2,& \text{if i=j∈Pξ1},\\ 0,& \text{otherwise}.\end{array}$$

$$\mathcal{J}=\left[\begin{array}{c}{\mathcal{J}}^x\\ {\mathcal{J}}^y\\ {\mathcal{J}}^z\end{array}\right],\mathcal{K}=\left[\begin{array}{c}{\mathcal{K}}^x\\ {\mathcal{K}}^y\\ {\mathcal{K}}^z\end{array}\right],$$

$${\nabla }_{\xi }{E}_{p_a}^z\equiv \left[\begin{array}{cccc}\frac{\partial E_{p_a}^z}{\partial {\xi }_1}& \frac{\partial E_{p_a}^z}{\partial {\xi }_2}& \dots & \frac{\partial E_{p_a}^z}{\partial {\xi }_N}\end{array}\right],$$

$$\begin{split}&\textbf{d}^{x}=\left[d^{x}\right]_{p}=\\ &\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern 2.5pt\kern-2.0pt\left[\kern 0.0pt\kern-2.5pt\kern-5.55557pt\vbox{\kern-0.86108pt\vbox{\vbox{ \halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep& \kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&& \kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr 2.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle 1$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle 2$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle p_{1}$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle\dots$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle p_{2}$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle\dots$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle p_{N}$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle\dots$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle P\\$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle\frac{1}{\Delta x_{1}}$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle\frac{1}{\Delta x_{2}}$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle\frac{1}{\Delta x_{p_{1}}+\delta_{1}}$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle\dots$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle\frac{1}{\Delta x_{p_{2}}+\delta_{2}}$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle\dots$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle\frac{1}{\Delta x_{p_{N}}+\delta_{N}}$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle\dots$\hfil\kern 2.0pt&2.0pt\hfil$\scriptstyle\frac{1}{\Delta x_{P}}$\hfil\kern 2.0pt\crcr}}}}\right]$}}\end{split}$$

$$\mathfrak{D}_{\bm{\xi}}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern 2.5pt\kern-5.0pt\left[\kern 0.0pt\kern-2.5pt\kern-5.55557pt\vbox{\kern-0.86108pt\vbox{\vbox{ \halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep& \kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&& \kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr 5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 1$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 2$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle p_{1}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\dots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle p_{2}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\dots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle p_{N}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\dots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle P\\1$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0\\2$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0\\p_{1}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\frac{\partial}{\partial\xi_{1}}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0\\\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots\\p_{2}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\frac{\partial}{\partial\xi_{2}}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0\\\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots\\p_{N}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\frac{\partial}{\partial\xi_{N}}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\\\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots\\P$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt\crcr}}}}\right]$}}$$

$${\mathfrak{D}}_{\mathbf{\xi }}{\text{D}}^x={\mathfrak{D}}_{\mathbf{\xi }}{\left[D^x\right]}_{ij}=\{\begin{array}{l}-1/{\Delta x_i}^2,i=j=p_1,p_2,...,p_N\\ 0,\text{otherwise}.\end{array}$$

$$\begin{array}{rl}\mathcal{J}& =\Delta t{\text{M}}_{\epsilon }^{-1}\left(\text{C}e\circ \left[\begin{array}{ccc}\text{0}& \text{0}& \text{0}\\ \text{0}& \text{0}& {\mathfrak{D}}_{\mathbf{\xi }}{\text{D}}^x\\ \text{0}& {\mathfrak{D}}_{\mathbf{\xi }}{\text{D}}^x& \text{0}\end{array}\right]\right)\text{H}\\ & ={\left[\begin{array}{c}{\mathcal{J}}^x{\mathcal{J}}^y{\mathcal{J}}^z\end{array}\right]}^T\end{array},$$

$$\begin{gathered} {\mathcal{J}}^z(t)={\left[{\mathcal{J}}^z(t)\right]}_p \\ =\{\begin{array}{l}\frac{-\Delta t}{{\epsilon }_0{\epsilon }_r{\Delta x_p}^2}\left(H_p^y(t)-H_{p-1}^y(t)\right),p=p_1,p_2,...,p_N\\ 0,\text{otherwise}.\end{array} \end{gathered}$$

$$\frac{\partial \widehat{E_{p_a}^z}}{\partial {\xi }_1}(f)={\mathcal{H}}_{p_a{p}_1}(f)\cdot \widehat{{\mathcal{J}}_{p_1}^{z,1}}(f),$$

$${\mathcal{H}}_{p_a{p}_n}(f)=\frac{\mathfrak{F}\left[E_p^z(t)\right]}{\mathfrak{F}\left[E_{p_a}^z(t)\right]},p_n=p_1,p_2,...,p_N,$$

Full text

I Introduction

As the design complexity of microwave circuits and systems increases, the sensitivity analysis with respect to a large number of parameters becomes an essential part in the design process [1, 2, 3]. Therefore, numerical techniques are in need to find not only the field solution to a given problem, but also its sensitivity under fabrication tolerances and other statistical uncertainties in geometric and material parameters [4, 5].

Notably, sensitivity of any field-based output function of interest can be derived if the corresponding derivatives of the electromagnetic field components are available. For example, the derivatives of the scattering parameters with respect to a design parameter $\xi$, written as $\partial S_{i,j}/\partial\xi$, can be derived by appyling the chain rule to the derivatives of the incident, reflected and transmitted fields at the ports of a microwave structure [6].

The standard approach to compute field derivatives is the central finite difference (CFD) method with second-order accuracy. For example, the electric field derivatives of the form $\partial E^{x,y,z}/\partial\xi$ over nominal value $\xi_{0}$ is approximated as follows:

[TABLE]

However, this CFD approximation method is subject to subtractive cancellation errors when a small step size $h$ is used [7]. In techniques such as FDTD, the geometric sensitivity analysis is often performed by assigning $h$ to the size of the cells corresponding to the geometry, and the step size dilemma is encountered in choosing a small $h$ to minimize the truncation error in (1) and the dispersion error in mesh distortion, while avoiding subtractive cancellation errors. Moreover, consider a case where the field derivatives with respect to an $N$-dimensional vector of parameters $\bm{\xi}=\left[\xi_{1},\xi_{2},\dots,\xi_{N}\right]^{T}$ are of interest. Then, CFD-based methods require at least $2N$ simulations to evaluate these derivatives. Hence, the computation overhead of CFD methods becomes significant when a large number of parameters is considered.

Alternatively, adjoint variable methods (AVM) have been combined with FDTD in [8] [9]. The AVM-FDTD method computes the gradient of an objective function with respect to $\bm{\xi}$ by running one additional FDTD simulation, which is a time-reversed version of the original system [10][11]. Recently, the AVM-FDTD method was extended to compute second-order derivatives [12] by applying FDTD to the wave equation rather than the first-order Maxwell’s equations. Indeed, second and higher-order derivatives are very useful in electromagnetic studies such as parametric modelling and optimization of a system [13, 14, 15]. For example, the electric field components can be expanded in a Taylor series with respect to a parameter $\xi$ as:

[TABLE]

Evidently, parametric design optimization as well as uncertainty quantification analysis can be performed using (2) [2, 16].

To compute field derivatives up to any order accurately, an FDTD-based approach to calculate derivatives without subtractive cancellation errors is proposed in [17] using the complex-step derivative (CSD) approximation method. Despite its accuracy, robustness and the compatibility to exiting FDTD codes, the resulting CSD-FDTD method produces a significant computation overhead in memory and execution time when multiple design parameters are considered. On the other hand, [18] introduced an alternative approach, by directly differentiating the FDTD update equations. As discussed in [18], the $M$-th order field derivatives can be found recursively from the $(M-1)$-th order field derivative. Hence, the computation overhead scales with the order of derivatives. For example, a total of $M$ simulations is needed to compute the $M$-th order field derivatives with respect to one parameter.

This paper builds on and significantly extends the work reported in [18] and introduces a new framework to compute high-order multi-parametric field derivatives accurately and efficiently. With the proposed approach:

• •

First-order field derivatives with respect to $N$ design parameters are found in a single FDTD simulation.

• •

Field derivatives up to order $M$ with respect to one parameter require two FDTD simulations.

• •

Field derivatives up to order $M$ with respect to $N$ parameters are found in $(N+1)$ FDTD simulations.

This is accomplished by mapping the FDTD equations with respect to field derivatives to the original FDTD simulation, driven by “equivalent sources”, placed at Yee cells related to the parameters under consideration. It should be noted that, all the electromagnetic fields in this paper are computed with original FDTD simulations without any mesh distortion.

II FDTD Update Equations for Field Derivatives

II-A Matrix form of FDTD update equations

The FDTD update equations are represented in matrix form throughout this paper. In a three dimensional FDTD domain consisting of $N_{x}\times N_{y}\times N_{z}=P$ cells, the electric and magnetic field components in $x,y,z$-directions at all Yee cells are assembled in two vectors E and H:

[TABLE]

where $\textbf{E}^{x,y,z}$ and $\textbf{H}^{x,y,z}$ are six vectors representing the electric and magnetic fields at all cells. For example:

[TABLE]

The index notation $p$ of $\textbf{E}^{x,y,z}$ and $\textbf{H}^{x,y,z}$ vectors is a transformation of the triad $(i,j,k)$ used in 3-D FDTD simulations to indicate the location of a Yee cell. The index $p$ is mapped to the triad $(i,j,k)$:

[TABLE]

where $\mathbb{P}$ is the set of all the cell indices in the simulation domain.

In addition to the field vectors, three vectors $\textbf{d}^{x,y,z}$ are defined to represent the reciprocal of the size of each cell. In an FDTD simulation of the geometry under consideration, $\textbf{d}^{x,y,z}$ are vectors of length $M$ with entries $1/\Delta x_{p}$, $1/\Delta y_{p}$ and $1/\Delta z_{p}$ respectively:

[TABLE]

In this paper, it is assumed that the size of a set of cells depends on a specific geometric parameter. That is, the $\Delta x,y,z$ of these cells are functions of a parameter. Consider the layout of a 3-D microstrip filter with three cascaded stubs, shown in Fig. 1. The lengths of the stubs are set to nominal values $\bar{\xi_{1}},\bar{\xi_{2}}$ and $\bar{\xi_{3}}$ at first. If the length of the first stub is a variable $\xi_{1}$ and deviated by $\delta\xi_{1}$, then $\xi_{1}=\bar{\xi_{1}}+\delta\xi_{1}$. Keeping $\delta\xi_{1}\rightarrow 0$, only the size of the cells coloured red in the inset of Fig. 1 would be affected by this perturbation. Hence, their new length $\tilde{\Delta y_{p}}$ in the $y$-direction is written as a function of $\xi_{1}$:

[TABLE]

These cells correspond to the end of the first stub and form a subset region $P_{\xi_{1}}$ of the FDTD simulation domain, i.e., $P_{\xi_{1}}\subset\mathbb{P}$. In this example, the vector $\textbf{d}^{y}$ is modified and written as:

[TABLE]

In this particular example, the $x$ and $y$-component of the size of the cells of subset $P_{\xi_{1}}$ are independent of $\xi_{1}$. Hence, $\textbf{d}^{x,z}$ are not modified.

II-B FDTD update equations for first-order field derivatives

Consider the 3-D microstrip filter shown in Fig.1. The standard FDTD update equations in matrix form are [19]:

[TABLE]

where $\textbf{A}_{e}$ and $\textbf{A}_{h}$ are the update matrices for E and H respectively:

[TABLE]

The $\circ$ operator denotes the element-wise multiplication of matrices; $\Delta t$ is the time-step used in the FDTD simulation; $\textbf{M}_{\mu}$ and $\textbf{M}_{\varepsilon}$ are the diagonal matrices representing material properties. $\textbf{C}e$ and $\textbf{C}h$ are two sparse matrices defining the discrete curl operators in matrix form with non-zero entities of either $+1$ or $1$. D is a matrix referring to the size of cells and defined as:

[TABLE]

The block matrices $\textbf{D}^{x,y,z}$ in D are diagonal matrices assembled from vectors $\textbf{d}^{x,y,z}$. Following the derivation in the previous subsection, $P_{\xi_{1}}$ is the subset including all the cells with the size depends on small change in $\xi_{1}$. $\textbf{D}^{y}$ is written as:

[TABLE]

Again, since the small change or perturbation in $\xi_{1}$ affects the cell size in the $y$-direction only, $\textbf{D}^{x}=\text{diag}(\textbf{d}^{x})$ and $\textbf{D}^{z}=\text{diag}(\textbf{d}^{z})$ are two diagonal matrices with non-zero, constant entities $1/\Delta x_{p}$ and $1/\Delta z_{p}$ respectively.

To find the field derivatives with respect to $\xi_{1}$, differentiation is performed directly to (37), yielding:

[TABLE]

Eq. (14) is the FDTD update equation with respect to first-order field derivatives. This differentiated form of FDTD update equations builds the connection between the computation of fields $\left\{\textbf{E},\textbf{H}\right\}$ and their derivatives. Indeed, $\left\{\textbf{E},\textbf{H}\right\}$ are discrete space-time fields, but also functions of $\xi_{1}$, assumed to be smooth for their derivatives to be well defined [20]. Hence:

[TABLE]

and their derivatives are derived as:

[TABLE]

where $h_{i,j,k}$ are pulse functions. For example:

[TABLE]

Further examination of (14) reveals that the electric and magnetic field derivatives are updated in time, using the original update matrices $\textbf{A}_{e,h}$. In addition, the fields derived from previous time steps appear as “equivalent sources” $\bm{\mathcal{J}}$ and $\bm{\mathcal{K}}$ in the update equations for field derivatives.

II-C Equivalent sources for field derivatives

As stated in (15), the equivalent sources $\bm{\mathcal{J}}$ and $\bm{\mathcal{K}}$ for the first-order field derivatives with respect to $\xi_{1}$ are products of $\partial\textbf{A}_{e,h}/\partial\xi_{1}$ and fields E and H. The derivatives of update matrices $\textbf{A}_{e,h}$ are further expressed as:

[TABLE]

The derivative of the matrix D with respect to $\xi_{1}$ is a sparse matrix:

[TABLE]

Since the perturbation in $\xi_{1}$ affects the size of the cells in the $y$-direction only and all entries in $\textbf{D}^{x,z}$ are independent of this perturbation. Hence, $\partial\textbf{D}^{x}/\partial\xi_{1}=\textbf{0}$ and $\partial\textbf{D}^{z}/\partial\xi_{1}=\textbf{0}$. On the other hand, the entities in $\textbf{D}^{y}$ are either constant [math], $1/\Delta y_{p}$ and variable $1/\tilde{\Delta y_{p}}$. Notably, only the derivative of $1/\tilde{\Delta y_{p}}$ with respect to $\xi_{1}$ is not zero:

[TABLE]

The following equation summarizes the derivative of $\textbf{D}^{y}$ with respect to $\xi_{1}$:

[TABLE]

The derivations in (22) and (23) indicate that the derivatives of update matrices $\textbf{A}_{e,h}$ are analytically available. Therefore, the equivalent source exciting the first-order field derivatives can be found with one standard FDTD simulation, exploiting the fields computed by this FDTD simulation and the analytical expression of the derivatives of $\textbf{A}_{e,h}$.

In this 3-D example, the vectors $\bm{\mathcal{J}}$ and $\bm{\mathcal{K}}$ are further written as:

[TABLE]

where $\bm{\mathcal{J}}^{x,y,z}$ and $\bm{\mathcal{K}}^{x,y,z}$ are the vectors of equivalent sources for first-order field derivatives in $x,y,z$-directions respectively.

Fig. 2 illustrates the full-wave computation of the equivalent sources in FDTD. When the fields reach the $P_{\xi_{1}}$ cells coloured red in the grid, the field derivatives are excited. Notably, the derivatives of $\textbf{A}_{e,h}$ and the fields probed at these cells are used to compute the equivalent source for the field derivatives, which further propagate throughout the grid.

III Multi-parametric Sensitivity and Equivalent Sources

III-A Jacobian matrix of field derivatives

Section II proposed a new approach to find the relation of fields and field derivatives in FDTD. In addition, the equivalent source for field derivatives with respect to parameter $\xi$ is located only at the cells of $P_{\xi}$.

This section extends this approach to the multiparametric case. Consider the derivative of the electric field in the $z$-direction sampled at a cell of index $p_{a}$ with respect to an $N$-dimesional vector of design parameters $\bm{\mathrm{\xi}}=\left[\xi_{1},\xi_{2},\cdots\xi_{N}\right]^{T}$:

[TABLE]

which is a $1\times N$ Jacobian matrix.

Assume perturbations $\delta_{1},\delta_{2},...,\delta_{N}$ of $\xi_{1},\xi_{2}...,\xi_{N}$ affect the $x$-component of the cell-size for cells $p_{1},p_{2},...,p_{N}$ respectively. Hence, the vector $\textbf{d}^{x}$ in this 3-D FDTD simulation is written as:

[TABLE]

In addition, an operator $\mathfrak{D}$ is defined to differentiate the entries of FDTD update matrices that are related to the parameters $\bm{\xi}$. In this example $\mathfrak{D}_{\bm{\xi}}$ is:

[TABLE]

Therefore, the derivative of $\textbf{D}^{x}=\text{diag}(\textbf{d}^{x})$ with respect to $\bm{\xi}$ is found as:

[TABLE]

Then, a compact form of the equivalent sources for field derivatives with respect to $\bm{\xi}$ is derived. Consider, for example, the equivalent sources for electric field derivatives $\bm{\mathcal{J}}$:

[TABLE]

where the vector $\bm{\mathcal{J}}^{z}$ contains equivalent sources for the $\textbf{E}^{z}$ derivatives. Substituting (27) into (28) , the vector $\bm{\mathcal{J}}^{z}$ is defined in the time domain as:

[TABLE]

The non-zero entries in $\bm{\mathcal{J}}^{z}$: $\mathcal{J}^{z}_{p_{1}},\mathcal{J}^{z}_{p_{2}},...,\mathcal{J}^{z}_{p_{N}}$ are the equivalent sources for field derivatives with respect to $\xi_{1},\xi_{2},...,\xi_{N}$ respectively. Eq. (29) indicates that the equivalent sources for field derivatives with respect to multiple parameters can be computed simultaneously with the original FDTD simulation.

It was shown before that the derivatives of $\text{E}^{z}$ with respect to $\xi_{1}$ are excited by $\mathcal{J}^{z}_{p_{1}}$ located at cell $p_{1}$; therefore, the field derivative at an observation port $p_{a}$, such as the input or output port of a system, can be found by using the following relation in the frequency domain:

[TABLE]

where $\mathcal{H}_{p_{a}p_{1}}$ is the transfer function from cell $p_{1}$ to port $p_{a}$ and superscript $\hat{}$ denotes a Fourier transform. $\hat{E}^{z}$ and $\hat{\mathcal{J}}^{z}$ are the $z$-component of the Fourier-transformed electric fields and of the equivalent source for electric field derivatives. In general, the transfer functions between port $p_{a}$ and the cells at $p_{1},p_{2},\dots,p_{N}$ are found by exciting the grid with a Gaussian source placed at $p_{a}$ and recording the fields at $p_{a}$ and $p_{n}$ over time :

[TABLE]

where $\mathfrak{F}$ denotes Fourier transform. With these transfer functions and equivalent sources available, $\nabla_{\bm{\mathrm{\xi}}}{\hat{E}^{z}_{p_{a}}}$ is written as:

[TABLE]

[TABLE]

[TABLE]

In (32), the following relation holds for each element of the matrix $\nabla_{\bm{\mathrm{\xi}}}{\hat{E}^{z}}$:

[TABLE]

Eq. (32) indicates that the Jacobian matrix of the field at a cell of interest, which is $p_{a}$ here, can be computed by running one FDTD simulation and exploiting the transfer functions and equivalent sources directly.

In general design problems, where the perturbation of a parameter $\xi_{n}$ is relevant to multiple cells, these affected cells form a subset $P_{\xi_{n}}\subset P$. For example, Fig. 3 illustrates three subsets of the cells affected by the perturbation in the length of the stubs of the example shown in Fig. 1. Therefore, in such cases, the source is distributed over a surface or a volume and (35) is modified as follows:

[TABLE]

III-B Implementation

To numerically implement the proposed method, the FDTD update equations and equivalent sources derived in (14) are discretized. For example, the electric fields in the $z$-direction are updated as:

[TABLE]

where ${A_{e}^{x}}$ and ${A_{e}^{y}}$ are the update coefficients:

[TABLE]

$M_{\varepsilon}$, $d^{x}$, $d^{y}$, $d^{z}$ are three-dimensional matrices represents the permittivity and cell size. In addition, the equivalent sources for the derivatives of the $z$-component of the electric field, with respect to $\xi_{n}$, are updated as:

[TABLE]

If the cell located at $(i,j,k)\in P_{\xi_{n}}$, the derivatives of the corresponding update coefficients are as follows:

[TABLE]

Meanwhile, the derivatives of update coefficients in $y$- and $z$-directions are zero throughout the simulation domain:

[TABLE]

The discretized update equations for magnetic fields and equivalent sources $\mathcal{K}$ share the same structure as of (37) and (39) with $M_{\varepsilon}$ replaced by $M_{\mu}$.

The steps for this method are then presented in the Algorithm 1. This algorithm includes three parts: (1) Mapping perturbations of design parameters to FDTD update coefficients, (2) FDTD simulation and (3) post-processing. The computation overhead for the first and third parts is negligible compared to an FDTD run.

IV High-order field derivatives

Beyond first-order field derivatives, this section further illustrates how high-order derivatives can be computed by applying the approach proposed in Section III. By using the equivalent sources for field derivatives and the transfer functions, field derivatives with respect to a parameter $\xi_{1}$, up to any order can be computed with only two FDTD simulations.

IV-A Equivalent sources for high-order field derivatives

Differentiating (14) with respect to $\xi_{1}$ once more yields the FDTD update equation with respect to second-order field derivatives:

[TABLE]

[TABLE]

Apparently, (42) shares the same structure as that of (14). The second-order field derivatives are updated through the same update matrices $\textbf{A}_{e,h}$. Equivalent sources $\bm{\mathcal{J}}^{(m)}$ and $\bm{\mathcal{K}}^{(m)}$ are used to excite the second-order electric and magnetic field derivatives respectively. The superscript $m$ denotes the equivalent sources for $m$-th order field derivatives. In (43), zero-th and first order field derivatives are used to compute $\bm{\mathcal{J}}^{(2)}$ and $\bm{\mathcal{K}}^{(2)}$. Repeatedly differentiating (42), the equivalent source terms that appear in the $m$-th order version of the previous equations, exciting $m$-th order field derivatives with respect to $\xi_{1}$ are:

[TABLE]

In brief, the equivalent sources for $m$-th order field derivatives depend on up to $(m-1)$-th order field derivatives and the $m$-th order derivatives of the update matrices $\textbf{A}_{e,h}$.

Following the assumption used in (23), the perturbation of $\xi_{1}$ affects $d^{y}_{p},p\in P_{\xi_{1}}$ only. Again, $P_{\xi_{1}}$ is the subset including all the affected cells. Hence, the $m$-th order derivative of the update matrix $\textbf{A}_{e}$ is derived as:

[TABLE]

where

[TABLE]

Using (22), the analytical expression of the non-zero terms in (47) is:

[TABLE]

Hence, only the terms corresponding to the cells of $P_{\xi_{1}}$ are non-zero in $\partial^{m}\textbf{A}_{e,h}/\partial\xi^{m}$. This indicates that the high-order equivalent sources are localized at the $P_{\xi_{1}}$ cells only. For example:

[TABLE]

Since the equivalent sources $\bm{\mathcal{J}}^{(m)}$ are localized at the $P_{\xi_{1}}$ cells, they are dependent on $\dfrac{\partial^{m-1}H^{x,y,z}_{p}}{\partial{\xi_{1}}^{m-1}},p\in P_{\xi_{1}}$, which are the $(m-1)$-th order magnetic field derivatives located at the $P_{\xi_{1}}$ cells. Therefore, the product terms in (45) can be represented as follows:

[TABLE]

where the superscript $\star$ indicates that only the entries corresponding to the cells that would be affected by perturbations of parameters, are non-zero in a matrix. Thus, the computation of the $m$-th order equivalent sources is focused on finding up to $(m-1)$-th order field derivatives at the affected cells only.

Indeed, the $(m-1)$-th order field derivatives are not available directly after running the original FDTD simulation. However, they are excited by $\bm{\mathcal{J}}^{m-1}$ and $\bm{\mathcal{K}}^{m-1}$. Specifically, using the relation stated in (36), ${\partial^{m-1}H^{x,y,z}_{p}}/{\partial\xi^{m-1}}$ is found in the frequency domain as:

[TABLE]

where $\hat{H}$ and $\hat{\mathcal{K}}$ are Fourier-transformed magnetic field and equivalent source; $\mathcal{H}_{pp},p\in P_{\xi_{1}}$ is the transfer function of the cell of $P_{\xi_{1}}$ itself. Assuming that an incident field is placed at $p$, then its self transfer function is found as:

[TABLE]

Similarly, the $(m-1)$-th order electric field derivatives at the $P_{\xi_{1}}$ cells are derived as:

[TABLE]

Therefore, (50) can be further written in the frequency domain as:

[TABLE]

Combining (51), (53), (54) and (45), the equivalent sources up to any order can be computed in the frequency domain as:

[TABLE]

Meanwhile, the equivalent source for high-order magnetic field derivatives is found as:

[TABLE]

Therefore, by recursively exploiting equivalent sources of previous orders, only $\bm{\mathcal{J}}^{(1)}$ and $\bm{\mathcal{K}}^{(1)}$ are needed in order to compute equivalent sources of higher orders $\bm{\mathcal{J}}^{(m)},\bm{\mathcal{K}}^{(m)}\forall m>1$. Here, $\bm{\mathcal{J}}^{(1)}$ and $\bm{\mathcal{K}}^{(1)}$ are referred as the initial states for high-order equivalent sources and are computed on the fly with a single FDTD run.

IV-B High-order field derivatives

Using the equivalent sources for $m$-th order field derivatives with respect to $\xi_{1}$, i.e., $\bm{\mathcal{J}}^{(m)}$ and $\bm{\mathcal{K}}^{(m)}$, the $m$-th order field derivatives at the observation point indexed $p_{a}$ are computed using the relation of (36):

[TABLE]

In (57), the $m$-th order field derivative is found by propagating the equivalent sources from the cells of $P_{\xi_{1}}$ to the observation point through transfer functions between them.

Fig. 4 illustrates the computation of high-order field derivatives by using this proposed method. In addition to transfer functions between observation point and $P_{\xi_{1}}$ cells, the transfer functions of the $P_{\xi_{1}}$ cells themselves are used to compute high-order equivalent sources.

Based on this method, only two FDTD simulations are necessary in order to compute these transfer functions and the initial state of equivalent sources. The first FDTD simulation is run with Gaussian excitation placed at the cells of $P_{\xi_{1}}$ to obtain $\mathcal{H}_{p_{i}p_{i}},p_{i}\in P_{\xi_{1}}$. Afterwards, the second FDTD simulation is run with a Gaussian excitation placed at $p_{a}$. After the grid is excited, the transfer function $\mathcal{H}_{p_{i}p_{a}}$ as well as $\bm{\mathcal{J}}^{(1)}$, $\bm{\mathcal{K}}^{(1)}$ are derived following the procedure stated in Algorithm 1. The equivalent sources for first order field derivatives, $\bm{\mathcal{J}}^{(1)}$ and $\bm{\mathcal{K}}^{(1)}$, are the fundamental terms to compute higher-order equivalent sources.

Algorithm 2 summarizes the procedure, including the two FDTD runs and the post-processing part. The first FDTD run is used to find the transfer functions at the subset of affected cells and the second is used to compute the initial state of equivalent sources as well as the field solution to the original problem. In the post-processing part, high-order equivalent sources are propagated to the observation port through the transfer function.

IV-C Partial mixed derivatives and Hessian matrix

In addition to the Jacobian matrix and high-order field derivatives, the Hessian matrix is another important part of sensitivity analysis and optimization in electromagnetic design. Consider the Hessian matrix of $E^{z}_{p_{a}}$, which is the $z$-component of the electric field at port $p_{a}$, with respect to $N$ parameters $\bm{\xi}=[\xi_{1},\xi_{2},\dots,\xi_{N}]^{T}$:

[TABLE]

The diagonal of this Hessian matrix consists of second-order field derivatives with respect to $\xi_{1},\xi_{2},\dots,\xi_{N}$, respectively. These second-order field derivatives are computed using the method proposed in the previous section. In addition, the upper and lower triangular parts of this matrix are the partial mixed derivatives with respect to any two variables out of $\bm{\xi}$.

Direct differentiation is performed to (37) with respect to $\xi_{u}$ and $\xi_{v}$, $u\neq v$, yielding:

[TABLE]

[TABLE]

Eq. (59) includes the FDTD update equations with respect to partial mixed field derivatives. Evidently, the partial mixed field derivatives are excited by $\bm{\mathcal{J}}^{uv}$ and $\bm{\mathcal{K}}^{uv}$. Assuming that the perturbations of $\xi_{u}$ and $\xi_{v}$ are mapped to the $d^{y}$ of the cells in the subset $P_{\xi_{u}}$ and $P_{\xi_{v}}$ respectively, and $P_{\xi_{u}}\cap P_{\xi_{v}}=\emptyset$. Therefore, the $\frac{\partial^{2}\textbf{A}_{e}}{\partial\xi_{u}\partial\xi_{v}}$ and $\frac{\partial^{2}\textbf{A}_{h}}{\partial\xi_{u}\partial\xi_{v}}$ terms in (60) are assumed to be zero.

As discussed in the previous section, in the derivatives of update matrices, only the terms corresponding to the cells of $P_{\xi_{u}}$ and $P_{\xi_{v}}$ are non-zero. Specifically, in (60):

[TABLE]

where

[TABLE]

Therefore, the product of the field derivative with respect to $\xi_{u}$ and the derivative of the update matrix with respect to $\xi_{v}$ is localized at the cells of $P_{\xi_{u}}$:

[TABLE]

Notably, the first term in (63), which is the derivative of the update matrix with respect to $\xi_{u}$ or $\xi_{v}$ anchors the field derivatives needed at the cells of $P_{\xi_{u}}$ or $P_{\xi_{v}}$, which are then used to compute the equivalent sources. Substituting (63) into (60), $\bm{\mathcal{J}}^{uv}$ is found in frequency domain as:

[TABLE]

where $\bm{\mathcal{K}}^{u}$ and $\bm{\mathcal{K}}^{v}$ are the equivalent sources for the magnetic field derivative with respect to $\xi_{u}$ and $\xi_{v}$ respectively; $\mathcal{H}_{p_{u}p_{v}}(f)$ is the transfer functions between two subsets of affected cells. Similarly,

[TABLE]

Eq. (64) and (65) indicate that the equivalent sources for the partial mixed field derivatives with respect to $\xi_{u}$ and $\xi_{v}$, are localized at the cells of subsets $P_{\xi_{u}}$ and $P_{\xi_{v}}$. For example,

[TABLE]

Using the relation stated in (36), the partial mixed field derivative with respect to $\xi_{u}$ and $\xi_{v}$, at the port $p_{a}$ is found by propagating the equivalent source from the cells of $P_{\xi_{u}}$ and $P_{\xi_{v}}$ to $p_{a}$:

[TABLE]

Algorithm 3 summarizes the computation of the partial-mixed field derivatives using equivalent sources. For $N$ design parameters under consideration, a total of $N+1$ FDTD runs are needed. The $N$ runs are used to compute the transfer functions between the affected subsets corresponding to each parameter. The last FDTD run is used to compute the equivalent sources. Finally, these equivalent sources recorded at different subsets are propagated to the observation point through transfer functions.

V Numerical Results

Two numerical examples, including a dielectric multilayer and a three-stub microstrip line are presented in this section to demonstrate using the proposed methods to compute Jacobian matrix, high-order and partial mixed field derivatives. The results computed by this paper are marked proposed, while results computed by the centered finite-difference method are marked CFD-FDTD.

V-A Dielectric multilayer

The first example is a 1-D dielectric multilayer geometry consisting of three polycarbonate slabs with $d_{1}=0.0051,d_{2}=0.0093,d_{3}=0.0136$ m, shown in Fig. 5. The relative permittivity of each slab is $\varepsilon_{r1}=2.2,\varepsilon_{r2}=3.0,\varepsilon_{r3}=4.0$. Two $0.025$ m air-gaps are located between the slabs. The structure is discretized with $1030$ cells. The cell size $\Delta z=0.424$mm. $4096$ time steps ($\Delta t=0.441$ps) are run and a Gaussian pulse $g(t)=\exp(-(t-t_{0})^{2}/T_{s}^{2})$, with $t_{0}=6Ts$ and $T_{s}=20$ ps is used to excite the grid. The reflection and transmission coefficient $S_{11,12}$ are shown in Fig. 5 as well.

V-A1 Jacobian matrix of the reflection coefficient

The Algorithm 1 proposed in Section III is used to compute the first-order field derivatives with respect to design parameters $d_{1,2,3}$ and $\varepsilon_{r1,r2,r2}$. In this example, the derivatives of the reflected fields observed at the input port of this structure are computed. Hence, the derivatives of the reflection coefficient $S_{11}$ with respect to these parameters can be found by applying the chain rule.

As illustrated in the bottom part of Fig. 5, the input port is at the cell indexed $10$. The three dielectric slabs are distributed over the cells indexed ${(402:414),(509:531),(626:658)}$ respectively.The perturbations of $d_{1},d_{2},d_{3}$ affect the size of the cells indexed ${414,531,658}$, corresponding to the edge of each slab. These cells are coloured grey in Fig. 5. Following the Algorithm 1, the Gaussian excitation is placed at cell ${10}$ to excite the grid. The incident fields and the fields probed at cells ${414,531,658}$ are used to compute the transfer functions between the affected cells and the input port. While updating fields in the main FDTD simulation, the equivalent sources $\mathcal{J}$ and $\mathcal{K}$ are computed in parallel using (28). The equivalent sources for the field derivative with respect to $d_{1,2,3}$ are found by probing $\mathcal{J}$ at cells ${414,531,658}$. That is, at each time step, $\mathcal{J}_{d_{1}}(t)=\mathcal{J}(414)$ $\mathcal{J}_{d_{2}}(t)=\mathcal{J}(531)$ and $\mathcal{J}_{d_{3}}(t)=\mathcal{J}(658)$.

In post-processing, using (36), the field derivatives at the input port with respect to $d_{1}$ are computed by the product of the transfer function between cells ${10}$ and ${414}$, and $\mathcal{J}_{d_{1}}(f)$ in the frequency domain. This post-processing step is repeated to find the field derivatives with respect to other design parameters as well. The derivatives of $S_{11}$ with respect to $d_{1},d_{2},d_{3}$ are reported in Fig. 6 (a). The results computed by using CFD-FDTD are also presented to verify the correctness of this proposed method.

Similarly, the field derivatives with respect to material parameters $\varepsilon_{r1,r2,r2}$ can be computed in parallel. In Fig. 5, the cells affected by the perturbations of $\varepsilon_{r1,r2,r2}$ are filled with slash lines, forming subsets $P_{\varepsilon_{r1}},P_{\varepsilon_{r2}},P_{\varepsilon_{r3}}$. Evidently, the equivalent sources for field derivatives with respect to material parameters, are distributed over each slab. Again, (36) is used to compute the field derivatives as the equivalent sources distribute over a volume. Fig. 6 (b) presents the derivative of $S_{11}$ with respect to $\varepsilon_{r1,r2,r3}$. In summary, Fig. 6 is the Jacobian matrix of the reflection coefficient of this structure over a wide frequency range. Notably, all the information needed to compute this matrix is available with one FDTD simulation, while 12 FDTD simulations are needed to compute this matrix via CFD-FDTD, excluding the runs to ensure the convergence of the result.

V-A2 High-order and partial mixed field derivatives

The high-order derivatives of the reflected fields are also computed using Algorithm 2. For example, the equivalent sources for the second-order field derivatives with respect to $d_{3}$ are derived by using the transfer function of the cell indexed ${658}$ itself, and the equivalent source for the first-order field derivatives, i.e., $\mathcal{J}_{d_{3}}(f)$. Fig. 7 (a) presents the second-order derivatives of $S_{11}$ with respect to $d_{3}$, computed by using the proposed method and CFD-FDTD. Indeed, higher-order derivatives can be found by recursively using the equivalent sources of lower-order.

Furthermore, Fig. 7 (b) shows the third-order derivatives of $S_{11}$ with respect to $d_{2}$. It should be noted that the subtractive cancellation errors associated with CFD constraint its feasibility in finding third or higher-order derivatives. Even in a 1-D FDTD simulation example, it is difficult to find an appropriate step-size in CFD to compute a convergent result. Alternatively, the Multi-Complex Step Derivative (MCSD)-FDTD method proposed in [17] is used to verify the third-order field derivatives in Fig. 7 (b).

In addition to high-order field derivatives, the partial mixed field derivatives with respect to any two design parameters are computed using Algorithm 3. Fig. 7 (c) shows the partial mixed derivative of $S_{11}$ with respect to $d_{1}$ and $d_{3}$, in comparison with the CFD-FDTD result. Two full-wave FDTD simulations are needed to derive this second-order derivative using Algorithm 3. The first FDTD is run to compute the transfer function between the cells affected by the perturbations of $d_{1}$ and of $d_{3}$. The second FDTD is run sampling the equivalent sources for the field derivatives with respect to $d_{1}$ and $d_{3}$. Thus, the equivalent source for this partial mixed field derivatives can be computed by using (60). In contrast, at least $4$ FDTD simulations are needed to compute this derivative via CFD-FDTD.

V-B Cascaded microstrip filter

A microstrip filter with three cascaded stubs, shown in Fig.1 is further studied in this section to demonstrate the proposed methods in three-dimensional FDTD simulations. The structure is discretized with a mesh of $80\times 150\times 16$ cells, with $\Delta x=0.4064$ mm, $\Delta y=0.4233$ mm,and $\Delta z=0.265$ mm. The metallic layout is located above a $0.795$ mm thick substrate with relative permittivity $\varepsilon_{r}=2.2$. The width of microstrips $w=2.4$ mm and the length of the three stubs $\xi_{1}=\xi_{2}=\xi_{3}=12.2$ mm.

The design parameters of interest in this example are the lengths of each stub: $\xi_{1},\xi_{2},\xi_{3}$. A Gaussian pulse $g(t)=\exp(-(t-t_{0})^{2}/T_{s}^{2})$, with $t_{0}=3Ts$ and $T_{s}=15$ ps is used to excite the grid, and $4000$ time steps ($\Delta t=0.441$ps) are run.

V-B1 First-order field derivatives and time domain analysis

In this example, the perturbations of $\xi_{1},\xi_{2},\xi_{3}$ are mapped to the $y$-dimension of the cells of subsets $P_{\xi_{1}},P_{\xi_{2}},P_{\xi_{3}}$ respectively. Here, $P_{\xi_{1}}:(67,40:45,4),P_{\xi_{2}}:(67,85:90,4)$ and $P_{\xi_{3}}:(67,120:125,4)$. These cells are marked red in Fig. 3. To find the derivatives of the reflected fields with respect to each parameter, Gaussian excitation is placed at the input port of this filter. Using (28), the equivalent sources for field derivatives with respect to $\xi_{1},\xi_{2},\xi_{3}$ are computed and recorded at the cells of $P_{\xi_{1}},P_{\xi_{2}},P_{\xi_{3}}$. Fig. 8 shows the time-domain form of this Gaussian excitation and these equivalent sources probed at the cells of $P_{\xi_{1}},P_{\xi_{2}},P_{\xi_{3}}$. By using these equivalent sources and the transfer functions between the input port and the affected cells, the field derivatives at the input port, with respect to $\xi_{1},\xi_{2},\xi_{3}$ can be computed. Fig. 9 shows the derivatives of the reflection coefficient $S_{11}$ of this filter with respect to $\xi_{1}$ and $\xi_{2}$. The results computed by using CFD-FDTD are included to confirm the correctness of this method. Notably, the agreement of the CFD results with those of the proposed method is improved with more time steps. The noise associated in the CFD-FDTD results can be eliminated if the FDTD time-steps are increased from 4000 to 6000.

The time-domain data of the field derivatives computed by the CFD-FDTD method are presented in Fig.10 to explain this phenomenon. Compared to the fields determined by one FDTD, it takes more time steps for the field derivatives with CFD-FDTD method to reach steady state. Meanwhile, the equivalent source probed at $P_{\xi_{1}}$ is computed by the FDTD fields in each time step. Thus, the equivalent sources converge along with the FDTD field solution and can further produce field derivatives without the noise present in CFD. This indicates that, not only does this method compute field derivatives with guaranteed accuracy, it also requires fewer time steps.

V-B2 High-order and partial mixed field derivatives

Finally, the high-order and partial mixed field derivatives of this filter are computed using Algorithms 2 and 3. For example, Fig. 11 (a) shows the second-order derivatives of $S_{11}$ with respect to $\xi_{3}$ compared to the result computed by CFD-FDTD. In addition, Fig. 11 (b) shows the derivatives of $S_{11}$ with respect to $\xi_{2}$ and $\xi_{3}$. These second-order and partial mixed derivatives are the elements of the reflection coefficient Hessian matrix of this structure. Based on the proposed algorithms, a total of 4 FDTD runs are needed to compute the full Hessian matrix, including 3 runs with a Gaussian excitation placed $P_{\xi_{1}}$, $P_{\xi_{2}}$, $P_{\xi_{3}}$ and 1 run with the Gaussian excitation placed at the input port to excite the grid. On the other hand, at least 21 full-wave FDTD runs are needed by using CFD-FDTD. The following table summarizes the number of full-wave FDTD simulation runs with the proposed method, to compute first-order, high-order and partial mixed field derivatives with respect to $N$ design parameters up to $M$-th order.

VI Conclusions

A novel and comprehensive framework to compute field derivatives with respect to multiple design parameters with FDTD up to any order has been presented. The proposed approach utilize standard FDTD simulations, mapping geometrical and material perturbation to equivalent sources. Accuracy-wise, this method is more robust than finite-difference methods, since it is not prone to subtraction errors; moreover, it is more flexible than adjoint variable methods and more scaleable with respect to the number of design variables than complex step methods. Hence, it is a useful addition to the toolbox of time-domain computer-aided analysis and design methods for microwave design.