From Time-Collocated to Leapfrog Fundamental Schemes for ADI and CDI FDTD Methods

Abstract

The leapfrog schemes have been developed for unconditionally stable alternating-direction implicit (ADI) finite-difference time-domain (FDTD) method, and recently the complying-divergence implicit (CDI) FDTD method. In this paper, the formulations from time-collocated to leapfrog fundamental schemes are presented for ADI and CDI FDTD methods. For the ADI FDTD method, the time-collocated fundamental schemes are implemented using implicit E-E and E-H update procedures, which comprise simple and concise right-hand sides (RHS) in their update equations. From the fundamental implicit E-H scheme, the leapfrog ADI FDTD method is formulated in conventional form, whose RHS are simplified into the leapfrog fundamental scheme with reduced operations and improved efficiency. For the CDI FDTD method, the time-collocated fundamental scheme is presented based on locally one-dimensional (LOD) FDTD method with complying divergence. The formulations from time-collocated to leapfrog schemes are provided, which result in the leapfrog fundamental scheme for CDI FDTD method. Based on their fundamental forms, further insights are given into the relations of leapfrog fundamental schemes for ADI and CDI FDTD methods. The time-collocated fundamental schemes require considerably fewer operations than all conventional ADI, LOD and leapfrog ADI FDTD methods, while the leapfrog fundamental schemes for ADI and CDI FDTD methods constitute the most efficient implicit FDTD schemes to date.

1. Introduction

One of the most popular methods in computational electromagnetics is the finite-difference time-domain (FDTD) method [1,2]. It is a leapfrog explicit scheme that involves time-staggered electric (E) and magnetic (H) fields update with time step size limited by the Courant–Friedrichs–Lewy (CFL) stability constraint. To overcome the CFL constraint, unconditionally stable FDTD methods have been developed. They include alternating direction implicit (ADI) FDTD [3,4,5,6], split-step (SS) and locally one-dimensional (LOD) FDTD [7,8,9,10,11,12,13,14,15,16], Crank–Nicolson (CN)-based approximate FDTD methods [17,18,19,20,21], etc. (Note that [17] has pointed out a potential instability in [22].) Many of the unconditionally stable methods have been adapted directly from the classical implicit schemes for parabolic, elliptic and hyperbolic equations [23,24,25,26,27,28,29]. In their conventional update procedures involving time-collocated fields, they have matrix operators at not only the left-hand sides (LHS) (making them implicit), but also the right-hand sides (RHS). Despite the unconditional stability, these update procedures often lead to complicated equations with substantial floating-point operations (flops) that would lower the overall efficiency gain. To reduce the flops and improve the efficiency, we have developed the fundamental implicit schemes that feature simple, concise and matrix-operator-free RHS [30]. The fundamental schemes enable efficient implementations for ADI, LOD, SS, CN and many classical implicit schemes. Moreover, they provide a unified formulation for simplification and interlinking of the conventional implicit FDTD schemes through physical and auxiliary field vectors [31].

Thus far, most conventional and fundamental implicit schemes involve time-collocated electric and magnetic fields at full- and half-integer (or arbitrary intermediate) time-indices. Alternatively, there exists a leapfrog scheme for the ADI FDTD method that involves time-staggered electric and magnetic fields [32,33,34,35]. Further developments of the leapfrog scheme have been carried out including several recent works [36,37,38,39,40,41,42,43,44,45,46]. Although the leapfrog ADI FDTD method is (slightly) more efficient than the conventional time-collocated one, the RHS of its update procedures are still complicated for comprising some products of matrix operator partitions. Unlike the leapfrog explicit FDTD method, analyses of the divergence properties show that the leapfrog ADI FDTD method does not have complying divergence and violates Gauss’s law just like many other implicit schemes [47]. To address the divergence issue, the divergence-preserving ADI FDTD method [48] and the fundamental LOD FDTD method with complying divergence have been developed [11]. While satisfying Gauss’s law, the latter method exploits the efficient fundamental scheme and finds useful applications. Meanwhile, both complying-divergence implicit (CDI) FDTD methods involve time-collocated electric and magnetic fields. Recently, a leapfrog scheme for the CDI FDTD method has been introduced that involves time-staggered fields [49]. It has implicit update procedures in the form of fundamental scheme with RHS free of operators, and explicit update procedures that are compatible with those of leapfrog explicit FDTD method. The leapfrog CDI FDTD method has featured many advantages including unconditional stability, complying divergence, simplicity and efficient leapfrog update procedures.

While the fundamental schemes (with matrix-operator-free RHS) have been useful for simplifications of many implicit FDTD methods, the leapfrog schemes (with time-staggered fields) have been achieved only for ADI and CDI FDTD methods so far. Moreover, most fundamental schemes have been applied for time-collocated fields, so their relations and transformations to leapfrog schemes remain unclear. In this paper, we present the formulations from time-collocated to leapfrog fundamental schemes for ADI and CDI FDTD methods. In Section 2, for the ADI FDTD method, the time-collocated fundamental schemes are implemented using implicit E-E and E-H update procedures, with the latter involving magnetic field in the second implicit update procedure. In these implementations, the update equations comprise simple and concise RHS that contain only the intrinsic matrix operator partitions, but without their complicated products that exist in the conventional form. From the fundamental implicit E-H scheme, the leapfrog ADI FDTD method is formulated in conventional form, whose RHS still involve considerable operations due to the products of matrix operator partitions. Using auxiliary variables, the RHS are simplified resulting in the leapfrog fundamental scheme for ADI FDTD method with reduced operations and improved efficiency. In Section 3, for the CDI FDTD method, the time-collocated fundamental LOD FDTD method with complying divergence is presented. It turns out that the method may be more aptly regarded as the time-collocated fundamental scheme for CDI FDTD method, in view of the close connection to the recent leapfrog CDI FDTD method. In particular, the formulations from time-collocated to leapfrog schemes are provided, which result in the leapfrog fundamental scheme for CDI FDTD method.

2. From Time-Collocated to Leapfrog Fundamental Schemes for ADI FDTD Method

2.1. Time-Collocated Fundamental Schemes for ADI FDTD Method

The conventional ADI FDTD method is a second-order implicit FDTD scheme and comprises two update procedures as [5]

$$\left(\mathbf{I}{}_6-\frac{\Delta t}{2}\mathbf{A}\right)\mathbf{u}{}_{\mathrm{a}}^{n+\frac{1}{2}}=\left(\mathbf{I}{}_6+\frac{\Delta t}{2}\mathbf{B}\right)\mathbf{u}{}_{\mathrm{a}}^n$$

(1)

$$\left(\mathbf{I}{}_6-\frac{\Delta t}{2}\mathbf{B}\right)\mathbf{u}{}_{\mathrm{a}}^{n+1}=\left(\mathbf{I}{}_6+\frac{\Delta t}{2}\mathbf{A}\right)\mathbf{u}{}_{\mathrm{a}}^{n+\frac{1}{2}}.$$

(2)

Here, $\mathbf{u}{}_{\mathrm{a}}$’s denote the physical field vectors that are time-collocated at full- ($n,n+1$) and half-integer ($n+\frac{1}{2}$) time indices: $\mathbf{u}{}_{\mathrm{a}}={\left[\mathbf{E}{}_{\mathrm{a}}\mathbf{H}{}_{\mathrm{a}}\right]}^T={\left[E_{\mathrm{a}x},E_{\mathrm{a}y},E_{\mathrm{a}z},H_{\mathrm{a}x},H_{\mathrm{a}y},H_{\mathrm{a}z}\right]}^T$. The subscripts ‘a’ signify the quantities for ADI FDTD that may be different from the subsequent ones (for CDI FDTD or others). $\mathbf{I}{}_6$ is the $6\times 6$ identity matrix, $\mathbf{A}$ and $\mathbf{B}$ are the matrix operators given specifically by

$$\begin{array}{c}\mathbf{A}=\left[\begin{array}{cc}\mathbf{0}& \mathbf{A}{}_{12}\\ \mathbf{A}{}_{21}& \mathbf{0}\end{array}\right],\mathbf{B}=\left[\begin{array}{cc}\mathbf{0}& \mathbf{B}{}_{12}\\ \mathbf{B}{}_{21}& \mathbf{0}\end{array}\right]\end{array}$$

(3)

$$\begin{array}{c}\mathbf{A}{}_{12}=\left[\begin{array}{ccc}0& 0& \frac{1}{ϵ}\frac{\partial }{\partial y}\\ \frac{1}{ϵ}\frac{\partial }{\partial z}& 0& 0\\ 0& \frac{1}{ϵ}\frac{\partial }{\partial x}& 0\end{array}\right],\mathbf{A}{}_{21}=\left[\begin{array}{ccc}0& \frac{1}{\mu }\frac{\partial }{\partial z}& 0\\ 0& 0& \frac{1}{\mu }\frac{\partial }{\partial x}\\ \frac{1}{\mu }\frac{\partial }{\partial y}& 0& 0\end{array}\right]\end{array}$$

(4)

$$\begin{array}{c}\mathbf{B}{}_{12}=\left[\begin{array}{ccc}0& \frac{-1}{ϵ}\frac{\partial }{\partial z}& 0\\ 0& 0& \frac{-1}{ϵ}\frac{\partial }{\partial x}\\ \frac{-1}{ϵ}\frac{\partial }{\partial y}& 0& 0\end{array}\right],\mathbf{B}{}_{21}=\left[\begin{array}{ccc}0& 0& \frac{-1}{\mu }\frac{\partial }{\partial y}\\ \frac{-1}{\mu }\frac{\partial }{\partial z}& 0& 0\\ 0& \frac{-1}{\mu }\frac{\partial }{\partial x}& 0\end{array}\right]\end{array}$$

(5)

($\mathbf{0}$ is the $3\times 3$ zero matrix). Note that (1) and (2) are implicit update equations due to the presence of LHS matrix operators that are no longer identity (or diagonal). Moreover, their RHS also have matrix operators that would result in substantial operations as well. This can be seen from the implementation of (1) and (2) using their update equations in terms of matrix operator partitions (4) and (5) as

$$\left(\mathbf{I}{}_3-\frac{\Delta t^2}{4}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\right)\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}=\mathbf{E}{}_{\mathrm{a}}^n+\frac{\Delta t^2}{4}\mathbf{A}{}_{12}\mathbf{B}{}_{21}\mathbf{E}{}_{\mathrm{a}}^n+\frac{\Delta t}{2}\left(\mathbf{A}{}_{12}+\mathbf{B}{}_{12}\right)\mathbf{H}{}_{\mathrm{a}}^n$$

(6)

