Computation of High-Order Electromagnetic Field Derivatives with FDTD and the Complex-Step Derivative Approximation

TL;DR

This paper presents a novel method integrating the complex-step derivative approximation into FDTD simulations to accurately compute electromagnetic field derivatives of any order with respect to design parameters, enhancing sensitivity analysis and optimization.

Contribution

The paper introduces a finite-difference free approach using complex-step derivatives within FDTD for precise, high-order electromagnetic field derivatives in a single simulation.

Findings

Accurate derivatives for small parameter perturbations

Enables sensitivity analysis and optimization

Applicable to multi-parametric electromagnetic modeling

Abstract

This paper introduces a new approach for the computation of electromagnetic field derivatives, up to any order, with respect to the material and geometric parameters of a given geometry, in a single Finite-Difference Time-Domain (FDTD) simulation. The proposed method is based on embedding the complex-step derivative (CSD) approximation into the standard FDTD update equations. Being finite-difference free, CSD provides accurate derivative approximations even for very small perturbations of the design parameters, unlike finite-difference approximations that are prone to subtractive cancellation errors. The availability of accurate approximations of field derivatives with respect to design parameters enables studies such as sensitivity analysis of multiple objective functions (as derivatives of those can be derived from field derivatives via the chain rule), uncertainty quantification, as well as multi-parametric modeling and optimization of electromagnetic structures. The theory, FDTD implementation and applications of this technique are presented.

Figures

Tables

Table 1. TABLE I: Operation Count for The Computation of High-order Derivatives (Add./Sub. : Additions/Subtractions; Mult. : Multiplications)

	CFD-FDTD$\star$		MCSD-FDTD		Iterative CSD-FDTD
	Add./Sub.	Mult.	Add./Sub.	Mult.	Add./Sub.	Mult.
$\frac{{\partial }^n\tilde{\text{E}}}{\partial {\xi }^n}$, $1⩽n⩽N$	$2N$		$2^N$	3N	$2N$
$\frac{{\partial }^{m_1+m_2+\cdots +m_N}\tilde{\text{E}}}{{\partial }^{m_1}{\xi }_1{\partial }^{m_2}{\xi }_2\dots {\partial }^{m_N}{\xi }_N}$, $m_i=0,1$	$3^N$		$2^N$	$3^N$	$2N$

Table 2. TABLE II: Comparison on FDTD-based Field Derivative Computation Methods

	MCSD-FDTD	CFD-FDTD	AVM-FDTD
Jacobian Complexity	$O\left(N\right)$	$O\left(N\right)$	$O\left(1\right)$
Hessian Complexity	$O\left(N^2\right)$	$O\left(N^2\right)$	$O\left(N\right)$
Accuracy	Machine accuracy	Subtraction-error, step-size needs testing	Comparable to CFD
Implementation	Simple once complex numbers are defined	Simple; runs standard FDTD	Non-trivial
Advantage	Simple implementation, scalability, accuracy	Simple implementation	Efficient for multiple parameters

Equations

$$\frac{F({\xi }_0+h)-F({\xi }_0-h)}{2h}=\frac{\partial F}{\partial \xi }({\xi }_0)+\frac{h^2}{6}\frac{{\partial }^{3}F}{\partial {\xi }^3}({\xi }_0)+\mathcal{O}(h^4)$$

$$\begin{gathered} F({\xi }_0+jh)=F({\xi }_0)+jh\frac{\partial F}{\partial \xi }({\xi }_0)-\frac{h^2}{2}\frac{{\partial }^{2}F}{\partial {\xi }^2}({\xi }_0) \\ -j\frac{h^3}{6}\frac{{\partial }^{3}F}{\partial {\xi }^3}({\xi }_0)+\mathcal{O}(h^4). \end{gathered}$$

$$\text{Re}\left\{F({\xi }_0+jh)\right\}=F({\xi }_0)-\frac{h^2}{2}\frac{{\partial }^{2}F}{\partial {\xi }^2}({\xi }_0)+\mathcal{O}(h^4)$$

$$\frac{\text{Im}\left\{F({\xi }_0+jh)\right\}}{h}=\frac{\partial F}{\partial \xi }({\xi }_0)-\frac{h^2}{6}\frac{{\partial }^{3}F}{\partial {\xi }^3}({\xi }_0)+\mathcal{O}(h^4)$$

$$y_{\mathrm{in}}=j\tan \frac{2\pi l}{\lambda }\equiv jb$$

$$\frac{\partial b}{\partial l}=\frac{2\pi }{\lambda }\frac{1}{{\cos }^2\frac{2\pi l}{\lambda }}=\frac{2\pi }{\lambda }\frac{1}{{\cos }^2\frac{\pi }{4}}=\frac{4\pi }{\lambda }$$

$$\frac{\partial b}{\partial l}\approx \text{Im}\frac{\tan \frac{2\pi (l+jh)}{\lambda }}{h}$$

$$\frac{\partial b}{\partial l}$$

$${\mathbb{C}}^n:=\{z_1+z_2{j}_n|z_1,z_2\in {\mathbb{C}}^{n-1}\}$$

$$j_1{j}_1$$

$$j_1{j}_2$$

$$z$$

$$=\left(x_1+j_1{y}_1\right)+j_2\left(x_2+j_1{y}_2\right)$$

$$=x_1+j_1{y}_1+j_2{x}_2+j_1{j}_2{y}_2$$

$$\begin{gathered} F\left({\xi }_0+\left(j_1+j_2\right)h\right)=F({\xi }_0)+\left(j_1+j_2\right)h\frac{\partial F}{\partial \xi }({\xi }_0) \\ -2\left(1-j_1{j}_2\right)\frac{h^2}{2}\frac{{\partial }^{2}F}{\partial {\xi }^2}({\xi }_0)-4\left(j_1+j_2\right)\frac{h^3}{6}\frac{{\partial }^{3}F}{\partial {\xi }^3}({\xi }_0) \\ +8\left(1-j_1{j}_2\right)\frac{h^4}{24}\frac{{\partial }^{4}F}{\partial {\xi }^4}({\xi }_0)+\mathcal{O}\left(h^5\right) \end{gathered}$$

$$\begin{gathered} \frac{{\partial }^{2}F}{\partial {\xi }^2}({\xi }_0)=\frac{{\text{Im}}_{12}\left\{F\left({\xi }_0+\left(j_1+j_2\right)h\right)\right\}}{h^2} \\ +\frac{h^2}{3}\frac{{\partial }^{4}F}{\partial {\xi }^4}({\xi }_0)+\mathcal{O}(h^4) \end{gathered}$$

$$\begin{gathered} \frac{{\partial }^{2}F}{\partial {\xi }^2}({\xi }_0)=\frac{F({\xi }_0+h)-2F({\xi }_0)+F({\xi }_0-h)}{h^2} \\ -\frac{h^2}{12}\frac{{\partial }^{4}F}{\partial {\xi }^4}({\xi }_0)+\mathcal{O}(h^4) \end{gathered}$$