$$\mathbf{H}{}_{\mathrm{a}}^{n+\frac{1}{2}}=\mathbf{H}{}_{\mathrm{a}}^n+\frac{\Delta t}{2}\mathbf{A}{}_{21}\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}+\frac{\Delta t}{2}\mathbf{B}{}_{21}\mathbf{E}{}_{\mathrm{a}}^n$$

(7)

$$\left(\mathbf{I}{}_3-\frac{\Delta t^2}{4}\mathbf{B}{}_{12}\mathbf{B}{}_{21}\right)\mathbf{E}{}_{\mathrm{a}}^{n+1}=\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}+\frac{\Delta t^2}{4}\mathbf{B}{}_{12}\mathbf{A}{}_{21}\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}+\frac{\Delta t}{2}\left(\mathbf{A}{}_{12}+\mathbf{B}{}_{12}\right)\mathbf{H}{}_{\mathrm{a}}^{n+\frac{1}{2}}$$

(8)

$$\mathbf{H}{}_{\mathrm{a}}^{n+1}=\mathbf{H}{}_{\mathrm{a}}^{n+\frac{1}{2}}+\frac{\Delta t}{2}\mathbf{A}{}_{21}\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}+\frac{\Delta t}{2}\mathbf{B}{}_{21}\mathbf{E}{}_{\mathrm{a}}^{n+1}$$

(9)

where $\mathbf{I}{}_3$ is the $3\times 3$ identity matrix. As is evident from the RHS of (6)–(9), there are multiple sums and complicated products of the matrix operator partitions, e.g., $\mathbf{A}{}_{12}\mathbf{B}{}_{21}$ & $\mathbf{B}{}_{12}\mathbf{A}{}_{21}$, etc. These would make the calculations more cumbersome and increase the flops count considerably.

2.1.1. Fundamental Implicit E-E Scheme for ADI FDTD Method

In the conventional ADI FDTD method, all terms and operations at the right-hand sides are calculated as they are, without identifying nor omitting any possible repeated ones in the update equations. To reduce the flops and improve the efficiency, one can introduce some auxiliary variables whose judicious exploitation would bypass the repeated terms and operations in the intermediate steps without incurring extra memory space. To that end, we have developed the fundamental ADI FDTD method, or in short, FADI FDTD method, which comprises the following update procedures:

$$\mathbf{v}{}^n=\mathbf{u}{}_{\mathrm{a}}^n-\mathbf{v}{}^{n-\frac{1}{2}}$$

(10)

$$\left(\frac{1}{2}\mathbf{I}{}_6-\frac{\Delta t}{4}\mathbf{A}\right)\mathbf{u}{}_{\mathrm{a}}^{n+\frac{1}{2}}=\mathbf{v}{}^n$$

(11)

$$\mathbf{v}{}^{n+\frac{1}{2}}=\mathbf{u}{}_{\mathrm{a}}^{n+\frac{1}{2}}-\mathbf{v}{}^n$$

(12)

$$\left(\frac{1}{2}\mathbf{I}{}_6-\frac{\Delta t}{4}\mathbf{B}\right)\mathbf{u}{}_{\mathrm{a}}^{n+1}=\mathbf{v}{}^{n+\frac{1}{2}}.$$

(13)

Through the auxiliary field vectors $\mathbf{v}$’s at full- and half-integer time indices: $\mathbf{v}={\left[\mathbf{e}\mathbf{h}\right]}^T$, the RHS of (10)–(13) become free of any matrix operator and contain only vectors. If there exist nonzero initial fields, the auxiliary field vector is initialized with

$$\begin{array}{c}\mathbf{v}{}^{-\frac{1}{2}}=\left(\frac{1}{2}\mathbf{I}{}_6-\frac{\Delta t}{4}\mathbf{B}\right)\mathbf{u}{}_{\mathrm{a}}^0.\end{array}$$

(14)

Upon substituting the matrix operator partitions (4) and (5), the implementation of (10)–(13) can be carried out using the update equations as follows [50]:

$$\mathbf{e}{}^n=\mathbf{E}{}_{\mathrm{a}}^n-\mathbf{e}{}^{n-\frac{1}{2}}$$

(15)

$$\left(\frac{1}{2}\mathbf{I}{}_3-\frac{\Delta t^2}{8}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\right)\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}=\mathbf{e}{}^n+\frac{\Delta t}{2}\mathbf{A}{}_{12}\mathbf{h}{}^n$$

(16)

$$\mathbf{h}{}^{n+\frac{1}{2}}=\mathbf{h}{}^n+\frac{\Delta t}{2}\mathbf{A}{}_{21}\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}$$

(17)

$$\mathbf{e}{}^{n+\frac{1}{2}}=\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}-\mathbf{e}{}^n$$

(18)

$$\left(\frac{1}{2}\mathbf{I}{}_3-\frac{\Delta t^2}{8}\mathbf{B}{}_{12}\mathbf{B}{}_{21}\right)\mathbf{E}{}_{\mathrm{a}}^{n+1}=\mathbf{e}{}^{n+\frac{1}{2}}+\frac{\Delta t}{2}\mathbf{B}{}_{12}\mathbf{h}{}^{n+\frac{1}{2}}$$

(19)

$$\mathbf{h}{}^{n+1}=\mathbf{h}{}^{n+\frac{1}{2}}+\frac{\Delta t}{2}\mathbf{B}{}_{21}\mathbf{E}{}_{\mathrm{a}}^{n+1}.$$

(20)

Notice that the RHS of (15)–(20) become simpler, concise and contain only each of the intrinsic matrix operator partitions, i.e., $\mathbf{A}{}_{12}$, $\mathbf{A}{}_{21}$, $\mathbf{B}{}_{12}$, $\mathbf{B}{}_{21}$. There is no more sum nor product of the matrix operator partitions like before in (6)–(9), e.g., no more $\mathbf{A}{}_{12}\mathbf{B}{}_{21}$ & $\mathbf{B}{}_{12}\mathbf{A}{}_{21}$, etc. Furthermore, the main implicit update procedures correspond to the physical electric fields $\mathbf{E}{}_{\mathrm{a}}$’s in (16) and (19). Therefore, (15)–(20) can be aptly regarded as the fundamental implicit E-E scheme for the ADI FDTD method, where ‘E-E’ signifies the electric fields in both first and second implicit update procedures, cf. (16) and (19).

To recover the physical magnetic fields $\mathbf{H}{}_{\mathrm{a}}$’s, one can simply sum the auxiliary ones as

$$\begin{array}{c}\hfill \mathbf{H}{}_{\mathrm{a}}^{n+\frac{1}{2}}=\mathbf{h}{}^n+\mathbf{h}{}^{n+\frac{1}{2}},\begin{array}{c}\hfill \mathbf{H}{}_{\mathrm{a}}^{n+1}=\mathbf{h}{}^{n+\frac{1}{2}}+\mathbf{h}{}^{n+1}.\end{array}\end{array}$$

(21)

In most cases, such output of $\mathbf{H}{}_{\mathrm{a}}$ is not always required but only when necessary after certain simulation duration, since having $\mathbf{E}{}_{\mathrm{a}}$ is rather sufficient to determine many parameters. Thus, the main iterations of (15)–(20) can be executed continually and efficiently without the output (21) except when needed.

2.1.2. Fundamental Implicit E-H Scheme for ADI FDTD Method

Alternative to the above implementation, one can carry out the implicit procedures using the update equations firstly for physical electric field and secondly for physical magnetic field as follows [51]:

$$\mathbf{e}{}^n=\mathbf{e}{}^{n-\frac{1}{2}}+\frac{\Delta t}{2}\mathbf{B}{}_{12}\mathbf{H}{}_{\mathrm{a}}^n$$

(22)

$$\left(\frac{1}{2}\mathbf{I}{}_3-\frac{\Delta t^2}{8}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\right)\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}=\mathbf{e}{}^n+\frac{\Delta t}{2}\mathbf{A}{}_{12}\mathbf{h}{}^n$$

(23)

$$\mathbf{e}{}^{n+\frac{1}{2}}=\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}-\mathbf{e}{}^n$$

(24)

$$\mathbf{h}{}^{n+\frac{1}{2}}=\mathbf{h}{}^n+\frac{\Delta t}{2}\mathbf{A}{}_{21}\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}$$

(25)

$$\left(\frac{1}{2}\mathbf{I}{}_3-\frac{\Delta t^2}{8}\mathbf{B}{}_{21}\mathbf{B}{}_{12}\right)\mathbf{H}{}_{\mathrm{a}}^{n+1}=\mathbf{h}{}^{n+\frac{1}{2}}+\frac{\Delta t}{2}\mathbf{B}{}_{21}\mathbf{e}{}^{n+\frac{1}{2}}$$

(26)

$$\mathbf{h}{}^{n+1}=\mathbf{H}{}_{\mathrm{a}}^{n+1}-\mathbf{h}{}^{n+\frac{1}{2}}.$$

(27)

Again like (15)–(20), the RHS of (22)–(27) are simple, concise and contain only each of the intrinsic matrix operator partitions without their sum nor product. However, unlike (16) and (19) previously, the main implicit update procedures (23) and (26) correspond to $\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}$ and $\mathbf{H}{}_{\mathrm{a}}^{n+1}$, respectively. Therefore, (22)–(27) can be aptly regarded as the fundamental implicit E-H scheme for ADI FDTD method, where ‘E-H’ signifies the electric and magnetic fields in the first and second implicit update procedures, respectively, cf. (23) and (26).

To recover the remaining physical electric and magnetic fields, one can sum the auxiliary ones when necessary as

$$\begin{array}{c}\hfill \mathbf{H}{}_{\mathrm{a}}^{n+\frac{1}{2}}=\mathbf{h}{}^n+\mathbf{h}{}^{n+\frac{1}{2}},\mathbf{E}{}_{\mathrm{a}}^{n+1}=\mathbf{e}{}^{n+\frac{1}{2}}+\mathbf{e}{}^{n+1}.\end{array}$$

(28)

For nonzero initial fields, the time-staggered auxiliary fields should be initialized with

$$\begin{array}{c}\left[\begin{array}{cc}\mathbf{I}{}_3& \frac{\Delta t}{2}\mathbf{A}{}_{12}\\ -\frac{\Delta t}{2}\mathbf{B}{}_{21}& \mathbf{I}{}_3\end{array}\right]\left[\begin{array}{c}\mathbf{e}{}^{-\frac{1}{2}}\\ \mathbf{h}{}^0\end{array}\right]=\left[\begin{array}{cc}\frac{1}{2}\mathbf{I}{}_3-\frac{\Delta t^2}{8}\mathbf{A}{}_{12}\mathbf{A}{}_{21}& \frac{\Delta t}{2}\mathbf{A}{}_{12}\\ \mathbf{0}& \frac{1}{2}\mathbf{I}{}_3+\frac{\Delta t^2}{8}\mathbf{B}{}_{21}\mathbf{B}{}_{12}\end{array}\right]\left[\begin{array}{c}\mathbf{E}{}_{\mathrm{a}}^{-\frac{1}{2}}\\ \mathbf{H}{}_{\mathrm{a}}^0\end{array}\right].\end{array}$$

(29)

This is expressed in terms of $\mathbf{E}{}_{\mathrm{a}}^{-\frac{1}{2}}$ and $\mathbf{H}{}_{\mathrm{a}}^0$ for the time-staggered initial physical fields in place of (14) for the time-collocated ones. Besides the fundamental implicit E-H and E-E schemes presented in the current and previous subsections, other alternative implicit procedures are also possible, such as fundamental implicit H-E and H-H schemes, etc.

2.2. From Time-Collocated to Leapfrog Schemes for ADI FDTD Method

While the leapfrog ADI FDTD method has been derived before from the conventional ADI FDTD method, it is not clear how the derivation should proceed from the fundamental scheme that involves intermediate auxiliary fields. Here, we shall formulate the leapfrog scheme from the fundamental implicit E-H scheme for ADI FDTD method, since both schemes involve time-staggered physical electric and magnetic fields. To that end, substituting (24) at one time step backward into (22), as well as (25) into (27) at one time step backward yields the auxiliary electric and magnetic fields at n as

$$\begin{array}{c}\mathbf{e}{}^n=\mathbf{E}{}_{\mathrm{a}}^{n-\frac{1}{2}}-\mathbf{e}{}^{n-1}+\frac{\Delta t}{2}\mathbf{B}{}_{12}\mathbf{H}{}_{\mathrm{a}}^n\end{array}$$

(30)

$$\begin{array}{c}\mathbf{h}{}^n=\mathbf{H}{}_{\mathrm{a}}^n-\mathbf{h}{}^{n-1}-\frac{\Delta t}{2}\mathbf{A}{}_{21}\mathbf{E}{}_{\mathrm{a}}^{n-\frac{1}{2}}.\end{array}$$

(31)

Upon using (30) and (31), (23) can be manipulated as follows:

$$\begin{array}{cc}\hfill (\frac{1}{2}\mathbf{I}{}_3& -\frac{\Delta t^2}{8}\mathbf{A}{}_{12}\mathbf{A}{}_{21})\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & =\mathbf{e}{}^n+\frac{\Delta t}{2}\mathbf{A}{}_{12}\mathbf{h}{}^n\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill & =\mathbf{E}{}_{\mathrm{a}}^{n-\frac{1}{2}}-\mathbf{e}{}^{n-1}+\frac{\Delta t}{2}\mathbf{B}{}_{12}\mathbf{H}{}_{\mathrm{a}}^n+\frac{\Delta t}{2}\mathbf{A}{}_{12}\left(\mathbf{H}{}_{\mathrm{a}}^n-\mathbf{h}{}^{n-1}-\frac{\Delta t}{2}\mathbf{A}{}_{21}\mathbf{E}{}_{\mathrm{a}}^{n-\frac{1}{2}}\right)\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill & =\mathbf{E}{}_{\mathrm{a}}^{n-\frac{1}{2}}-\frac{\Delta t^2}{4}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\mathbf{E}{}_{\mathrm{a}}^{n-\frac{1}{2}}-\left(\mathbf{e}{}^{n-1}+\frac{\Delta t}{2}\mathbf{A}{}_{12}\mathbf{h}{}^{n-1}\right)+\frac{\Delta t}{2}\left(\mathbf{A}{}_{12}+\mathbf{B}{}_{12}\right)\mathbf{H}{}_{\mathrm{a}}^n\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill & =\left(\frac{1}{2}\mathbf{I}{}_3-\frac{\Delta t^2}{8}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\right)\mathbf{E}{}_{\mathrm{a}}^{n-\frac{1}{2}}+\frac{\Delta t}{2}\left(\mathbf{A}{}_{12}+\mathbf{B}{}_{12}\right)\mathbf{H}{}_{\mathrm{a}}^n.\hfill \end{array}$$

(35)

The last equation is due to the center bracketed terms in (34) that may be considered as (32) at one time step backward, i.e., $\mathbf{e}{}^{n-1}+\frac{\Delta t}{2}\mathbf{A}{}_{12}\mathbf{h}{}^{n-1}=\left(\frac{1}{2}\mathbf{I}{}_3-\frac{\Delta t^2}{8}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\right)\mathbf{E}{}_{\mathrm{a}}^{n-\frac{1}{2}}$. Multiplying two across both sides of (35), we obtain the leapfrog update equation for physical electric field:

$$\begin{array}{c}\left(\mathbf{I}{}_3-\frac{\Delta t^2}{4}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\right)\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}=\left(\mathbf{I}{}_3-\frac{\Delta t^2}{4}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\right)\mathbf{E}{}_{\mathrm{a}}^{n-\frac{1}{2}}+\Delta t\left(\mathbf{A}{}_{12}+\mathbf{B}{}_{12}\right)\mathbf{H}{}_{\mathrm{a}}^n.\end{array}$$

(36)

To find the corresponding equation for magnetic field, we substitute (22) into (24), as well as (27) at one time step backward into (25) to yield the auxiliary electric and magnetic fields at $n+\frac{1}{2}$ as

$$\begin{array}{c}\mathbf{e}{}^{n+\frac{1}{2}}=\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}-\mathbf{e}{}^{n-\frac{1}{2}}-\frac{\Delta t}{2}\mathbf{B}{}_{12}\mathbf{H}{}_{\mathrm{a}}^n\end{array}$$

(37)

$$\begin{array}{c}\mathbf{h}{}^{n+\frac{1}{2}}=\mathbf{H}{}_{\mathrm{a}}^n-\mathbf{h}{}^{n-\frac{1}{2}}+\frac{\Delta t}{2}\mathbf{A}{}_{21}\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}.\end{array}$$

(38)

Upon using (37) and (38), (26) can be manipulated as follows:

$$\begin{array}{cc}\hfill (\frac{1}{2}\mathbf{I}{}_3& -\frac{\Delta t^2}{8}\mathbf{B}{}_{21}\mathbf{B}{}_{12})\mathbf{H}{}_{\mathrm{a}}^{n+1}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & =\mathbf{h}{}^{n+\frac{1}{2}}+\frac{\Delta t}{2}\mathbf{B}{}_{21}\mathbf{e}{}^{n+\frac{1}{2}}\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill & =\mathbf{H}{}_{\mathrm{a}}^n-\mathbf{h}{}^{n-\frac{1}{2}}+\frac{\Delta t}{2}\mathbf{A}{}_{21}\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}+\frac{\Delta t}{2}\mathbf{B}{}_{21}\left(\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}-\mathbf{e}{}^{n-\frac{1}{2}}-\frac{\Delta t}{2}\mathbf{B}{}_{12}\mathbf{H}{}_{\mathrm{a}}^n\right)\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill & =\mathbf{H}{}_{\mathrm{a}}^n-\frac{\Delta t^2}{4}\mathbf{B}{}_{21}\mathbf{B}{}_{12}\mathbf{H}{}_{\mathrm{a}}^n-\left(\mathbf{h}{}^{n-\frac{1}{2}}+\frac{\Delta t}{2}\mathbf{B}{}_{21}\mathbf{e}{}^{n-\frac{1}{2}}\right)+\frac{\Delta t}{2}\left(\mathbf{A}{}_{21}+\mathbf{B}{}_{21}\right)\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill & =\left(\frac{1}{2}\mathbf{I}{}_3-\frac{\Delta t^2}{8}\mathbf{B}{}_{21}\mathbf{B}{}_{12}\right)\mathbf{H}{}_{\mathrm{a}}^n+\frac{\Delta t}{2}\left(\mathbf{A}{}_{21}+\mathbf{B}{}_{21}\right)\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}.\hfill \end{array}$$

(42)

Multiplying two across both sides of (42), we arrive at the leapfrog update equation for physical magnetic field:

$$\begin{array}{c}\left(\mathbf{I}{}_3-\frac{\Delta t^2}{4}\mathbf{B}{}_{21}\mathbf{B}{}_{12}\right)\mathbf{H}{}_{\mathrm{a}}^{n+1}=\left(\mathbf{I}{}_3-\frac{\Delta t^2}{4}\mathbf{B}{}_{21}\mathbf{B}{}_{12}\right)\mathbf{H}{}_{\mathrm{a}}^n+\Delta t\left(\mathbf{A}{}_{21}+\mathbf{B}{}_{21}\right)\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}.\end{array}$$

(43)

2.3. Leapfrog Fundamental Scheme for ADI FDTD Method

Equations (36) and (43) are the update equations of leapfrog ADI FDTD method in the conventional form [32,33]. Compared to the conventional ADI FDTD method with time-collocated fields in (6)–(9), the leapfrog ADI FDTD method has only time-staggered fields $\mathbf{E}{}_{\mathrm{a}}$ and $\mathbf{H}{}_{\mathrm{a}}$ at half- and full-integer time-indices, respectively. However, the RHS of (36) and (43) still involve considerable operations especially due to the complicated products of matrix operator partitions $\mathbf{A}{}_{12}\mathbf{A}{}_{21}$ and $\mathbf{B}{}_{21}\mathbf{B}{}_{12}$. To simplify the RHS, we resort to the principle of fundamental schemes as before [30,31]. In particular, introducing the auxiliary variables $\mathbf{e}{}_{\mathrm{a}}$ and $\mathbf{h}{}_{\mathrm{a}}$, the update procedures can be written as

$$\left(\mathbf{I}{}_3-\frac{\Delta t^2}{4}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\right)\mathbf{e}{}_{\mathrm{a}}^n=\Delta t\left(\mathbf{A}{}_{12}+\mathbf{B}{}_{12}\right)\mathbf{H}{}_{\mathrm{a}}^n$$

(44)

$$\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}=\mathbf{E}{}_{\mathrm{a}}^{n-\frac{1}{2}}+\mathbf{e}{}_{\mathrm{a}}^n$$

(45)

$$\left(\mathbf{I}{}_3-\frac{\Delta t^2}{4}\mathbf{B}{}_{21}\mathbf{B}{}_{12}\right)\mathbf{h}{}_{\mathrm{a}}^{n+\frac{1}{2}}=\Delta t\left(\mathbf{A}{}_{21}+\mathbf{B}{}_{21}\right)\mathbf{E}{}_{\mathrm{a}}^{n+\frac{1}{2}}$$