$$\frac{{\partial }^{2}b}{\partial l^2}=\frac{8{\pi }^2}{{\lambda }^2}\frac{\sin \frac{2\pi l}{\lambda }}{{\cos }^3\frac{2\pi l}{\lambda }}=\frac{16{\pi }^2}{{\lambda }^2}$$

$$\frac{{\partial }^{2}b}{\partial l^2}\approx {\text{Im}}_{12}\frac{\tan \frac{2\pi (l+(j_1+j_2)h)}{\lambda }}{h^2}$$

$$\frac{{\partial }^{2}b}{\partial l^2}\approx \text{Im}\frac{{\zeta }_2{\zeta }_3-{\zeta }_1{\zeta }_4}{{\zeta }_3^2+{\zeta }_4^2}$$

$$\begin{array}[]{lclcl}\zeta_{1}&=&\sin{\alpha}\cosh^{2}{\beta}&+&\dfrac{j}{2}\cos{\alpha}\sinh{2\beta}\\[8.5359pt] \zeta_{2}&=&\dfrac{1}{2}\cos{\alpha}\sinh{2\beta}&-&j\sin{\alpha}\sinh^{2}{\beta}\\[8.5359pt] \zeta_{3}&=&\cos{\alpha}\cosh^{2}{\beta}&-&\dfrac{j}{2}\sin{\alpha}\sinh{2\beta}\\[8.5359pt] \zeta_{4}&=&-\dfrac{1}{2}\sin{\alpha}\sinh{2\beta}&-&j\cos{\alpha}\sinh^{2}{\beta}\\[8.5359pt] \end{array}$$

$$\frac{{\partial }^{m_1+m_2+\cdots m_N}F}{\partial {\xi }_1^{m_1}\partial {\xi }_2^{m_2}\cdots \partial {\xi }_N^{m_N}}$$

$$\begin{array}[]{l}\dfrac{\partial^{m_{1}+m_{2}+\cdots m_{N}}F}{\partial\xi_{1}^{m_{1}}\partial\xi_{2}^{m_{2}}\cdots\partial\xi_{N}^{m_{N}}}\approx\\[8.5359pt] \dfrac{\text{Im}_{m_{1},\cdots,m_{N}}F\left(\xi_{1}+h_{1}\cdot\displaystyle\sum_{k=1}^{m_{1}}j_{k},\cdots,\xi_{N}+h_{N}\cdot\displaystyle\sum_{k=P-m_{N}+1}^{P}j_{k}\right)}{h_{1}^{m_{1}}h_{2}^{m_{2}}\cdots h_{N}^{m_{N}}}\end{array}$$

$$E_{i^{\prime},m,k}^{x,n+1}=E_{i^{\prime},m,k}^{x,n}+\dfrac{\Delta t}{\varepsilon_{0}\left(\varepsilon_{r}+\left(j_{1}+j_{2}\right)h_{1}\right)}\times\\[11.38092pt] \Bigg{(}\dfrac{H_{i^{\prime},m^{\prime},k}^{z,n^{\prime}}-H_{i^{\prime},m^{\prime}-1,k}^{z,n^{\prime}}}{\Delta y}-\dfrac{H_{i^{\prime},m,k^{\prime}}^{y,n^{\prime}}-H_{i^{\prime},m,k^{\prime}-1}^{y,n^{\prime}}}{\Delta z}\Bigg{)}$$

$$\frac{{\partial }^2{E}_{i^{\prime },m,k}^{x,n}}{\partial {\epsilon }_r^2}\approx \frac{{\text{Im}}_{12}{E}_{i^{\prime },m,k}^{x,n}}{{h_1}^2}$$

$$\begin{array}[]{l}\dfrac{\Delta t}{\varepsilon_{0}\left(\varepsilon_{r}+\left(j_{1}+j_{2}\right)h_{1}\right)}=\\[8.5359pt] \dfrac{\Delta t}{\varepsilon_{0}\varepsilon_{r}\left(1+4h_{1}^{2}\right)}\dfrac{\left(\varepsilon_{r}^{2}+2h_{1}^{2}\right)-(j_{1}+j_{2})h_{1}\varepsilon_{r}-j_{1}j_{2}h_{1}^{2}}{\varepsilon_{r}^{2}+4h_{1}^{2}}\equiv\\[11.38092pt] \alpha^{R}+j_{1}\alpha^{I_{1}}+j_{2}\alpha^{I_{2}}+j_{1}j_{2}\alpha^{I_{12}}\end{array}$$

$$E_{i^{\prime},m,k}^{I_{12},x,n+1}=E_{i^{\prime},m,k}^{I_{12},x,n}\\[14.22636pt] +\alpha^{I_{12}}\Bigg{(}\dfrac{H_{i^{\prime},m^{\prime},k}^{R,z,n^{\prime}}-H_{i^{\prime},m^{\prime}-1,k}^{R,z,n^{\prime}}}{\Delta y}-\dfrac{H_{i^{\prime},m,k^{\prime}}^{R,y,n^{\prime}}-H_{i^{\prime},m,k^{\prime}-1}^{R,y,n^{\prime}}}{\Delta z}\Bigg{)}\\[14.22636pt] +\alpha^{I_{2}}\Bigg{(}\dfrac{H_{i^{\prime},m^{\prime},k}^{I_{1},z,n^{\prime}}-H_{i^{\prime},m^{\prime}-1,k}^{I_{1},z,n^{\prime}}}{\Delta y}-\dfrac{H_{i^{\prime},m,k^{\prime}}^{I_{1},y,n^{\prime}}-H_{i^{\prime},m,k^{\prime}-1}^{I_{1},y,n^{\prime}}}{\Delta z}\Bigg{)}\\[14.22636pt] +\alpha^{I_{1}}\Bigg{(}\dfrac{H_{i^{\prime},m^{\prime},k}^{I_{2},z,n^{\prime}}-H_{i^{\prime},m^{\prime}-1,k}^{I_{2},z,n^{\prime}}}{\Delta y}-\dfrac{H_{i^{\prime},m,k^{\prime}}^{I_{2},y,n^{\prime}}-H_{i^{\prime},m,k^{\prime}-1}^{I_{2},y,n^{\prime}}}{\Delta z}\Bigg{)}\\[14.22636pt] +\alpha^{R}\Bigg{(}\dfrac{H_{i^{\prime},m^{\prime},k}^{I_{12},z,n^{\prime}}-H_{i^{\prime},m^{\prime}-1,k}^{I_{12},z,n^{\prime}}}{\Delta y}-\dfrac{H_{i^{\prime},m,k^{\prime}}^{I_{12},y,n^{\prime}}-H_{i^{\prime},m,k^{\prime}-1}^{I_{12},y,n^{\prime}}}{\Delta z}\Bigg{)}$$