(46)

$$\mathbf{H}{}_{\mathrm{a}}^{n+1}=\mathbf{H}{}_{\mathrm{a}}^n+\mathbf{h}{}_{\mathrm{a}}^{n+\frac{1}{2}}.$$

(47)

Compared to the conventional form in (36) and (43), the RHS of (44)–(47) no longer contain the complicated product of matrix operator partitions, i.e., no more $\mathbf{A}{}_{12}\mathbf{A}{}_{21}$ & $\mathbf{B}{}_{21}\mathbf{B}{}_{12}$. These lead to simplifications of the implicit update procedures with reduced flops and improved efficiency. Like all the fundamental schemes, their RHS contain only each of the intrinsic matrix operator partitions, i.e., $\mathbf{A}{}_{12}$, $\mathbf{A}{}_{21}$, $\mathbf{B}{}_{12}$, $\mathbf{B}{}_{21}$, thus one may regard (44)–(47) as the leapfrog fundamental scheme for ADI FDTD method, or in short, leapfrog FADI FDTD method. Note that for this leapfrog scheme with time-staggered electric and magnetic fields, the remaining physical fields cannot be recovered easily without additional auxiliary ones like for the time-collocated schemes using (21) and (28).

Our earlier investigations of the leapfrog ADI FDTD method have found several issues that are seldom discussed by many. Although being leapfrog like explicit scheme, the method does not have complying divergence and violates Gauss’s law for both 2-D and 3-D [47]. Moreover, the method formulated for lossy media suffers from the field leakage problem that would yield non-zero field values even in/thru perfect electric conducting (PEC) media, i.e., when the conductivity approaches infinity. This may be explained by the infeasibility of magnetic field to incorporate the electric conductivity effect in the second implicit update procedure, cf. (43) or (46) and (47). To alleviate the issue, one simple approach is by using high permittivity values to achieve similar total reflection characteristics for incident waves into non-penetrable targets like PEC [52]. Alternative implicit schemes for lossy media may be adapted to resolve the field leakage problem [53]. In addition, there have also been further developments of implicit schemes involving time-collocated fields, such as higher order, dispersionless and parameter optimized methods [54,55,56,57,58,59,60]. Efforts have been under way to carry out the corresponding developments for leapfrog scheme involving time-staggered fields, especially using the leapfrog fundamental form with simpler, more concise and efficient RHS.

3. From Time-Collocated to Leapfrog Fundamental Schemes for CDI FDTD Method

3.1. Time-Collocated Fundamental Scheme for CDI FDTD Method

Besides the ADI FDTD method above, we shall also consider alternative unconditionally stable methods. The conventional LOD FDTD method comprises two update procedures as [7,8,9]

$$\left(\mathbf{I}{}_6-\frac{\Delta t}{2}\mathbf{A}{}_{}\right)\mathbf{u}{}_{\mathrm{l}}^{n+\frac{1}{2}}=\left(\mathbf{I}{}_6+\frac{\Delta t}{2}\mathbf{A}{}_{}\right)\mathbf{u}{}_{\mathrm{l}}^n$$

(48)

$$\left(\mathbf{I}{}_6-\frac{\Delta t}{2}\mathbf{B}{}_{}\right)\mathbf{u}{}_{\mathrm{l}}^{n+1}=\left(\mathbf{I}{}_6+\frac{\Delta t}{2}\mathbf{B}{}_{}\right)\mathbf{u}{}_{\mathrm{l}}^{n+\frac{1}{2}}$$

(49)

where the field vectors are subscripted with ‘l’, $\mathbf{A}$ and $\mathbf{B}$ are the same matrix operators as before, cf. (3). The LOD FDTD method is only first-order accurate in time and is also called split-step FDTD of first order (SS1). It can be extended readily for higher order, compact, multi-stage, GPU/parallel/cloud computations, etc. [7,8,9,10,11,61,62,63,64,65,66,67,68].

While the LHS of (48) and (49) are the same as those of (1) and (2), their RHS differ with (slightly) reduced operations by adopting the same types of LHS operators. This is evident from the update equations that can be expressed in terms of matrix operator partitions as

$$\left(\mathbf{I}{}_3-\frac{\Delta t^2}{4}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\right)\mathbf{E}{}_{\mathrm{l}}^{n+\frac{1}{2}}=\mathbf{E}{}_{\mathrm{l}}^n+\frac{\Delta t^2}{4}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\mathbf{E}{}_{\mathrm{l}}^n+\Delta t\mathbf{A}{}_{12}\mathbf{H}{}_{\mathrm{l}}^n$$

(50)

$$\mathbf{H}{}_{\mathrm{l}}^{n+\frac{1}{2}}=\mathbf{H}{}_{\mathrm{l}}^n+\frac{\Delta t}{2}\mathbf{A}{}_{21}\left(\mathbf{E}{}_{\mathrm{l}}^{n+\frac{1}{2}}+\mathbf{E}{}_{\mathrm{l}}^n\right)$$

(51)

$$\left(\mathbf{I}{}_3-\frac{\Delta t^2}{4}\mathbf{B}{}_{12}\mathbf{B}{}_{21}\right)\mathbf{E}{}_{\mathrm{l}}^{n+1}=\mathbf{E}{}_{\mathrm{l}}^{n+\frac{1}{2}}+\frac{\Delta t^2}{4}\mathbf{B}{}_{12}\mathbf{B}{}_{21}\mathbf{E}{}_{\mathrm{l}}^{n+\frac{1}{2}}+\Delta t\mathbf{B}{}_{12}\mathbf{H}{}_{\mathrm{l}}^{n+\frac{1}{2}}$$

(52)

$$\mathbf{H}{}_{\mathrm{l}}^{n+1}=\mathbf{H}{}_{\mathrm{l}}^{n+\frac{1}{2}}+\frac{\Delta t}{2}\mathbf{B}{}_{21}\left(\mathbf{E}{}_{\mathrm{l}}^{n+\frac{1}{2}}+\mathbf{E}{}_{\mathrm{l}}^{n+1}\right).$$

(53)

Equations (50)–(53) with fewer RHS terms incur (slightly) less flops count compared to (6)–(9), but there are still complicated products of matrix operator partitions $\mathbf{A}{}_{12}\mathbf{A}{}_{21}$ and $\mathbf{B}{}_{12}\mathbf{B}{}_{21}$.

To further reduce the flops, we apply the principle of fundamental schemes for (48) and (49) to arrive at [30,31]

$$\left(\frac{1}{2}\mathbf{I}{}_6-\frac{\Delta t}{4}\mathbf{A}{}_{}\right)\mathbf{v}{}_{\mathrm{l}}^{n+\frac{1}{2}}=\mathbf{u}{}_{\mathrm{l}}^n$$

(54)

$$\mathbf{u}{}_{\mathrm{l}}^{n+\frac{1}{2}}=\mathbf{v}{}_{\mathrm{l}}^{n+\frac{1}{2}}-\mathbf{u}{}_{\mathrm{l}}^n$$

(55)

$$\left(\frac{1}{2}\mathbf{I}{}_6-\frac{\Delta t}{4}\mathbf{B}{}_{}\right)\mathbf{v}{}_{\mathrm{l}}^{n+1}=\mathbf{u}{}_{\mathrm{l}}^{n+\frac{1}{2}}$$

(56)

$$\mathbf{u}{}_{\mathrm{l}}^{n+1}=\mathbf{v}{}_{\mathrm{l}}^{n+1}-\mathbf{u}{}_{\mathrm{l}}^{n+\frac{1}{2}}.$$

(57)

Again through the auxiliary field vectors $\mathbf{v}{}_{\mathrm{l}}$’s at full- and half-integer time indices: $\mathbf{v}{}_{\mathrm{l}}={\left[\mathbf{e}{}_{\mathrm{l}}\mathbf{h}{}_{\mathrm{l}}\right]}^T$, the RHS of (54)–(57) become free of any matrix operator and contain only vectors. Equations (54)–(57) constitute our fundamental LOD FDTD method, or in short, FLOD FDTD method, which is comparable to the FADI FDTD method in (10)–(13). Based on the matrix operator partitions (4) and (5), the implementation of FLOD FDTD method can be carried out using the update equations as follows:

$$\left(\frac{1}{2}\mathbf{I}{}_3-\frac{\Delta t^2}{8}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\right)\mathbf{e}{}_{\mathrm{l}}^{n+\frac{1}{2}}=\mathbf{E}{}_{\mathrm{l}}^n+\frac{\Delta t}{2}\mathbf{A}{}_{12}\mathbf{H}{}_{\mathrm{l}}^n$$

(58)

$$\mathbf{H}{}_{\mathrm{l}}^{n+\frac{1}{2}}=\mathbf{H}{}_{\mathrm{l}}^n+\frac{\Delta t}{2}\mathbf{A}{}_{21}\mathbf{e}{}_{\mathrm{l}}^{n+\frac{1}{2}}$$

(59)

$$\mathbf{E}{}_{\mathrm{l}}^{n+\frac{1}{2}}=\mathbf{e}{}_{\mathrm{l}}^{n+\frac{1}{2}}-\mathbf{E}{}_{\mathrm{l}}^n$$

(60)

$$\left(\frac{1}{2}\mathbf{I}{}_3-\frac{\Delta t^2}{8}\mathbf{B}{}_{21}\mathbf{B}{}_{12}\right)\mathbf{h}{}_{\mathrm{l}}^{n+1}=\mathbf{H}{}_{\mathrm{l}}^{n+\frac{1}{2}}+\frac{\Delta t}{2}\mathbf{B}{}_{21}\mathbf{E}{}_{\mathrm{l}}^{n+\frac{1}{2}}$$

(61)

$$\mathbf{E}{}_{\mathrm{l}}^{n+1}=\mathbf{E}{}_{\mathrm{l}}^{n+\frac{1}{2}}+\frac{\Delta t}{2}\mathbf{B}{}_{12}\mathbf{h}{}_{\mathrm{l}}^{n+1}$$

(62)

$$\mathbf{H}{}_{\mathrm{l}}^{n+1}=\mathbf{h}{}_{\mathrm{l}}^{n+1}-\mathbf{H}{}_{\mathrm{l}}^{n+\frac{1}{2}}.$$

(63)

Equations (58)–(63) for FLOD FDTD method are analogous to the fundamental implicit E-H scheme for FADI FDTD method. One can also derive the corresponding equations analogous to the above fundamental implicit E-E or other schemes.

To improve the temporal accuracy to second order, one needs to perform proper input and output processings for $n=0$ and when output is required, respectively, cf. [9]. Alternatively in addition to the second-order accuracy, the following processings will also lead to more desirable characteristics of complying divergence [11]:

–

Input processing:

$$\begin{array}{c}\left(\mathbf{I}{}_6-\frac{\Delta t}{2}\mathbf{B}{}_{}\right)\mathbf{u}{}_{\mathrm{l}}^0=\mathbf{u}{}_{\mathrm{c}}^0\end{array}$$

(64)

–

Output processing:

$$\begin{array}{c}\mathbf{u}{}_{\mathrm{c}}^{n+1}=\left(\mathbf{I}{}_6-\frac{\Delta t}{2}\mathbf{B}{}_{}\right)\mathbf{u}{}_{\mathrm{l}}^{n+1}.\end{array}$$

(65)

In particular, (64) and (65) result in the output fields $\mathbf{u}{}_{\mathrm{c}}$’s that feature complying divergence and satisfy Gauss’s law. These equations together with (54)–(57) have been previously called the fundamental LOD FDTD method with complying divergence, or in short, FLOD-CD FDTD method. In view of the close connection to the recent leapfrog CDI FDTD method below, the method may be more aptly regarded as the time-collocated fundamental scheme for CDI FDTD method. The subscripts ‘c’ signify the quantities for CDI FDTD with ‘CD’ denoting ‘complying divergence’.

3.2. From Time-Collocated to Leapfrog Schemes for CDI FDTD Method

Due to their connection that is not apparent, we shall provide the formulations from time-collocated to leapfrog schemes for the CDI FDTD method. From (49), (56) and (65), we have

$$\mathbf{u}{}_{\mathrm{c}}^{n+1}=\left(\mathbf{I}{}_6+\frac{\Delta t}{2}\mathbf{B}{}_{}\right)\mathbf{u}{}_{\mathrm{l}}^{n+\frac{1}{2}}$$

(66)

$$\begin{array}{cc}\hfill & =\left(\mathbf{I}{}_6+\frac{\Delta t}{2}\mathbf{B}{}_{}\right)\left(\frac{1}{2}\mathbf{I}{}_6-\frac{\Delta t}{4}\mathbf{B}{}_{}\right)\mathbf{v}{}_{\mathrm{l}}^{n+1}.\hfill \end{array}$$

(67)

This equation gives the expressions of auxiliary fields in terms of the physical ones at $n+1$ as

$$\begin{array}{c}\mathbf{e}{}_{\mathrm{l}}^{n+1}=2{\left(\mathbf{I}{}_3-\frac{\Delta t}{4}\mathbf{B}{}_{12}\mathbf{B}{}_{21}\right)}^{-1}\mathbf{E}{}_{\mathrm{c}}^{n+1}\end{array}$$

(68)

$$\begin{array}{c}\mathbf{h}{}_{\mathrm{l}}^{n+1}=2{\left(\mathbf{I}{}_3-\frac{\Delta t}{4}\mathbf{B}{}_{21}\mathbf{B}{}_{12}\right)}^{-1}\mathbf{H}{}_{\mathrm{c}}^{n+1}.\end{array}$$

(69)

Besides the physical fields at integer time indices (65), one can also define the corresponding ones at half-integer time indices along with (56) as

$$\mathbf{u}{}_{\mathrm{c}}^{n+\frac{1}{2}}=\left(\mathbf{I}{}_6-\frac{\Delta t}{2}\mathbf{A}{}_{}\right)\mathbf{u}{}_{\mathrm{l}}^{n+\frac{1}{2}}$$

(70)

$$\begin{array}{cc}\hfill & =\left(\mathbf{I}{}_6-\frac{\Delta t}{2}\mathbf{A}{}_{}\right)\left(\frac{1}{2}\mathbf{I}{}_6-\frac{\Delta t}{4}\mathbf{B}{}_{}\right)\mathbf{v}{}_{\mathrm{l}}^{n+1}.\hfill \end{array}$$

(71)

Upon substituting (68) and (69) into (71), the electric field at $n+\frac{1}{2}$ is found to be

$$\begin{array}{cc}\hfill \mathbf{E}{}_{\mathrm{c}}^{n+\frac{1}{2}}& =\frac{1}{2}\left(\mathbf{I}+\frac{\Delta t}{4}\mathbf{A}{}_{12}\mathbf{B}{}_{21}\right)\mathbf{e}{}_{\mathrm{l}}^{n+1}-\frac{\Delta t}{4}\left(\mathbf{A}{}_{12}+\mathbf{B}{}_{12}\right)\mathbf{h}{}_{\mathrm{l}}^{n+1}\hfill \\ \hfill & =\left(\mathbf{I}+\frac{\Delta t}{4}\mathbf{A}{}_{12}\mathbf{B}{}_{21}\right){\left(\mathbf{I}-\frac{\Delta t}{4}\mathbf{B}{}_{12}\mathbf{B}{}_{21}\right)}^{-1}\mathbf{E}{}_{\mathrm{c}}^{n+1}\hfill \end{array}$$

(72)

$$\begin{array}{cc}\hfill & -\frac{\Delta t}{2}\left(\mathbf{A}{}_{12}+\mathbf{B}{}_{12}\right){\left(\mathbf{I}-\frac{\Delta t}{4}\mathbf{B}{}_{21}\mathbf{B}{}_{12}\right)}^{-1}\mathbf{H}{}_{\mathrm{c}}^{n+1}.\hfill \end{array}$$

(73)

Taking (73) at one time step backward, the electric field at $n-\frac{1}{2}$ simply reads

$$\begin{array}{cc}\hfill \mathbf{E}{}_{\mathrm{c}}^{n-\frac{1}{2}}=\left(\mathbf{I}+\frac{\Delta t}{4}\mathbf{A}{}_{12}\mathbf{B}{}_{21}\right){\left(\mathbf{I}-\frac{\Delta t}{4}\mathbf{B}{}_{12}\mathbf{B}{}_{21}\right)}^{-1}& {\mathbf{E}}_{\mathrm{c}}^n\hfill \\ \hfill -& \frac{\Delta t}{2}\left(\mathbf{A}{}_{12}+\mathbf{B}{}_{12}\right){\left(\mathbf{I}-\frac{\Delta t}{4}\mathbf{B}{}_{21}\mathbf{B}{}_{12}\right)}^{-1}\mathbf{H}{}_{\mathrm{c}}^n.\hfill \end{array}$$

(74)

Also taking (65) at one time step backward, we have its inverse relation

$$\begin{array}{c}\mathbf{u}{}_{\mathrm{l}}^n={\left(\mathbf{I}-\frac{\Delta t}{2}\mathbf{B}{}_{}\right)}^{-1}\mathbf{u}{}_{\mathrm{c}}^n.\end{array}$$

(75)

Substituting this relation into (70) or (48), the physical fields at half-integer time indices can be expressed by

$$\mathbf{u}{}_{\mathrm{c}}^{n+\frac{1}{2}}=\left(\mathbf{I}+\frac{\Delta t}{2}\mathbf{A}{}_{}\right)\mathbf{u}{}_{\mathrm{l}}^n$$

(76)

$$\begin{array}{cc}\hfill & =\left(\mathbf{I}+\frac{\Delta t}{2}\mathbf{A}{}_{}\right){\left(\mathbf{I}-\frac{\Delta t}{2}\mathbf{B}{}_{}\right)}^{-1}\mathbf{u}{}_{\mathrm{c}}^n.\hfill \end{array}$$

(77)

Then the electric field at $n+\frac{1}{2}$ can be written in the form alternative to (73) as

$$\begin{array}{cc}\hfill \mathbf{E}{}_{\mathrm{c}}^{n+\frac{1}{2}}=\left(\mathbf{I}+\frac{\Delta t}{4}\mathbf{A}{}_{12}\mathbf{B}{}_{21}\right){\left(\mathbf{I}-\frac{\Delta t}{4}\mathbf{B}{}_{12}\mathbf{B}{}_{21}\right)}^{-1}& {\mathbf{E}}_{\mathrm{c}}^n\hfill \\ \hfill +& \frac{\Delta t}{2}\left(\mathbf{A}{}_{12}+\mathbf{B}{}_{12}\right){\left(\mathbf{I}-\frac{\Delta t}{4}\mathbf{B}{}_{21}\mathbf{B}{}_{12}\right)}^{-1}\mathbf{H}{}_{\mathrm{c}}^n.\hfill \end{array}$$

(78)

Subtracting (74) from (78), we arrive at the leapfrog update equation for electric field:

$$\begin{array}{c}\mathbf{E}{}_{\mathrm{c}}^{n+\frac{1}{2}}-\mathbf{E}{}_{\mathrm{c}}^{n-\frac{1}{2}}=\Delta t\left(\mathbf{A}{}_{12}+\mathbf{B}{}_{12}\right){\left(\mathbf{I}-\frac{\Delta t}{4}\mathbf{B}{}_{21}\mathbf{B}{}_{12}\right)}^{-1}\mathbf{H}{}_{\mathrm{c}}^n.\end{array}$$

(79)

To derive the corresponding equation for magnetic field, we use the inverse relation of (70):

$$\begin{array}{c}\mathbf{u}{}_{\mathrm{l}}^{n+\frac{1}{2}}={\left(\mathbf{I}-\frac{\Delta t}{2}\mathbf{A}{}_{}\right)}^{-1}\mathbf{u}{}_{\mathrm{c}}^{n+\frac{1}{2}}\end{array}$$

(80)

and substitute into (66) as

$$\begin{array}{c}\mathbf{u}{}_{\mathrm{c}}^{n+1}=\left(\mathbf{I}+\frac{\Delta t}{2}\mathbf{B}{}_{}\right){\left(\mathbf{I}-\frac{\Delta t}{2}\mathbf{A}{}_{}\right)}^{-1}\mathbf{u}{}_{\mathrm{c}}^{n+\frac{1}{2}}.\end{array}$$

(81)

From this equation, the magnetic field at $n+1$ is found to be

$$\begin{array}{cc}\hfill \mathbf{H}{}_{\mathrm{c}}^{n+1}=\left(\mathbf{I}+\frac{\Delta t}{4}\mathbf{B}{}_{21}\mathbf{A}{}_{12}\right){\left(\mathbf{I}-\frac{\Delta t}{4}\mathbf{A}{}_{21}\mathbf{A}{}_{12}\right)}^{-1}& {\mathbf{H}}_{\mathrm{c}}^{n+\frac{1}{2}}\hfill \\ \hfill +\frac{\Delta t}{2}& \left(\mathbf{A}{}_{21}+\mathbf{B}{}_{21}\right){\left(\mathbf{I}-\frac{\Delta t}{4}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\right)}^{-1}\mathbf{E}{}_{\mathrm{c}}^{n+\frac{1}{2}}.\hfill \end{array}$$