$$\begin{gathered} H_{i^{\prime },m,k^{\prime }}^{y,n+1/2}=H_{i^{\prime },m,k^{\prime }}^{y,n-1/2}+\frac{\Delta t}{{\mu }_0\Delta z}\left(E_{i^{\prime },m,k+1}^{x,n}-E_{i^{\prime },m,k}^{x,n}\right) \\ -\frac{\Delta t}{{\mu }_0\Delta x^{-}}\left(E_{i+1,m,k^{\prime }}^{z,n}-\text{cancelto}0E_{i,m,k^{\prime }}^{z,n}\right) \end{gathered}$$

$$\begin{array}[]{l}\text{Re}\left\{\tilde{\textbf{E}}_{i,m,k}^{n}\>\right\}\approx{\textbf{E}}_{i,m,k}^{n}\\[5.69046pt] \text{Im}_{1}\left\{\tilde{\textbf{E}}_{i,m,k}^{n}\>\right\}/h_{1}\approx\text{Im}_{2}\left\{\tilde{\textbf{E}}_{i,m,k}^{n}\right\}/h_{1}\approx\partial\textbf{E}_{i,m,k}^{n}/\partial\varepsilon_{r}\\[5.69046pt] \text{Im}_{3}\left\{\tilde{\textbf{E}}_{i,m,k}^{n}\>\right\}/h_{2}\approx\text{Im}_{4}\left\{\tilde{\textbf{E}}_{i,m,k}^{n}\>\right\}/h_{2}\approx\partial\textbf{E}_{i,m,k}^{n}/\partial l\\[5.69046pt] \text{Im}_{12}\left\{\tilde{\textbf{E}}_{i,m,k}^{n}\>\right\}/h_{1}^{2}\approx\partial^{2}\textbf{E}_{i,m,k}^{n}/\partial\varepsilon_{r}^{2}\\[5.69046pt] \text{Im}_{13}\left\{\tilde{\textbf{E}}_{i,m,k}^{n}\>\right\}/(h_{1}h_{2})\approx\partial^{2}\textbf{E}_{i,m,k}^{n}/\partial\varepsilon_{r}\partial l\\[5.69046pt] \text{Im}_{34}\left\{\tilde{\textbf{E}}_{i,m,k}^{n}\>\right\}/h_{2}^{2}\approx\partial^{2}\textbf{E}_{i,m,k}^{n}/\partial l^{2}\\[5.69046pt] \text{Im}_{123}\left\{\tilde{\textbf{E}}_{i,m,k}^{n}\>\right\}/(h_{1}^{2}h_{2})\approx\partial^{3}\textbf{E}_{i,m,k}^{n}/\partial\varepsilon_{r}^{2}\partial l\\[5.69046pt] \text{Im}_{134}\left\{\tilde{\textbf{E}}_{i,m,k}^{n}\>\right\}/(h_{1}h_{2}^{2})\approx\partial^{3}\textbf{E}_{i,m,k}^{n}/\partial\varepsilon_{r}\partial l^{2}\\[5.69046pt] \text{Im}_{1234}\left\{\tilde{\textbf{E}}_{i,m,k}^{n}\>\right\}/(h_{1}^{2}h_{2}^{2})\approx\partial^{4}\textbf{E}_{i,m,k}^{n}/\partial\varepsilon_{r}^{2}\partial l^{2}\end{array}$$

$${\epsilon }_r^{i,m,k}=z_1^{i,m,k}+j_2{z}_2^{i,m,k},z_1,z_2\in \mathbb{C}$$