(82)

Using (54) in (76), we have

$$\begin{array}{c}\mathbf{u}{}_{\mathrm{c}}^{n+\frac{1}{2}}=\left(\mathbf{I}+\frac{\Delta t}{2}\mathbf{A}{}_{}\right)\left(\frac{1}{2}\mathbf{I}-\frac{\Delta t}{4}\mathbf{A}{}_{}\right)\mathbf{v}{}_{\mathrm{l}}^{n+\frac{1}{2}}.\end{array}$$

(83)

This equation gives the expressions of auxiliary fields in terms of the physical ones at $n+\frac{1}{2}$ as

$$\begin{array}{c}\mathbf{e}{}_{\mathrm{l}}^{n+\frac{1}{2}}=2{\left(\mathbf{I}-\frac{\Delta t}{4}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\right)}^{-1}\mathbf{E}{}_{\mathrm{c}}^{n+\frac{1}{2}}\end{array}$$

(84)

$$\begin{array}{c}\mathbf{h}{}_{\mathrm{l}}^{n+\frac{1}{2}}=2{\left(\mathbf{I}-\frac{\Delta t}{4}\mathbf{A}{}_{21}\mathbf{A}{}_{12}\right)}^{-1}\mathbf{H}{}_{\mathrm{c}}^{n+\frac{1}{2}}.\end{array}$$

(85)

Again, taking (65) at n and using (54) yields

$$\mathbf{u}{}_{\mathrm{c}}^n=\left(\mathbf{I}-\frac{\Delta t}{2}\mathbf{B}{}_{}\right)\mathbf{u}{}_{\mathrm{l}}^n$$

(86)

$$\begin{array}{cc}\hfill & =\left(\mathbf{I}-\frac{\Delta t}{2}\mathbf{B}{}_{}\right)\left(\frac{1}{2}\mathbf{I}-\frac{\Delta t}{4}\mathbf{A}{}_{}\right)\mathbf{v}{}_{\mathrm{l}}^{n+\frac{1}{2}}.\hfill \end{array}$$

(87)

Substituting (84) and (85) into (87), the magnetic field at n is found to be

$$\begin{array}{cc}\hfill \mathbf{H}{}_{\mathrm{c}}^n& =\frac{1}{2}\left(\mathbf{I}+\frac{\Delta t}{4}\mathbf{B}{}_{21}\mathbf{A}{}_{12}\right)\mathbf{h}{}_{\mathrm{l}}^{n+\frac{1}{2}}-\frac{\Delta t}{4}\left(\mathbf{A}{}_{21}+\mathbf{B}{}_{21}\right)\mathbf{e}{}_{\mathrm{l}}^{n+\frac{1}{2}}\hfill \\ \hfill & =\left(\mathbf{I}+\frac{\Delta t}{4}\mathbf{B}{}_{21}\mathbf{A}{}_{12}\right){\left(\mathbf{I}-\frac{\Delta t}{4}\mathbf{A}{}_{21}\mathbf{A}{}_{12}\right)}^{-1}\mathbf{H}{}_{\mathrm{c}}^{n+\frac{1}{2}}\hfill \end{array}$$

(88)

$$\begin{array}{cc}\hfill & -\frac{\Delta t}{2}\left(\mathbf{A}{}_{21}+\mathbf{B}{}_{21}\right){\left(\mathbf{I}-\frac{\Delta t}{4}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\right)}^{-1}\mathbf{E}{}_{\mathrm{c}}^{n+\frac{1}{2}}.\hfill \end{array}$$

(89)

Subtracting (89) from (82), we obtain the leapfrog update equation for magnetic field:

$$\begin{array}{c}\mathbf{H}{}_{\mathrm{c}}^{n+1}-\mathbf{H}{}_{\mathrm{c}}^n=\Delta t\left(\mathbf{A}{}_{21}+\mathbf{B}{}_{21}\right){\left(\mathbf{I}-\frac{\Delta t}{4}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\right)}^{-1}\mathbf{E}{}_{\mathrm{c}}^{n+\frac{1}{2}}.\end{array}$$

(90)

3.3. Leapfrog Fundamental Scheme for CDI FDTD Method

Equations (79) and (90) can be implemented in the following update procedures:

$$\left(\mathbf{I}-\frac{\Delta t^2}{4}\mathbf{B}{}_{21}\mathbf{B}{}_{12}\right)\mathbf{h}{}_{\mathrm{c}}^n=\mathbf{H}{}_{\mathrm{c}}^n$$

(91)

$$\mathbf{E}{}_{\mathrm{c}}^{n+\frac{1}{2}}=\mathbf{E}{}_{\mathrm{c}}^{n-\frac{1}{2}}+\Delta t\left(\mathbf{A}{}_{12}+\mathbf{B}{}_{12}\right)\mathbf{h}{}_{\mathrm{c}}^n$$

(92)

$$\left(\mathbf{I}-\frac{\Delta t^2}{4}\mathbf{A}{}_{12}\mathbf{A}{}_{21}\right)\mathbf{e}{}_{\mathrm{c}}^{n+\frac{1}{2}}=\mathbf{E}{}_{\mathrm{c}}^{n+\frac{1}{2}}$$

(93)

$$\mathbf{H}{}_{\mathrm{c}}^{n+1}=\mathbf{H}{}_{\mathrm{c}}^n+\Delta t\left(\mathbf{A}{}_{21}+\mathbf{B}{}_{21}\right)\mathbf{e}{}_{\mathrm{c}}^{n+\frac{1}{2}}.$$

(94)

Here, $\mathbf{E}{}_{\mathrm{c}}$ and $\mathbf{H}{}_{\mathrm{c}}$ are the physical variables for leapfrog CDI FDTD fields, while $\mathbf{e}{}_{\mathrm{c}}$ and $\mathbf{h}{}_{\mathrm{c}}$ are the intermediate auxiliary ones subscripted with ‘c’ (which are different from those subscripted with ‘l’, ‘a’, or unsubscripted). Equations (91)–(94) constitute the leapfrog fundamental scheme for CDI FDTD method, which features complying divergence and satisfies Gauss’s law [49]. Notice that the implicit update procedures in (91) and (93) have the simplest fundamental form with RHS completely free of matrix operators. The explicit update procedures in (92) and (94) are compatible with those of leapfrog explicit FDTD method, except that they now involve operations with auxiliary fields instead of physical ones.

Based on their equations in fundamental forms, both leapfrog fundamental schemes for ADI and CDI FDTD methods can be seen to be related closely. In particular, their update procedures appear to be convertible readily by shifting the RHS operators of (44) to (92) and (46) to (94), while also changing their operands accordingly between physical and auxiliary fields. Such insights into the relations would not be evident if one refers only to the conventional form of leapfrog ADI FDTD method in (36) and (43) [without the fundamental form in (44)–(47)], or the time-collocated fundamental scheme for CDI FDTD method in (54)–(65). Meanwhile, pertaining to the field leakage problem of leapfrog ADI FDTD method [52], the previous simple approach using high permittivity values for non-penetrable targets is also applicable for the leapfrog CDI FDTD method. Moreover, with the presence of both electric and magnetic fields in the implicit and explicit update equations, it is feasible to incorporate the electric conductivity effect in both update procedures of the leapfrog CDI FDTD method. Besides, other variants of implicit schemes for lossy media [53] may be adapted for the method as well. With its many desirable characteristics, the leapfrog fundamental CDI FDTD method is very promising for further developments following many previous works, e.g., multiconductor transmission lines [69,70,71,72].

4. Discussion

To facilitate subsequent discussions, let us write out the component equations in detail for various unconditionally stable implicit FDTD schemes. Based on (36) and (43), we have the detailed update procedures of leapfrog ADI FDTD method in the conventional form:

-

First procedure for $\mathbf{E}{}_{\mathrm{a}}$ from $n-\frac{1}{2}$ to $n+\frac{1}{2}$:

$$\begin{array}{c}\hfill E_{\mathrm{a}x}^{n+\frac{1}{2}}-\frac{\Delta t^2}{4ϵ}\frac{\partial }{\partial y}\frac{1}{\mu }\frac{\partial }{\partial y}{E}_{\mathrm{a}x}^{n+\frac{1}{2}}=E_{\mathrm{a}x}^{n-\frac{1}{2}}-\frac{\Delta t^2}{4ϵ}\frac{\partial }{\partial y}\frac{1}{\mu }\frac{\partial }{\partial y}{E}_{\mathrm{a}x}^{n-\frac{1}{2}}+\frac{\Delta t}{ϵ}\frac{\partial }{\partial y}{H}_{\mathrm{a}z}^n-\frac{\Delta t}{ϵ}\frac{\partial }{\partial z}{H}_{\mathrm{a}y}^n\end{array}$$

(95)

$$\begin{array}{c}\hfill E_{\mathrm{a}y}^{n+\frac{1}{2}}-\frac{\Delta t^2}{4ϵ}\frac{\partial }{\partial z}\frac{1}{\mu }\frac{\partial }{\partial z}{E}_{\mathrm{a}y}^{n+\frac{1}{2}}=E_{\mathrm{a}y}^{n-\frac{1}{2}}-\frac{\Delta t^2}{4ϵ}\frac{\partial }{\partial z}\frac{1}{\mu }\frac{\partial }{\partial z}{E}_{\mathrm{a}y}^{n-\frac{1}{2}}+\frac{\Delta t}{ϵ}\frac{\partial }{\partial z}{H}_{\mathrm{a}x}^n-\frac{\Delta t}{ϵ}\frac{\partial }{\partial x}{H}_{\mathrm{a}z}^n\end{array}$$

(96)

$$\begin{array}{c}\hfill E_{\mathrm{a}z}^{n+\frac{1}{2}}-\frac{\Delta t^2}{4ϵ}\frac{\partial }{\partial x}\frac{1}{\mu }\frac{\partial }{\partial x}{E}_{\mathrm{a}z}^{n+\frac{1}{2}}=E_{\mathrm{a}z}^{n-\frac{1}{2}}-\frac{\Delta t^2}{4ϵ}\frac{\partial }{\partial x}\frac{1}{\mu }\frac{\partial }{\partial x}{E}_{\mathrm{a}z}^{n-\frac{1}{2}}+\frac{\Delta t}{ϵ}\frac{\partial }{\partial x}{H}_{\mathrm{a}y}^n-\frac{\Delta t}{ϵ}\frac{\partial }{\partial y}{H}_{\mathrm{a}x}^n\end{array}$$

(97)

-

Second procedure for $\mathbf{H}{}_{\mathrm{a}}$ from n to $n+1$:

$$\begin{array}{c}\hfill H_{\mathrm{a}x}^{n+1}-\frac{\Delta t^2}{4\mu }\frac{\partial }{\partial y}\frac{1}{ϵ}\frac{\partial }{\partial y}{H}_{\mathrm{a}x}^{n+1}=H_{\mathrm{a}x}^n-\frac{\Delta t^2}{4\mu }\frac{\partial }{\partial y}\frac{1}{ϵ}\frac{\partial }{\partial y}{H}_{\mathrm{a}x}^n+\frac{\Delta t}{\mu }\frac{\partial }{\partial z}{E}_{\mathrm{a}y}^{n+\frac{1}{2}}-\frac{\Delta t}{\mu }\frac{\partial }{\partial y}{E}_{\mathrm{a}z}^{n+\frac{1}{2}}\end{array}$$

(98)

$$\begin{array}{c}\hfill H_{\mathrm{a}y}^{n+1}-\frac{\Delta t^2}{4\mu }\frac{\partial }{\partial z}\frac{1}{ϵ}\frac{\partial }{\partial z}{H}_{\mathrm{a}y}^{n+1}=H_{\mathrm{a}y}^n-\frac{\Delta t^2}{4\mu }\frac{\partial }{\partial z}\frac{1}{ϵ}\frac{\partial }{\partial z}{H}_{\mathrm{a}y}^n+\frac{\Delta t}{\mu }\frac{\partial }{\partial x}{E}_{\mathrm{a}z}^{n+\frac{1}{2}}-\frac{\Delta t}{\mu }\frac{\partial }{\partial z}{E}_{\mathrm{a}x}^{n+\frac{1}{2}}\end{array}$$

(99)

$$\begin{array}{c}\hfill H_{\mathrm{a}z}^{n+1}-\frac{\Delta t^2}{4\mu }\frac{\partial }{\partial x}\frac{1}{ϵ}\frac{\partial }{\partial x}{H}_{\mathrm{a}z}^{n+1}=H_{\mathrm{a}z}^n-\frac{\Delta t^2}{4\mu }\frac{\partial }{\partial x}\frac{1}{ϵ}\frac{\partial }{\partial x}{H}_{\mathrm{a}z}^n+\frac{\Delta t}{\mu }\frac{\partial }{\partial y}{E}_{\mathrm{a}x}^{n+\frac{1}{2}}-\frac{\Delta t}{\mu }\frac{\partial }{\partial x}{E}_{\mathrm{a}y}^{n+\frac{1}{2}}.\end{array}$$

(100)

Upon resorting to the principle of fundamental schemes, we have the detailed update procedures of leapfrog FADI FDTD method from (44)–(47):

-

First procedure for $\mathbf{E}{}_{\mathrm{a}}$ from $n-\frac{1}{2}$ to $n+\frac{1}{2}$:

$$\begin{array}{c}{e}_{\mathrm{a}x}^{n+\frac{1}{2}}-\frac{\Delta t^2}{4ϵ}\frac{\partial }{\partial y}\frac{1}{\mu }\frac{\partial }{\partial y}{e}_{\mathrm{a}x}^{n+\frac{1}{2}}=\frac{\Delta t}{ϵ}\frac{\partial }{\partial y}{H}_{\mathrm{a}z}^n-\frac{\Delta t}{ϵ}\frac{\partial }{\partial z}{H}_{\mathrm{a}y}^n\end{array}$$

(101)

$$\begin{array}{c}{e}_{\mathrm{a}y}^{n+\frac{1}{2}}-\frac{\Delta t^2}{4ϵ}\frac{\partial }{\partial z}\frac{1}{\mu }\frac{\partial }{\partial z}{e}_{\mathrm{a}y}^{n+\frac{1}{2}}=\frac{\Delta t}{ϵ}\frac{\partial }{\partial z}{H}_{\mathrm{a}x}^n-\frac{\Delta t}{ϵ}\frac{\partial }{\partial x}{H}_{\mathrm{a}z}^n\end{array}$$

(102)

$$\begin{array}{c}{e}_{\mathrm{a}z}^{n+\frac{1}{2}}-\frac{\Delta t^2}{4ϵ}\frac{\partial }{\partial x}\frac{1}{\mu }\frac{\partial }{\partial x}{e}_{\mathrm{a}z}^{n+\frac{1}{2}}=\frac{\Delta t}{ϵ}\frac{\partial }{\partial x}{H}_{\mathrm{a}y}^n-\frac{\Delta t}{ϵ}\frac{\partial }{\partial y}{H}_{\mathrm{a}x}^n\end{array}$$

(103)

$$\begin{array}{c}{E}_{\mathrm{a}x}^{n+\frac{1}{2}}=E_{\mathrm{a}x}^{n-\frac{1}{2}}+e_{\mathrm{a}x}^{n+\frac{1}{2}}\end{array}$$

(104)

$$\begin{array}{c}{E}_{\mathrm{a}y}^{n+\frac{1}{2}}=E_{\mathrm{a}y}^{n-\frac{1}{2}}+e_{\mathrm{a}y}^{n+\frac{1}{2}}\end{array}$$

(105)

$$\begin{array}{c}{E}_{\mathrm{a}z}^{n+\frac{1}{2}}=E_{\mathrm{a}z}^{n-\frac{1}{2}}+e_{\mathrm{a}z}^{n+\frac{1}{2}}\end{array}$$

(106)

Second procedure for $\mathbf{H}{}_{\mathrm{a}}$ from n to $n+1$:

$$\begin{array}{c}{h}_{\mathrm{a}x}^{n+1}-\frac{\Delta t^2}{4\mu }\frac{\partial }{\partial y}\frac{1}{ϵ}\frac{\partial }{\partial y}{h}_{\mathrm{a}x}^{n+1}=\frac{\Delta t}{\mu }\frac{\partial }{\partial z}{E}_{\mathrm{a}y}^{n+\frac{1}{2}}-\frac{\Delta t}{\mu }\frac{\partial }{\partial y}{E}_{\mathrm{a}z}^{n+\frac{1}{2}}\end{array}$$

(107)

$$\begin{array}{c}{h}_{\mathrm{a}y}^{n+1}-\frac{\Delta t^2}{4\mu }\frac{\partial }{\partial z}\frac{1}{ϵ}\frac{\partial }{\partial z}{h}_{\mathrm{a}y}^{n+1}=\frac{\Delta t}{\mu }\frac{\partial }{\partial x}{E}_{\mathrm{a}z}^{n+\frac{1}{2}}-\frac{\Delta t}{\mu }\frac{\partial }{\partial z}{E}_{\mathrm{a}x}^{n+\frac{1}{2}}\end{array}$$

(108)

$$\begin{array}{c}{h}_{\mathrm{a}z}^{n+1}-\frac{\Delta t^2}{4\mu }\frac{\partial }{\partial x}\frac{1}{ϵ}\frac{\partial }{\partial x}{h}_{\mathrm{a}z}^{n+1}=\frac{\Delta t}{\mu }\frac{\partial }{\partial y}{E}_{\mathrm{a}x}^{n+\frac{1}{2}}-\frac{\Delta t}{\mu }\frac{\partial }{\partial x}{E}_{\mathrm{a}y}^{n+\frac{1}{2}}\end{array}$$

(109)

$$\begin{array}{c}{H}_{\mathrm{a}x}^{n+1}=H_{\mathrm{a}x}^n+h_{\mathrm{a}x}^{n+1}\end{array}$$

(110)

$$\begin{array}{c}{H}_{\mathrm{a}y}^{n+1}=H_{\mathrm{a}y}^n+h_{\mathrm{a}y}^{n+1}\end{array}$$

(111)

$$\begin{array}{c}{H}_{\mathrm{a}z}^{n+1}=H_{\mathrm{a}z}^n+h_{\mathrm{a}z}^{n+1}.\end{array}$$

(112)

To yield the physical fields that feature complying divergence and satisfy Gauss’s law, we have the detailed update procedures of leapfrog CDI FDTD method from (91)–(94):

-

First procedure for $\mathbf{E}{}_{\mathrm{c}}$ from $n-\frac{1}{2}$ to $n+\frac{1}{2}$:

$$\begin{array}{c}{h}_{\mathrm{c}x}^n-\frac{\Delta t^2}{4\mu }\frac{\partial }{\partial y}\frac{1}{ϵ}\frac{\partial }{\partial y}{h}_{\mathrm{c}x}^n=H_{\mathrm{c}x}^n\end{array}$$

(113)

$$\begin{array}{c}{h}_{\mathrm{c}y}^n-\frac{\Delta t^2}{4\mu }\frac{\partial }{\partial z}\frac{1}{ϵ}\frac{\partial }{\partial z}{h}_{\mathrm{c}y}^n=H_{\mathrm{c}y}^n\end{array}$$

(114)

$$\begin{array}{c}{h}_{\mathrm{c}z}^n-\frac{\Delta t^2}{4\mu }\frac{\partial }{\partial x}\frac{1}{ϵ}\frac{\partial }{\partial x}{h}_{\mathrm{c}z}^n=H_{\mathrm{c}z}^n\end{array}$$

(115)

$$\begin{array}{c}{E}_{\mathrm{c}x}^{n+\frac{1}{2}}=E_{\mathrm{c}x}^{n-\frac{1}{2}}+\frac{\Delta t}{ϵ}\frac{\partial }{\partial y}{h}_{\mathrm{c}z}^n-\frac{\Delta t}{ϵ}\frac{\partial }{\partial z}{h}_{\mathrm{c}y}^n\end{array}$$

(116)

$$\begin{array}{c}{E}_{\mathrm{c}y}^{n+\frac{1}{2}}=E_{\mathrm{c}y}^{n-\frac{1}{2}}+\frac{\Delta t}{ϵ}\frac{\partial }{\partial z}{h}_{\mathrm{c}x}^n-\frac{\Delta t}{ϵ}\frac{\partial }{\partial x}{h}_{\mathrm{c}z}^n\end{array}$$

(117)

$$\begin{array}{c}{E}_{\mathrm{c}z}^{n+\frac{1}{2}}=E_{\mathrm{c}z}^{n-\frac{1}{2}}+\frac{\Delta t}{ϵ}\frac{\partial }{\partial x}{h}_{\mathrm{c}y}^n-\frac{\Delta t}{ϵ}\frac{\partial }{\partial y}{h}_{\mathrm{c}x}^n\end{array}$$