$$z_{1}^{i,m,k}=\left\{\begin{array}[]{ll}\varepsilon_{r}^{\mathrm{substrate}}+j_{1}h,&(i,m,k)\>\text{in substrate}\\[8.5359pt] 1,&(i,m,k)\>\text{in air}\end{array}\right.$$

Full text

I Introduction

Quantifying the influence of material and fabrication tolerances has been a topic of interest since the early days of computer-aided analysis and design of microwave circuits and systems [1, 2, 3]. This topic is even more significant today, with the ever increasing complexity of electromagnetic structures. Nevertheless, research on computational electromagnetic solvers has primarily focused on enhancing their speed and accuracy for well-defined problems, rather than problems defined with some degree of uncertainty and the evaluation of their relevant sensitivities.

Variations in the geometry specifications and the materials of a structure can be considered in a sensitivity analysis [4, 5], by computing derivatives of an output function of interest (in the following referred to as an “objective function”), with respect to the design parameters. A standard approach for this analysis is offered by the finite-difference method. For an input parameter $\xi$, with nominal value $\xi_{0}$, and an objective function $F(\xi)$, the centered finite difference (CFD):

[TABLE]

is a second-order accurate approximation to the derivative of $F$ with respect to $\xi$ (i.e. the leading error term of this approximation is proportional to $h^{2}$ and the notation $\mathcal{O}(h^{4})$ implies that the rest of the terms are proportional to $h^{4}$ or higher powers of $h$). This expression relies on iterative calls of a solver to calculate the objective function at points $\xi_{0}\pm h$ of the parameter $\xi$, without any modification of the solver itself. However, as $h$ decreases, the absolute value of the difference $|F(\xi+h)-F(\xi-h)|$ may become smaller than machine precision (even if $F$ is computed analytically), or the error floor of the numerical method used to determine $F$ (if $F$ is numerically computed, which is the most practically interesting case). Then, the error of (1) diverges and this approximation ceases to be useful. This issue becomes even more important in the context of the Finite-Difference Time-Domain (FDTD) method, where small geometric and material perturbations are desired, to limit associated numerical dispersion errors. Moreover, CFDs require two simulations per parameter, accumulating computational overhead when sensitivities with respect to multiple parameters are considered.

To address this latter issue rather than the former, adjoint variable methods (AVM), originally formulated for sensitivity analysis of electric circuits [6], have been implemented in frequency and time-domain numerical electromagnetic techniques [7]. The key advantage of AVM is that it computes first-order sensitivities of an objective function with respect to multiple parameters, with just one additional simulation of the so-called adjoint problem. To do so though, the adjoint solutions of $N$ perturbed problems with respect to each of the parameters are approximated by the adjoint solution to the unperturbed problem, hence further compromising the accuracy of the method. In terms of FDTD, [8] presented a wave equation based formulation, which has been recently extended to the computation of second-order sensitivities of an objective function [9]. Using the wave equation rather than the first-order system of Maxwell’s equations, which are discretized in conventional FDTD, is a drawback of these methods; a first step towards addressing this was recently presented in [10].

To overcome the subtractive cancellation errors in finite-difference methods, first and higher-order derivative approximations can be computed by considering imaginary step perturbations $jh$ instead of the real steps employed in finite-difference expressions such as (1). Indeed, for a real function $F$, the Taylor expansion is:

[TABLE]

Evidently, (2) is complex-valued, whose real part leads to the approximation:

[TABLE]

and its imaginary part leads to:

[TABLE]

Therefore, the complex-step perturbation of a real function provides a second-order accurate approximation to the unperturbed value of the function as its real part and the derivative of the function with respect to the perturbed parameter as its imaginary part. Notably, (4) is a finite-difference free (hence, free of the associated subtraction errors) approximation called the complex-step derivative approximation (CSD [11]).

It was recently pointed out [12] that (4) can be directly embedded in numerical electromagnetic methods, such as FDTD, to enable the direct calculation of electromagnetic field derivatives by a standard FDTD code. The approach is simple: if the design parameters of interest are perturbed by an imaginary step $jh$, the fields $\bm{\mathrm{E}}^{n}$, $\bm{\mathrm{H}}^{n+1/2}$ produced by the update equations at the $n$-th time step become complex; then $\text{Re}\left\{\bm{\mathrm{E}}^{n}\right\}$, $\text{Re}\left\{\bm{\mathrm{H}}^{n+1/2}\right\}$ are the fields of the unperturbed problem, and $\text{Im}\left\{\bm{\mathrm{E}}^{n}\right\}/h$, $\text{Im}\left\{\bm{\mathrm{H}}^{n+1/2}\right\}/h$ are their derivatives with respect to the parameter that is perturbed by the imaginary step. Moreover, [13] introduced an FDTD-based method to compute partial and high-order derivatives with respect to multiple parameters, using an augmented version of CSD with multiple imaginary dimensions, the multi-complex step derivative (MCSD) approximation [14].

There is ample motivation for further exploring the potential of this approach. First, embedding the CSD and MCSD in existing FDTD codes is straightforward. Second, such a combination allows for the computation of field derivatives on the fly; sensitivities of any field-based objective function can be found by post-processing, simply applying the chain rule. Third, a method that allows for the computation of any field derivative of any order is not just a tool for sensitivity analysis, but also a technique for parametric modeling [15, 16] and uncertainty quantification [17, 18, 19, 20, 21]. Fourth, the high accuracy that complex step methods retain for small perturbations (unlike finite-difference based methods) is particularly important for techniques such as FDTD, where geometric and material perturbations (typically modeled by stretching or squeezing Yee cells) can contribute to numerical dispersion errors.

This paper builds on work reported in [12, 13, 16, 21] to provide a comprehensive presentation of general MCSD approximations implemented in FDTD, to formulate an accurate means of calculating field derivatives, along with the field solution to a problem. We extensively discuss the small, yet necessary modifications of standard FDTD codes to embed MCSD approximations. The computational complexity and overhead of this method, with respect to the standard FDTD, are thoroughly studied. We evaluate the accuracy of the method in a cavity case study with analytically known field derivatives and we compute high-order field derivatives in a 3-D microwave circuit example. Simple as it is, this example illustrates how the subtractive cancellation errors of standard finite-difference methods limit their ability to compute second, third and higher-order derivatives. With high-order derivatives available, parameteric expressions of electromagnetic fields and field-based functions, such as scattering parameters, are derived via Taylor expansions, demonstrating an important application of the proposed method to full-wave analysis based parametric modeling for electromagnetic design.

II The Complex-Step Derivative Approximation

A brief overview of the CSD and MCSD approximations is provided in this section, for the sake of completeness. We do so through specific examples. Further details and notes on the historical origins of this approximation can be found in [11, 14, 22].

II-A Single complex-step derivative approximation

Let us consider the normalized input admittance of an open-ended, lossless transmission-line stub of length $l$:

[TABLE]

At $l/\lambda=0.125$, the normalized susceptance $b=1$. The sensitivity of $b$ with respect to the stub length can be expressed as:

[TABLE]

The CSD approximation to this sensitivity is:

[TABLE]

Forward, backward and centered difference approximations can also be used to find this sensitivity as follows:

[TABLE]

Normalizing $\lambda=1$, based on the analytical expression (6), the errors of all four approximations can be found, as a function of $h$, for $l=0.125$. The results are shown in Fig. 1. These confirm that CFD and CSD are second-order accurate with initially identical error performance, as also shown by the expressions (1) and (4). Yet, the accuracy of CFD is clearly compromised by the subtraction error of finite differences, while the accuracy of CSD is practically as good as that of the analytical solution up to machine precision.

II-B Complex-step derivative approximation of second-order derivatives

To derive a complex-step approximation to a second-order derivative, the bi-complex numbers are introduced. In general, $n$-complex numbers $z\in\mathbb{C}^{n}$ are defined recursively as follows:

[TABLE]

For the bi-complex case ($n=2$), $j_{2}$ is the imaginary unit of a second imaginary dimension that is added to the one of the standard $\mathbb{C}$, for which the imaginary unit $j$ is written as $j_{1}$ in this generalized notation. For these two new imaginary units, and for any additional ones recursively introduced per (11),

[TABLE]

Then, building a bi-complex number from the real and imaginary parts of the standard complex numbers leads to:

[TABLE]

with: $\text{Re}(z)=x_{1}$, $\text{Im}_{1}(z)=y_{1}$, $\text{Im}_{2}(z)=x_{2}$, $\text{Im}_{12}(z)=y_{2}$. Note that there are three types of imaginary parts separately defined here, to distinguish $j_{1}$, $j_{2}$ and $j_{1}j_{2}$ terms. Then, a bi-complex perturbation $\left(j_{1}+j_{2}\right)h$ can lead to a second-order derivative approximation as indicated by re-visiting the Taylor expansion (2):

[TABLE]

Letting $\text{Im}_{12}$ denote the term preceded by $j_{1}j_{2}$, (14) leads to the approximation (in addition to approximations (3), (4)):

[TABLE]

Note that the leading error term is proportional to $h^{2}$. However, as shown in Fig. 1(b), the constant in front of $h^{2}$ is four times larger than the one in the centered finite-difference approximation given below:

[TABLE]

These formulas are applied to approximate the second order derivative of the normalized susceptance $b$ in (5), which is:

[TABLE]

for $l=0.125$, $\lambda=1$. The bi-complex step approximation to this second-order derivative is:

[TABLE]

Expanding the tangent, one can readily re-cast (18) in a conventional $\mathbb{C}$ notation:

[TABLE]

where $\text{Im}\equiv\text{Im}_{1}$, and:

[TABLE]

with $\alpha=2\pi l/\lambda$, $\beta=2\pi h/\lambda$. This approximation is compared to standard forward, backward and centered finite-differences in Fig. 1. Notably, subtraction errors dominate and, indeed, destroy the accuracy of finite-difference approximations, to a greater degree than in the case of first-order derivatives (Fig. 1), with the breakpoint moving to $h=10^{-5}$. On the other hand, the bi-complex step approximation is able to deliver accuracy that is practically equivalent to that of the analytical expression of the second-order derivative, as $h\rightarrow 0$. In general, the accuracy advantage of the complex step approximations increases with the order of the derivative under consideration. With many applications, such as optimization studies and parametric modeling, requiring the computation of second and higher-order derivatives, this observation further motivates the study of complex-step approximations. In the following, their generalized form is presented.

II-C Multi-complex step approximation for arbitrary order derivatives

The last example has paved the way for the introduction of generalized complex-step derivative approximations for an arbitrary number of parameters and derivatives up to any order, introduced as the multi-complex step derivative (MCSD) approximation in [14]. Let us consider a function $F\left(\bm{\mathrm{\xi}}\right)$, $\bm{\mathrm{\xi}}=\left[\xi_{1},\xi_{2},\cdots\xi_{N}\right]^{T}$ and a partial derivative:

[TABLE]

The MCSD approach is to use $m_{1}+m_{2}+\cdots m_{N}\equiv P$ imaginary dimensions, i.e. perform computations in $\mathbb{C}^{P}$, introducing $m_{k}$ imaginary perturbations of each parameter $\xi_{k}$, $k=1,\cdots N$, with magnitude $h_{k}$. Then, the MCSD approximation of (21) becomes:

[TABLE]

In this expression, $\text{Im}_{m_{1},\cdots,m_{N}}$ means that once all perturbations have been applied to $F$, the approximation is actually derived from the $j_{m_{1}}j_{m_{2}}\cdots j_{m_{N}}$ term, just as the second-order derivative approximation (18) was derived from the $j_{1}j_{2}$ term.

The general approximation (22) and its lower-dimension counterparts, presented through the examples of the previous subsections, can be applied to several numerical methods for the analysis, design and optimization of microwave circuits. In the following, we focus on how this approximation can be embedded into the FDTD technique, to compute arbitrary order field derivatives with respect to design parameters along with the full-wave solution to a given problem.

III FDTD-based Computation of Field Derivatives using the Complex-Step Derivative Approximation

III-A Formulation

The general strategy for the implementation of CSD in FDTD is to introduce complex perturbations into the standard FDTD update equations. We explain this through an example of a microstrip segment (Fig. 2), printed on a dielectric substrate. Here, first and second order field derivatives with respect to the length $l$ of the segment and the relative dielectric permittivity $\varepsilon_{r}$ of the substrate are sought for.

III-A1 Field derivatives with respect to material parameters

To compute derivatives with respect to $\varepsilon_{r}$, two imaginary dimensions $j_{1}$ and $j_{2}$ are introduced and $\varepsilon_{r}$ is set to $\varepsilon_{r}+\left(j_{1}+j_{2}\right)h_{1}$ in the FDTD update equations. For example, the update of the $x$-component of the electric field becomes:

[TABLE]

In (23) and the FDTD update equations that follow $\left(i,m,k\right)$ are Yee cell indices, $n$ is a time-step index and $p^{\prime}=p+1/2$ ($p=i,m,k,n$) is the half cell/time-step offset of space/time nodes in the FDTD mesh. Once a class of $n$-complex numbers is available (as discussed in the next subsection), this update equation is implemented in $\mathbb{C}^{2}$ and the MCSD approximation to the derivative $\partial^{2}E_{i^{\prime},m,k}^{x,n}/\partial\varepsilon_{r}^{2}$ is simply:

[TABLE]

To gain further insights into this equation and how field derivatives are deduced from it, one can decompose it into real and (multiple) imaginary parts. To that end, it is recognized that:

[TABLE]

where the super-scripts $R$, $I_{1}$, $I_{2}$, $I_{12}$ indicate the real and multi-imaginary ($j_{1}$, $j_{2}$, $j_{1}j_{2}$) terms, respectively. Then, the update equation for $\text{Im}_{12}E_{i^{\prime},m,k}^{x,n}$ $\equiv$ $E_{i^{\prime},m,k}^{I_{12},x,n}$ becomes:

[TABLE]

Inspection of (26) reveals that (23) is equivalent to standard FDTD-type finite difference equations, which couple the real and imaginary fields. Note also that $|\alpha^{R}|$, $|\alpha^{I_{1}}|$, $|\alpha^{I_{2}}|$, $|\alpha^{I_{12}}|$ $\leq$ $\Delta t/\varepsilon_{0}\varepsilon_{r}$. So, if these finite differences appeared in standard FDTD, they would correspond to media with $\tilde{\varepsilon}_{r}>\varepsilon_{r}$. As a result, these updates are subject to the standard Courant-Friedrichs-Lewy (CFL) stability condition of FDTD. Moreover, the memory requirements of MCSD-FDTD are greater than those of standard FDTD by a factor equal to $2^{n}$, for computations in $\mathbb{C}^{n}$. For $n=2$, for example, the real part of the fields is accompanied by the $I_{1}$, $I_{2}$ and $I_{12}$ terms.

III-A2 Field derivatives with respect to geometric properties

To derive the MCSD expressions for field derivatives with respect to the stub length $l$ in Fig. 2, the corresponding analysis with the finite-difference method is invoked. In standard finite-difference sensitivity analysis, a common approach is to locally stretch or compress the Yee cells by a geometric perturbation $\delta$ of the length of the cells at the edge of the stub [8], as shown in Fig. 3.

Accordingly, the $H_{y}$ node shown in the figure can be updated as:

[TABLE]

where $\Delta x^{-}=\Delta x-\delta$. Likewise, update equations on stretched cells are derived by replacing $\Delta x$ by $\Delta x^{+}=\Delta x+\delta$. Translating this approach into the context of MCSD, the same form of equations is used, with $\Delta x^{+,-}=\Delta x\pm\left(j_{3}+j_{4}\right)h_{2}$, where imaginary dimensions have been introduced through $j_{3,4}$, in addition to those corresponding to the derivatives with respect to $\varepsilon_{r}$. In other words, there is a one-to-one mapping between the perturbations applied to a finite-difference based sensitivity analysis and MCSD.

The relations between the derivatives of the electric field ${\textbf{E}}_{i,m,k}^{n}$ with respect to $l$, $\varepsilon_{r}$ and the solution ${\tilde{\textbf{E}}}_{i,m,k}^{n}$ of the problem with the complex-step perturbations, are summarized here:

[TABLE]

Note that the real part of the field arrays contains the solution to the problem and the imaginary dimensions contain derivatives from first to fourth order, including all entries to the Hessian matrix and all derivatives needed for a second-order Taylor series expansion of the fields with respect to ($l$, $\varepsilon_{r})$. Hence, MCSD-FDTD is more than a sensitivity analysis technique; it is a general framework for field derivative calculation, which can be readily employed in parametric modeling and uncertainty quantification studies.

III-B FDTD Computations in $\mathbb{C}^{n}$

A key step for embedding the MCSD approximation in existing FDTD codes is the implementation of a programming class for multi-complex numbers. With this class at hand, the FDTD arrays are defined as $n$-complex (with $n$ chosen according to the maximum order of the derivatives that are needed) and perturbations of the imaginary part of any design parameter can be introduced, without any other change in the structure of the code.

Since no standard multi-complex number data type is currently available, a custom class, equipped with rules for the initialization of and numerical operations on multi-complex numbers, is introduced in this paper. In particular, an $n$-complex library has been developed for C++, adhering to the C++11 standards [23]. In addition, a MATLAB version of this library has been constructed and will be made available along with this paper.

III-B1 The $n$-complex number class

The class of $n$-complex numbers has been developed under the core concept of class inheritance. In an arbitrary class defining an $n$-complex number, two private members z1 and z2 of the $(n-1)$-complex number class are included, per the recursive definition of (11). For example, a variable of the $bi$-complex data type ($\mathbb{C}^{n}$, $n=2$) consists of two private members z1 and z2 of the $single$-complex number data type (i.e. $z_{1},z_{2}\in\mathbb{C}$).

For computations involving $n$-complex numbers (as those in the FDTD update equations), basic mathematical and linear algebra operators are overloaded, to accept $n$-complex numbers as input arguments. Since the $n$-complex class is defined recursively, these functions can be used repeatedly without duplicate declarations. The above characteristics of the $n$-complex library allow users to define an arbitrary number of imaginary dimensions. In addition to basic operators, it is necessary to have new functions for the extraction of the real and multiple imaginary parts of an $n$-complex variable. For example, functions Im1(z), Im2(z), Im12(z) return $\text{Im}_{1}(z)\equiv y$, $\text{Im}_{2}(z)\equiv z$, $\text{Im}_{12}(z)\equiv\tau$ of a $bi-$complex variable $z=x+j_{1}y+j_{2}z+j_{12}\tau$, respectively.

III-B2 Using the $n$-complex number class in FDTD

The $n$-complex number class is employed to implement MCSD in standard FDTD codes as follows. First, all design parameters of interest are declared as $n$-complex, consisting of two $(n-1)$-complex number arrays $z_{1},z_{2}$, with the same number of cells as that of the computational domain. Let us consider an example where second order field derivatives with respect to the dielectric permittivity of the substrate of a microstrip structure are computed. Then, the array of the relative dielectric permittivities of the ($i,m,k$) Yee cell is defined as:

[TABLE]

where the two arrays $z_{1}$, $z_{2}$ are:

[TABLE]

and

[TABLE]

Since the $n$-complex number class is included in the programming environment, the electric and magnetic field components are transformed to the $n$-complex number data type if they are associated with perturbed cells. This self-defined $n$-complex class is compatible with any existing FDTD codes in both MATLAB and C++ environments. In the following, FDTD codes loaded with the $n$-complex number class will be referred to as MCSD-FDTD. In the next sections, MCSD-FDTD is evaluated in terms of its accuracy and applied to a 3-D simulation scenario.

III-B3 Defining the $n$-complex number class in matrix form

The alternative implementation

III-C Jacobian and Hessian Matrix Computation

This section further studies the cost of computing Jacobian and Hessian matrices via CSD- and MCSD-FDTD.

III-C1 Jacobian matrix

Consider a 3-D FDTD simulation discretized by $I\times M\times K$ cells, and $I+M+K=P$. Then, at each time step, a total of $6P$ field components are computed, where $p=1,2,\cdots P$ is a Yee cell index corresponding to the triad $(i,m,k)$. The first-order derivatives of these field components with respect to $N$ design variables $\bm{\mathrm{\xi}}=\left[\xi_{1},\xi_{2},\cdots\xi_{N}\right]^{T}$ are written in the form of a $P\times N$ Jacobian matrix:

[TABLE]

where

[TABLE]

To derive these matrices, a total of $N$ CSD-FDTD simulations are needed. In each independent CSD-FDTD simulation, one imaginary perturbation is assigned to a specific parameter $\xi_{i}$. In contrast, by using the CFD-FDTD method, at least $2N$ CFD-FDTD simulations are needed.

III-C2 Hessian matrix

In addition to first-order derivatives, the Hessian matrix of a field component can be further derived by using MCSD-FDTD with bi-complex numbers. To obtain the Hessian matrix of $E^{z}_{p}$ with respect to $N$ design parameters $\bm{\mathrm{\xi}}=\left[\xi_{1},\xi_{2},\cdots\xi_{N}\right]^{T}$, written as,

[TABLE]

at least $4N^{2}$ CFD-FDTD simulations are needed. Alternatively, a total of $N^{2}$ MCSD-FDTD simulations can be performed, with two imaginary perturbations $j_{1}h_{1}$, $j_{2}h_{1}$ assigned to $\xi_{i}$, $\xi_{j}$ ($i,j\in[1,N]$) respectively in each independent run.

III-D High-order Partial Mixed Derivatives

In addition to computing the second-order partial mixed field derivatives in the Hessian matrix, MCSD-FDTD is capable of producing $N$-order mixed partial derivatives of field components with respect to $N$ design variables $\bm{\mathrm{\xi}}=\left[\xi_{1},\xi_{2},\cdots\xi_{N}\right]^{T}$:

[TABLE]

$N$ imaginary perturbations $j_{1}h_{1},j_{2}h_{2},\dots,j_{N}h_{N}$ are assigned to each variable respectively. In addition to $N$-th order mixed partial derivative, $N^{\prime}$-th order partial mixed derivatives with respect to $\bm{\mathrm{\xi}}^{\prime}=\left[\xi_{1},\xi_{2},\cdots\xi_{N^{\prime}}\right]^{T}$, where $N^{\prime}=1,\cdots N-1$, are computed at the same time in this MCSD-FDTD simulation. If CFD-FDTD is used to find all these derivatives, at least $3^{N}$ real-valued FDTD simulations are needed.

III-E High-order Derivatives with an Iterative CSD-FDTD Method

The computational redundancy in finding high-order field derivatives with respect to one variable is found by examining (28), in which both $\text{Im}_{1}\left\{\tilde{\textbf{E}}_{i,m,k}^{n}\>\right\}/h_{1}$ and $\text{Im}_{2}\left\{\tilde{\textbf{E}}_{i,m,k}^{n}\>\right\}/h_{2}$ approximate $\partial\textbf{E}_{i,m,k}^{n}/\partial\varepsilon_{r}$. The additional computation overhead becomes heavier if more imaginary perturbations are assigned to one variable. For example, to find $\partial^{3}\textbf{E}_{i,m,k}^{n}/\partial{\varepsilon_{r}}^{3}$, three imaginary perturbations $j_{1}h_{1},j_{2}h_{1},j_{3}h_{1}$are assigned to $\varepsilon_{r}$. The MCSD approximation to the third-order derivative is simply:

[TABLE]

In addition to the third-order derivative, multiple approximations to first and second-order derivatives are redundantly computed:

[TABLE]

To alleviate the computational overhead of MCSD-FDTD, an alternative approach using only one imaginary perturbation per design parameter is proposed. Based on that, derivatives up to any order can be computed iteratively, with a “marching-in-order” approach. This iterative CSD-FDTD method is elucidated through the update equation of the $x$-component of the electric field:

[TABLE]

where

[TABLE]

Direct differentiation of (38) with respect to $\varepsilon_{r}$ yields:

[TABLE]

where

[TABLE]

Introducing an imaginary perturbation $j_{1}h_{1}$ to $\varepsilon_{r}$ in (38) and (40) turns them into complex-valued equations. Then, (40) is updated iteratively by the real and imaginary parts of $\mathcal{J}$ from (38), as:

[TABLE]

Higher order field derivatives can be further derived following the same scheme. In general, the added term ${\mathcal{K}}$ that appears in the $p$-th order field derivative is:

[TABLE]

Therefore, to find the $N-$th order field derivative with respect to a particular variable, $N$ additional terms are needed in ${\mathcal{K}}$, and these terms are available from the update equations of the $(N-1)$-th order field derivatives.

Finally, the operation count for MCSD-FDTD is summarized and compared to that of CFD-FDTD, relative to a single real-valued FDTD simulation, in Table I. Parallelization or multi-threading, both viable options to accelerate these computations, are not considered in this analysis. For one MCSD-FDTD simulation with field arrays in $\mathbb{C}^{n}$, there are $2^{N}$ times the additions/subtractions and $3^{N}$ times the multiplications of a real-valued FDTD simulation.

IV MCSD-FDTD: Validation

IV-A The relation between complex step and accuracy

The accuracy of MCSD-FDTD is first tested in a two-dimensional case study, where field derivatives are analytically available: the derivatives of electromagnetic field components within a rectangular, air-filled, metallic cavity, with respect to the dimensions of the cavity. The cavity is discretized by a uniform mesh of $150\times 100$ Yee cells with $\Delta x=\Delta y=\Delta=1$ mm. Transverse electric (TE) modes with ($H_{x}$, $H_{y}$, $E_{z}$) are considered. The electric field of the ($m$, $n$) mode is:

[TABLE]

with

[TABLE]

In the expressions above, $\alpha=15$ cm is the width and $b=10$ cm is the height of the cavity. Based on these, the derivatives of the ($m,n$) modal fields with respect to $\alpha$ and $b$ are found and compared to their numerically computed values, via MCSD-FDTD, in the time-domain. To that end, the electric field of the ($m,n$) mode is injected as an initial condition (at $t=0$) and its time evolution is simulated.

In this example, first and second order derivatives of the field components with respect to $\alpha$ and $b$ are computed. Therefore, the cells at the edges of the cavity (in both the $x$- and $y$-direction) are perturbed by complex steps $j_{1}h_{1}$ and $j_{2}h_{2}$, and set to $\Delta x^{\prime}\equiv\Delta x\left(1+j_{1}h_{1}\right)$, $\Delta y^{\prime}\equiv\Delta y\left(1+j_{2}h_{2}\right)$, with $h_{1}=h_{2}=10^{-5}$. The corresponding electric field update equation is:

[TABLE]

The results of MCSD-FDTD are compared to the analytical ones in Fig. 4, for the electric field of the ($1,1$) cavity mode. The real part of the electric field in MCSD-FDTD accurately reproduces the electric field of the unperturbed problem (sampled at the center of the cavity), as shown in Fig. 4(a). Moreover, Fig. 4(b) shows the second order derivative of the electric field at the same sampling point, with respect to $\alpha$ and $b$, found by MCSD-FDTD, along with its analytical expression. Note that this derivative grows with time, as:

[TABLE]

The temporal evolution of $\partial^{2}E_{z}/\partial\alpha\partial b$ is accurately represented by the MCSD-FDTD solution, in agreement with its analytical counterpart. To further elucidate the accuracy of the MCSD-based derivative calculation and to compare it to the conventional alternative of centered finite-differences (CFD), the analytical expression of $\partial^{2}E_{z}/\partial\alpha\partial b$ is employed to characterize the accuracy of MCSD and CFD with respect to the step size $h$. The following $l_{\infty}$ error norm is used:

[TABLE]

This registers the maximum error over space and time within 5,000 time steps of the simulation.

In Fig. 5, this relative error norm (48) is plotted with respect to the step size $h$. Notably, the general error trends that were found in the theoretical analysis of CFD and MCSD in section II appear here as well. The important difference is that these errors are now super-imposed to the FDTD dispersion errors. Hence, MCSD-FDTD cannot reach machine accuracy levels. Yet, it does outperform CFD by an order of magnitude, in terms of accuracy. It is also evident that CFD fails to converge, as $h$ is reduced; $h=10^{-4}$ appears to be a breakpoint, where subtraction errors start to dominate. On the other hand, MCSD has a stable performance. As derived in (15) and (16), the leading error term of second-order derivative is proportional to $h^{2}/3$ and $h^{2}/12$ by MCSD and by CFD respectively. Hence, the MCSD-FDTD is prone to higher relative error than CFD-FDTD when $h>10^{-4}$. As MCSD-FDTD produces converged results as step-size is reduced, it is allowed to use $h=10^{-5}$ in our computations in Fig. 4, accurately recovering the unperturbed solution through the real part of the fields.

In terms of execution time, one run of MCSD-FDTD took 16.4 secs. On the other hand, nine FDTD simulations, each taking 2.4 secs are needed to compute the same derivatives as MCSD-FDTD via CFD. Hence, the total execution time was 21.6 secs, excluding runs that would be needed to choose the step $h$ used in the centered finite differences.

IV-B FDTD numerical dispersion and MCSD-FDTD

The relation of the numerical error of field derivatives and FDTD dispersion error is further studied in this section. The Yee cell size, $\Delta$, of the uniform mesh of this cavity is varied from $0.4$ to $2$ mm. In Fig. 5, the relative error of the electric field calculated by FDTD compared to the analytical solution decreases with mesh refinement, along with the FDTD dispersion error. Notably, the relative error of the second-order mixed partial field derivative computed by MCSD-FDTD also decreases quadratically. In contrast, the subtractive cancellation error associated with the CFD method constrains its accuracy even as the FDTD dispersion error is minimized.

Based on Fig. 5 and Fig. 5, the advantage of MCSD-FDTD over CFD-FDTD for the computation of field derivatives is substantial. MCSD-FDTD provides guaranteed field derivatives with second-order accuracy, and it can be improved by reducing the step-size or the Yee cell size in FDTD. Hence, there is no need for additional simulations just to assess the accuracy of the derivative approximation as in CFD-FDTD. Note also that MCSD-FDTD with a Yee cell size $\Delta=2$mm has the same accuracy as CFD-FDTD with cell size $\Delta=1$ mm, in terms of second-order field derivatives.

V The 3-D MCSD-FDTD: Microstrip filter

A microstrip filter geometry, originally studied in [24] and reproduced here in Fig. 6, is now used to demonstrate a three-dimensional application of MCSD-FDTD and its distinct advantages over standard finite-difference methods. In particular, we focus on sensitivities of the $S$-parameters of the filter, with respect to the widths $w_{1,2}$, shown in Fig. 6. To that end, the derivatives of the $S$-parameters, are expressed in terms of the derivatives of the scattered and transmitted fields at the ports of the filter with respect to the design parameters. In turn, these derivatives are found applying MCSD-FDTD. The relevant computations are further explained in the following.

The structure is discretized with a mesh of $80\times 100\times 16$ cells, with $\Delta x=0.4064$ mm, $\Delta y=0.4233$ mm, $\Delta z=0.265$mm. $4000$ time steps ($\Delta t=0.441$ps) are used for the extraction of $S$- parameters in the frequency domain. A Gaussian pulse with half-width $T=15$ ps is used as a source excitation. For the MCSD-FDTD based computation of the sensitivities with respect to $w_{1,2}$, the method of section III-A2 is applied. For $w_{1}$, the $\Delta x$ of the cells with the largest $x$-index modeling the segment of width $w_{1}$ are perturbed by $j_{1}h_{1}\Delta x$. Likewise, the $\Delta y$ of the cells with the largest $y$-index modeling the segment of width $w_{2}$ are perturbed by $j_{2}h_{2}\Delta y$.

In particular, the Yee cell sizes in the $x$- and $y$-directions, $\Delta x$ and $\Delta y$ respectively, are initialized as three-dimensional $80\times 100\times 16$ real-valued arrays. The $\Delta x$ sub-array corresponding to the edge of the first microstrip along the $y$- direction is perturbed by $j_{1}h_{1}$, namely $\Delta x(36,1:46,4)\equiv\Delta x(1+j_{1}h_{1})$. The $\Delta x$ of the neighbouring cells become $\Delta x(37,1:46,4)\equiv\Delta x(1-j_{1}h_{1})$. Similarly, the sub-array $\Delta y$ corresponding to the edge of middle patch in $x$- direction is modified as: $\Delta y(16:66,50,4)\equiv\Delta y(1+j_{2}h_{2})$ and $\Delta y(16:66,51,4)\equiv\Delta y(1-j_{2}h_{2})$, to implement the perturbation in $w_{2}$. No perturbations are needed in the $z-$ direction. Hence, $\Delta z$ is real everywhere. With this choice, MCSD-FDTD can provide the field solution to the unperturbed problem, through the real part of the fields. This is shown in Fig. 6, which demonstrates the excellent agreement in the $S_{11}$ and $S_{21}$ found via standard FDTD applied to the unperturbed problem and MCSD-FDTD.

Furthermore, Fig. 7 shows first and second order derivatives of $S_{21}$, computed via MCSD-FDTD and CFD, for values of $h_{1}$, $h_{2}$ that vary between $10^{-7}$ and $10^{-3}$. These results indicate the significant convergence problem of CFD and the robustness of MCSD. The gap between the two becomes even more significant in the case of $\partial^{2}S_{21}/\partial w_{1}\partial w_{2}$, where the CFD results diverge before reaching any satisfactory level of convergence over the simulated frequency bandwidth. On the other hand, the MCSD-FDTD results remain practically unchanged as $h_{1}$ and $h_{2}$ vary.

With regards to execution time, one run of MCSD-FDTD took 142 seconds for this 3-D microstrip filter simulation. Alternatively, nine FDTD simulations, each taking 13.3 seconds are needed to compute the same derivatives as MCSD-FDTD via CFD. The total execution time is therefore $13.3\times 9=119.7$ seconds, excluding the number of additional runs needed to ensure the convergence of CFD.

The capability of MCSD-FDTD to compute field derivatives up to any order is further presented through the following example. We assign one more imaginary perturbation to the permittivity sub-array corresponding to the substrate of the filter, that is, $\varepsilon_{r}(1:80,1:100,1:3)\equiv\varepsilon_{r}(1+j_{3}h_{3})$. The third-order derivatives of $S$-parameters with respect to widths $w_{1,2}$ and substrate permittivity $\varepsilon_{r}$ are found following the generalized approximation in (22) and shown in Fig. 8 (a). On the other hand, these three imaginary perturbations can be assigned to a particular parameter for the computation of high-order derivatives with respect to one variable. Here, $\Delta z(1:80,1:100,3)\equiv\Delta z(1+j_{1}h_{1}+j_{2}h_{2}+j_{3}h_{3})$ and the neighbouring cells $\Delta x(1:80,1:100,4)\equiv\Delta z(1-j_{1}h_{1}-j_{2}h_{2}-j_{3}h_{3})$. Third-order derivatives of $S$-parameters with respect to substrate thickness $d$ are calculated and shown in Fig. 8 (b). Notably, the computation of this derivative through CFD-FDTD was not possible. The current state-of-art in AVM-FDTD has not gone past second-order derivatives.

Finally, the application of these high-order field derivatives to the parametric modelling of output functions of interest is presented. For example, a parametric model of $S_{21}$ with respect to substrate thickness $d$, can be derived from a Taylor expansion around the nominal value $d_{0}=0.794$ mm:

[TABLE]

Fig.9 shows the $S_{21}$ estimated by a Taylor expansion with thickness $d$ varying from $0.9d$ to $1.1d$, at frequencies of 4.2 GHz and 2.2 GHz. Full-wave simulation results are shown along with the parametric model. It is found that the insertion loss is more sensitive to thickness at 4.2 GHz, which is at the edge of the first pass-band of the filer, and higher order derivatives become important in building an accurate parametric model.

Notably, the effective bandwidth of the model is increased from $4\%$ to $18\%$ at 4.2 GHz by increasing the order of the Taylor expansion. This shows the importance of high-order derivative computations that is enabled by MCSD and the iterative CSD scheme. The wide effective range allows users to estimate the performance of electromagnetic design with varying parameters.

VI Conclusions

This paper presented a general framework for accurate computation of high-order field derivatives with respect to design parameters in FDTD simulations. Indeed, the accuracy and robustness of the MCSD approximation and the relative simplicity of its implementation in FDTD (independently of the order of derivative needed), just by including an appropriately defined class of multi-complex numbers, are the most significant advantages of this method. In a nutshell, MCSD-FDTD offers the versatility of finite-differences, as it allows for the computation of field derivatives of any order with FDTD, while its accuracy significantly surpasses that of finite-difference methods, especially for small perturbation steps.

A detailed analysis of the computational cost of CSD and MCSD was presented. Based on this analysis, first-order derivatives with respect to multiple design parameters can be more efficiently computed by applying CSD-FDTD to each parameter rather than a full MCSD-FDTD analysis. For high-order derivatives, CSD-FDTD can be iteratively applied in a marching-in-order scheme to alleviate the computational overhead of MCSD-FDTD. This is being said, MCSD-FDTD does provide a complete framework to compute field derivatives of any order. To summarize, a comparison between the proposed method and two popular alternatives, CFD and AVM respectively, is presented in Table II. Notably, AVM has so far been formulated towards the computation of derivatives of a given, field-based objective function rather than field derivatives themselves.

While this paper focused on FDTD, the underlying theory of complex step derivative approximations can clearly be combined with other full-wave and equivalent circuit analysis techniques (in frequency and time-domain), employed in electromagnetic design.