(118)

-

Second procedure for $\mathbf{H}{}_{\mathrm{c}}$ from n to $n+1$:

$$\begin{array}{c}{e}_{\mathrm{c}x}^{n+\frac{1}{2}}-\frac{\Delta t^2}{4ϵ}\frac{\partial }{\partial y}\frac{1}{\mu }\frac{\partial }{\partial y}{e}_{\mathrm{c}x}^{n+\frac{1}{2}}=E_{\mathrm{c}x}^{n+\frac{1}{2}}\end{array}$$

(119)

$$\begin{array}{c}{e}_{\mathrm{c}y}^{n+\frac{1}{2}}-\frac{\Delta t^2}{4ϵ}\frac{\partial }{\partial z}\frac{1}{\mu }\frac{\partial }{\partial z}{e}_{\mathrm{c}y}^{n+\frac{1}{2}}=E_{\mathrm{c}y}^{n+\frac{1}{2}}\end{array}$$

(120)

$$\begin{array}{c}{e}_{\mathrm{c}z}^{n+\frac{1}{2}}-\frac{\Delta t^2}{4ϵ}\frac{\partial }{\partial x}\frac{1}{\mu }\frac{\partial }{\partial x}{e}_{\mathrm{c}z}^{n+\frac{1}{2}}=E_{\mathrm{c}z}^{n+\frac{1}{2}}\end{array}$$

(121)

$$\begin{array}{c}{H}_{\mathrm{c}x}^{n+1}=H_{\mathrm{c}x}^n+\frac{\Delta t}{\mu }\frac{\partial }{\partial z}{e}_{\mathrm{c}y}^{n+\frac{1}{2}}-\frac{\Delta t}{\mu }\frac{\partial }{\partial y}{e}_{\mathrm{c}z}^{n+\frac{1}{2}}\end{array}$$

(122)

$$\begin{array}{c}{H}_{\mathrm{c}y}^{n+1}=H_{\mathrm{c}y}^n+\frac{\Delta t}{\mu }\frac{\partial }{\partial x}{e}_{\mathrm{c}z}^{n+\frac{1}{2}}-\frac{\Delta t}{\mu }\frac{\partial }{\partial z}{e}_{\mathrm{c}x}^{n+\frac{1}{2}}\end{array}$$

(123)

$$\begin{array}{c}{H}_{\mathrm{c}z}^{n+1}=H_{\mathrm{c}z}^n+\frac{\Delta t}{\mu }\frac{\partial }{\partial y}{e}_{\mathrm{c}x}^{n+\frac{1}{2}}-\frac{\Delta t}{\mu }\frac{\partial }{\partial x}{e}_{\mathrm{c}y}^{n+\frac{1}{2}}.\end{array}$$

(124)

Based on the detailed update procedures above, we first compare various unconditionally stable leapfrog implicit FDTD schemes in

Table 1. Comparison of Unconditionally Stable Leapfrog Implicit FDTD Schemes.

Implicit FDTD Scheme		Leapfrog ADI FDTD	Leapfrog FADI FDTD	Leapfrog CDI FDTD
Implicit update equations		(36) & (43), or (95)–(100)	(44) & (46), or (101)–(103) & (107)–(109)	(91) & (93), or (113)–(115) & (119)–(121)
Explicit update equations		–	(45) & (47), or (104)–(106) & (110)–(112)	(92) & (94), or (116)–(118) & (122)–(124)
Implicit	M/D	24	12	0
update	A/S	48	18	0
Explicit	M/D	0	0	12
update	A/S	0	6	24
Total	M/D	24	12	12
	A/S	48	24	24
	M/D+A/S	72	36	36
Temporal accuracy		Second-order	Second-order	Second-order
Complying divergence		No, violating Gauss’s law	No, violating Gauss’s law	Yes, satisfying Gauss’s law

, which include the leapfrog ADI, FADI and CDI FDTD methods. The table lists the flops count for one complete step involving the pertaining implicit and explicit update equations using second-order spatial central-differencing on Yee cells. The flops consist of the operations of multiplications/divisions (M/D) and additions/subtractions (A/S) at the RHS of update equations (with the multiplicative factors precomputed). From Table 1, one can see that the conventional leapfrog ADI FDTD method still requires substantial flops, while the leapfrog FADI and CDI FDTD methods require merely half the flops. These latter leapfrog schemes have been formulated in the fundamental forms, which have omitted the cumbersome second-order matrix operators at their RHS. Although all the presented leapfrog schemes are of second-order temporal accuracy with respect to time-staggered fields, only the leapfrog CDI FDTD method has complying divergence satisfying Gauss’s law.

For completeness, we also compare various unconditionally stable time-collocated implicit FDTD schemes in

Table 2. Comparison of Unconditionally Stable Time-Collocated Implicit FDTD Schemes.

Implicit FDTD Scheme		ADI FDTD	FADI FDTD Implicit E-E or E-H	LOD FDTD	FLOD-CD or CDI FDTD
Implicit update equations		(6) & (8)	E-E: (16) & (19) E-H: (23) & (26)	(50) & (52)	(58) & (61)
Explicit update equations		(7) & (9)	E-E: (15), (17), (18) & (20) E-H: (22), (24), (25) & (27)	(51) & (53)	(59), (60), (62) & (63)
Implicit	M/D	24	6	18	6
update	A/S	48	12	36	12
Explicit	M/D	12	6	6	6
update	A/S	24	18	24	18
Total	M/D	36	12	24	12
	A/S	72	30	60	30
	M/D+A/S	108	42	84	42
Temporal accuracy		Second-order	Second-order	First-order, (stable for NUTS)	Second-order with (64) and (65)
Complying divergence		No, violating Gauss’s law	No, violating Gauss’s law	No, violating Gauss’s law	Yes, satisfying Gauss’s law with (64) and (65)

, which include the ADI, FADI, LOD and FLOD-CD or CDI FDTD methods. The table lists the flops count for one complete step with the pertaining implicit and explicit update equations given in matrix forms that are readily expanded into detailed component equations as before. From Table 2, one can see that the conventional ADI and LOD FDTD methods require many more flops than the fundamental schemes. Note that their flops are slightly more for general inhomogeneous media as compared to the homogeneous case in [30]. Meanwhile, the FADI and FLOD-CD or CDI FDTD methods can reduce the flops considerably from the conventional ADI and LOD FDTD methods. Upon further comparison with Table 1, it is evident that the time-collocated fundamental schemes require fewer flops than the conventional leapfrog ADI FDTD method, while the leapfrog fundamental schemes for ADI and CDI FDTD methods constitute the most efficient implicit FDTD schemes so far. In fact, the time-collocated fundamental schemes may still be advantageous since with comparable flops (only 6 more additions), they can provide both electric and magnetic fields completely at every time instant and deal with PEC/PMC media without field leakage [52,53]. Considering the temporal accuracy, all the presented time-collocated schemes are of second-order except the LOD FDTD method. This method is only first-order accurate but may remain stable for non-uniform time step (NUTS) cases while other methods fail [73]. Moreover, it is readily adapted to feature complying divergence as in FLOD-CD or CDI FDTD method satisfying Gauss’s law (while other methods violate it).

To further illustrate the divergence property of various implicit FDTD methods, let us consider a simple cavity of PEC walls meshed with $50\times 50\times 9$ uniform grids. With the source being impressed at the cavity center, we implement the methods for various $\mathrm{CFLN}=\Delta t/\Delta t_{\mathrm{CFL}}$ ($\Delta t_{\mathrm{CFL}}$ is the CFL limit time step size) and calculate the normalized numerical divergence of $\mathbf{D}$ in dB. Figure 1 shows the normalized numerical divergence of $\mathbf{D}$ in the source plane using (a) time-collocated/leapfrog ADI FDTD for $\mathrm{CFLN}=1$; (b) leapfrog CDI FDTD for $\mathrm{CFLN}=1$; (c) time-collocated/leapfrog ADI FDTD for $\mathrm{CFLN}=4$; (d) leapfrog CDI FDTD for $\mathrm{CFLN}=4$. From Figure 1a,b for $\mathrm{CFLN}=1$, we see that the ADI FDTD method is not divergence-free even outside the source point. On the other hand, the leapfrog CDI FDTD method has divergence only at the source point while its divergence outside is mostly of numerical noise level (∼−200 dB). For $\mathrm{CFLN}=4$ in Figure 1c,d, one can see that the ADI FDTD method has larger spread of non-zero divergence outside the source point, whereas the leapfrog CDI FDTD method still has complying divergence everywhere.

Thus far, the methods discussed in this paper have found applications mostly for electromagnetic problems. They should be readily extendable to handle multi-physics problems such as thermal/heat conduction, quantum mechanics, circuits, etc. These can be exemplified by our previous works using time-collocated fundamental schemes of ADI and LOD FDTD methods [73,74,75]. Alternatively, the multi-physics problems can be treated using leapfrog fundamental schemes of ADI and CDI FDTD methods to exploit their advantages further.

5. Conclusions

In this paper, the formulations from time-collocated to leapfrog fundamental schemes have been presented for ADI and CDI FDTD methods. For the ADI FDTD method, the time-collocated fundamental schemes have been implemented using implicit E-E and E-H update procedures. Their update equations comprise simple and concise RHS that contain only the intrinsic matrix operator partitions, but without their complicated products that exist in the conventional form. From the fundamental implicit E-H scheme, the leapfrog ADI FDTD method has been formulated in conventional form, whose RHS still involve the products of matrix operator partitions. Using auxiliary variables, the RHS have been simplified into the leapfrog fundamental scheme for ADI FDTD method with reduced operations and improved efficiency. For the CDI FDTD method, the time-collocated fundamental scheme has been presented based on the LOD FDTD method with complying divergence. The formulations from time-collocated to leapfrog schemes have been provided, which result in the leapfrog fundamental scheme for CDI FDTD method. Based on their equations in fundamental forms, further insights have been given into the relations of leapfrog fundamental schemes for ADI and CDI FDTD methods. Comparisons among various implicit FDTD schemes have been discussed including flops count, temporal accuracy and divergence property. The time-collocated fundamental schemes require considerably fewer flops than all conventional ADI, LOD and leapfrog ADI FDTD methods, while the leapfrog fundamental schemes for ADI and CDI FDTD methods constitute the most efficient implicit FDTD schemes to date. While providing a unified formulation for simplification, the fundamental schemes feature many advantages including unconditional stability, simplicity and efficient time-collocated/leapfrog update, as well as complying divergence (for FLOD-CD or CDI FDTD). They are very promising and well-suited for further developments and applications in electromagnetic computations and simulations.