Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group
Journal Home About Issues in Progress Current Issue All Issues Feature Issues

Auxiliary differential equation (ADE) method based complying-divergence implicit FDTD method for simulating the general dispersive anisotropic material

Guoda Xie, Guilin Hou, Naixing Feng, Kaihong Song, Ming Fang, Yingsong Li, Xianliang Wu, and Zhixiang Huang

Open Access Open Access

Abstract

The preceding works introduced the leapfrog complying divergence implicit finite-difference time-domain (CDI-FDTD) method, which exhibits high accuracy and unconditional stability. In this study, the method is reformulated to simulate general electrically anisotropic and dispersive media. The auxiliary differential equation (ADE) method is employed to solve the equivalent polarization currents, which are then integrated into the CDI-FDTD method. The iterative formulae are presented, and the calculation method is similar to that of the traditional CDI-FDTD method. Additionally, the Von Neumann method is utilized to analyze the unconditional stability of the proposed method. To evaluate the performance of the proposed method, three numerical cases are conducted. These include calculating the transmission and reflection coefficients of a monolayer graphene sheet and a monolayer magnetized plasma, as well as the scattering properties of a cubic block plasma. The numerical results obtained by the proposed method demonstrate its accuracy and efficiency in simulating general anisotropic dispersive media, compared to both the analytical method and the traditional FDTD method.

© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

The study of the electromagnetic characteristics of anisotropic materials holds vital significance for numerous domains, including electronics, novel electronic nanodevices, optical materials, navigation, communication systems and biology. For instance, researchers incorporate anisotropic effects into partial element equivalent circuit (PEEC) models to address various electromagnetic compatibility issues [1]. The design of graphene frequency selective surfaces (GFSS) [2] and optical and microwave devices [3] also require an understanding of anisotropic materials. Additionally, some scholars concentrate on developing effective techniques for analyzing the electromagnetic properties of multilayered microstrip transmission lines containing anisotropic dielectric layers [4], anisotropic dielectric objects embedded in layered arbitrary anisotropic media [5], and studying the phenomenon of reflection of plane electrostatic waves on a plane interface between a magnetized plasma and an insulator [6]. Therefore, investigating the electromagnetic properties of anisotropic materials is not only critical for practical problem-solving but also an essential focus in the field of computational electromagnetics. The time-domain discretization of Maxwell's differential equations enables the acquisition of wideband responses through time-domain computations, rendering the approach simple and practical. The finite-difference time-domain (FDTD) technique has been extensively utilized for modeling dispersive media and structures made of a diverse range of isotropic and anisotropic materials, which, in turn, aids in exploring the propagation and scattering of electromagnetic waves in such materials. Notably, this method has been applied to study the behavior of electromagnetic waves in magnetized plasma [7–10], metasurfaces consisting of asymmetric scatterers [11], and invisibility materials [12]. However, the conventional explicit FDTD approaches are restricted by the Courant-Friedrichs-Lewy (CFL) time condition, which imposes a strict limit on the range of admissible time steps contingent on the spatial grid size [13–17]. Consequently, the computational expense of the numerical simulation is significantly amplified when employing this method to model electromagnetic problems featuring fine structures. As a remedy, alternative FDTD methods, which are unconditionally stable, such as the alternating direction implicit finite-difference time-domain (ADI-FDTD), and locally one-dimensional FDTD (LOD-FDTD) methods, have been advanced to surmount the CFL limitation and have been effectively utilized to emulate anisotropic dispersive models [18,19].

However, both of these techniques necessitate bifurcating one time step into two sub-time steps during the iterative process, which results in poor computational efficiency [20,21]. Recently, the one-step leapfrog ADI-FDTD (LADI-FDTD) approach has garnered substantial interest. This technique stems from the conventional ADI-FDTD method, but incorporates a one-step update mechanism analogous to the explicit FDTD method for both electric and magnetic fields, while preserving the same numerical precision as the conventional ADI-FDTD method [22–24]. As a consequence, the LADI-FDTD method has attracted persistent interest and development, and has been successfully employed to emulate the interplay between electromagnetic waves and dispersive medium. The studies presented in [25–27] and [28] respectively introduce LADI-FDTD methods predicated on the current density and electric field (J-E) constitutive relation and the electric flux and electric field (D-E) constitutive relation to manage frequency-dependent dispersive medium, and demonstrate the precise and efficient simulation of the LADI-FDTD method for dispersive medium. Nevertheless, research on the LADI-FDTD approach has revealed that it fails to satisfy the divergence property and contravenes Gauss's law [29].

A new unconditionally stable the leapfrog complying-divergence implicit finite-difference time-domain (CDI-FDTD) approach has recently been introduced in [30–32], exhibiting superior numerical accuracy compared to the LADI-FDTD method [33,34]. Moreover, a high-order CPML grounded on the CDI-FDTD approach has been documented for numerical implementation, resolving electromagnetic computational difficulties in open boundary condition [35]. Nevertheless, to the best of the author's knowledge, the potence of the CDI-FDTD method in simulating anisotropic dispersive media is limited by the paucity of research in this area, hence curtailing the applicability of this method in computational problems.

Additionally, in simulating dispersive medium, a crucial concern pertains to the conversion of the dispersion model function that characterizes the material's electromagnetic properties from the frequency domain to the time domain, as governed by the FDTD framework. Various methods exist for effecting this conversion, including the recursive convolution method (RC), the auxiliary differential equation method (ADE), the Z-transform method [36–40], etc. The ADE method stands out as one that is both easily implementable and highly accurate, capable of effectively treating linear and nonlinear dispersive medium [41]. In a previous study, the ADE-FDTD method was employed in [42] to simulate the propagation of electromagnetic waves in anisotropic dispersive medium starting from the D-E constitutive relationship. Specifically, the researchers applied this method to model the behavior of electromagnetic waves in one-dimensional magnetized plasma, and subsequently extended their findings to two-dimensional electromagnetic simulations. Moreover, several studies have explored the gyroscopic properties of graphene in the presence of a static magnetic field at microwave and optical frequencies, and have developed highly efficient techniques for rapidly simulating the interaction between electromagnetic waves and graphene structures using the ADE-FDTD method [43–47].

Drawing on the ADE method's demonstrated efficacy in treating anisotropic dispersive medium and the superior performance of the CDI-FDTD method, which offers unconditionally stable and accurate results, we have developed a novel computational approach that combines the CDI-FDTD method with the ADE method to simulate general anisotropic dispersive medium. Thereafter, the study investigates the numerical stability of the proposed frequency dependent ADE-CDI-FDTD method. Finally, three different numerical models are simulated to verify the accuracy and efficiency of this method.

The paper is structured as follows: In Section 2, we derive the numerical iteration formula for the frequency-dependent ADE-CDI-FDTD method used to simulate anisotropic media. In Section 3, a numerical stability analysis is presented for simulating graphene sheets subject to a static magnetic field. Finally, in Section 4, we provide results from three numerical cases that demonstrate the efficiency and accuracy of our proposed method.

2. Method

2.1 Derivation of the CDI-FDTD method for anisotropic dispersive medium

For a general anisotropic dispersive medium, its permittivity tensor can be represented as:

$$\widetilde \varepsilon (\omega ) = {\varepsilon _0} + {\varepsilon _0}\widetilde \chi (\omega ),$$
where, ε0 is the relative permittivity in vacuum, and the susceptibility function $\widetilde \chi$ is given by:
$$\widetilde \chi (\omega ) = \left[ {\begin{array}{{ccc}} {{\varepsilon_{xx}}(\omega )}&{{\varepsilon_{xy}}(\omega )}&{{\varepsilon_{xz}}(\omega )}\\ {{\varepsilon_{yx}}(\omega )}&{{\varepsilon_{yy}}(\omega )}&{{\varepsilon_{yz}}(\omega )}\\ {{\varepsilon_{zx}}(\omega )}&{{\varepsilon_{zy}}(\omega )}&{{\varepsilon_{zz}}(\omega )} \end{array}} \right],$$

Simultaneously, the Maxwell's equations in frequency domain are:

$$j\omega \widetilde \varepsilon \textrm{(}\omega \textrm{)}{\boldsymbol{E}}\textrm{(}\omega \textrm{)} = \nabla \times {\boldsymbol{H}}\textrm{(}\omega \textrm{)},$$
$$j\omega {\mu _0}{\boldsymbol{H}}\textrm{(}\omega \textrm{)} = \nabla \times {\boldsymbol{E}}\textrm{(}\omega \textrm{)}.$$

Substituting Eq. (1) into Eq. (3), we get:

$$j\omega \widetilde \varepsilon \textrm{(}\omega \textrm{)}{\boldsymbol{E}}\textrm{(}\omega \textrm{)} = j\omega {\varepsilon _0}{\boldsymbol{E}}\textrm{(}\omega \textrm{) + }{\boldsymbol{J}}\textrm{(}\omega \textrm{),}$$
where ${\boldsymbol{J}}\textrm{(}\omega \textrm{)} = j\omega {\varepsilon _0}\widetilde \chi \textrm{(}\omega \textrm{)}{\boldsymbol{E}}\textrm{(}\omega \textrm{)}$, ${\boldsymbol{E}} = {[{{E_x},{E_y},{E_z}} ]^T}$, ${\boldsymbol H} = {[{{H_x},{H_y},{H_z}} ]^T}$, ${\boldsymbol{J}} = {[{{J_x},{J_y},{J_z}} ]^T}$. According to Eq. (2), the components of current density J in x, y, and z directions can be written as:
$${J_x}(\omega ) = j\omega {\varepsilon _0}[{{\varepsilon_{xx}}{E_x}(\omega ) + {\varepsilon_{xy}}{E_y}(\omega ) + {\varepsilon_{xz}}{E_z}(\omega )} ],$$
$${J_y}(\omega ) = j\omega {\varepsilon _0}[{{\varepsilon_{yx}}{E_x}(\omega ) + {\varepsilon_{yy}}{E_y}(\omega ) + {\varepsilon_{yz}}{E_z}(\omega )} ],$$
$${J_z}(\omega ) = j\omega {\varepsilon _0}[{{\varepsilon_{zx}}{E_x}(\omega ) + {\varepsilon_{zy}}{E_y}(\omega ) + {\varepsilon_{zz}}{E_z}(\omega )} ].$$

According to Eq. (5), using the Fourier transform theory to transform Eqs. (3) and (4) from frequency domain to time domain, we get:

$$\frac{{\partial {\boldsymbol{E}}}}{{\partial t}} + \frac{1}{{{\varepsilon _0}}}{\boldsymbol{J}} = ({{\boldsymbol{A}}_{12}} + {{\boldsymbol{B}}_{12}}){\boldsymbol{H}}, $$
$$\frac{{\partial {\boldsymbol{H}}}}{{\partial t}} = ({{\boldsymbol{A}}_{21}} + {{\boldsymbol{B}}_{21}}){\boldsymbol{E}}, $$
where the coefficient matrices A12, B12, A21, and B21 are given by:
$${{\boldsymbol{A}}_{12}} = \left[ {\begin{array}{{ccc}} 0&0&{\frac{1}{{{\varepsilon _0}}}\frac{\partial }{{\partial y}}} \\ {\frac{1}{{{\varepsilon _0}}}\frac{\partial }{{\partial z}}}&0&0 \\ 0&{\frac{1}{{{\varepsilon _0}}}\frac{\partial }{{\partial x}}}&0 \end{array}} \right],\,{{\boldsymbol{B}}_{12}} = \left[ {\begin{array}{{ccc}} 0&{\frac{{ - 1}}{{{\varepsilon _0}}}\frac{\partial }{{\partial z}}}&0 \\ 0&0&{\frac{{ - 1}}{{{\varepsilon _0}}}\frac{\partial }{{\partial x}}} \\ {\frac{{ - 1}}{{{\varepsilon _0}}}\frac{\partial }{{\partial y}}}&0&0 \end{array}} \right]$$
$${{\boldsymbol{A}}_{21}} = \left[ {\begin{array}{{ccc}} 0&{\frac{1}{{{\mu _0}}}\frac{\partial }{{\partial z}}}&0 \\ 0&0&{\frac{1}{{{\mu _0}}}\frac{\partial }{{\partial x}}} \\ {\frac{1}{{{\mu _0}}}\frac{\partial }{{\partial y}}}&0&0 \end{array}} \right],\,{{\boldsymbol{B}}_{21}} = \left[ {\begin{array}{{ccc}} 0&0&{\frac{{ - 1}}{{{\mu _0}}}\frac{\partial }{{\partial y}}} \\ {\frac{{ - 1}}{{{\mu _0}}}\frac{\partial }{{\partial z}}}&0&0 \\ 0&{\frac{{ - 1}}{{{\mu_0}}}\frac{\partial }{{\partial x}}}&0 \end{array}} \right]$$

The numerical calculation formula obtained using the locally one-dimensional (LOD) scheme based on the splitting idea to handle the dispersive medium with formulas (7) and (8) can be expressed in a compact matrix form as follows:

$$({{\boldsymbol{I}}_{6 \times 6}} - \frac{{\Delta t}}{2}{\boldsymbol{A}}){{\boldsymbol{u}}^{n + 1/2}} = ({{\boldsymbol{I}}_{6 \times 6}} + \frac{{\Delta t}}{2}{\boldsymbol{A}}){{\boldsymbol{u}}^n} - \Delta t{{\boldsymbol{S}}^n}$$
$$({{\boldsymbol{I}}_{6 \times 6}} - \frac{{\Delta t}}{2}{\boldsymbol{B}}){{\boldsymbol{u}}^{n + 1}} = ({{\boldsymbol{I}}_{6 \times 6}} + \frac{{\Delta t}}{2}{\boldsymbol{B}}){{\boldsymbol{u}}^{n + 1/2}}, $$
where In × n is the n-dimensional identity matrix, ${\boldsymbol u} = {[{{\boldsymbol E},{\boldsymbol H}} ]^T}$, ${\boldsymbol{S}} = {\left[ {\frac{{{J_x}}}{{{\varepsilon_0}}},\frac{{{J_y}}}{{{\varepsilon_0}}},\frac{{{J_z}}}{{{\varepsilon_0}}},0,0,0} \right]^T}$, and the coefficient matrices A, B are given by:
$${\boldsymbol{A}} = \left[ {\begin{array}{{cc}} {{\mathbf{0}_{3 \times 3}}}&{{{\boldsymbol{A}}_{12}}} \\ {{{\boldsymbol{A}}_{21}}}&{{\mathbf{0}_{3 \times 3}}} \end{array}} \right],\,{\boldsymbol{B}} = \left[ {\begin{array}{{cc}} {{\mathbf{0}_{3 \times 3}}}&{{{\boldsymbol{B}}_{12}}} \\ {{{\boldsymbol{B}}_{21}}}&{{\mathbf{0}_{3 \times 3}}} \end{array}} \right].$$

To endow the LOD-FDTD method with both second-order temporal accuracy and a consistent divergence feature, it is necessary to apply the following processing steps separately at n = 0 and at the output [30]:

The first step is input processing:

$$({{\boldsymbol{I}}_{6 \times 6}} - \frac{{\Delta t}}{2}{\boldsymbol{B}}){{\boldsymbol{u}}^0} = {\boldsymbol{u}}_c^0, $$

The second step is output processing:

$${\boldsymbol{u}}_c^{n + 1} = ({{\boldsymbol{I}}_{6 \times 6}} - \frac{{\Delta t}}{2}{\boldsymbol{B}}){{\boldsymbol{u}}^{n + 1}}. $$
where ${{\boldsymbol u}_c} = {[{{{\boldsymbol{E}}_c},{{\boldsymbol{H}}_c}} ]^T}$, the subscript ‘c’ represents complying divergence [32]. Eq. (15) yields the result of obeying divergence and satisfying Gauss law at n + 1 time step. Similarly, the result of uc in n + 1/2 step is obtained as:
$${\boldsymbol{u}}_c^{n + 1/2} = ({{\boldsymbol{I}}_{6 \times 6}} - \frac{{\Delta t}}{2}{\boldsymbol{A}}){{\boldsymbol{u}}^{n + 1/2}}, $$

Substituting Eq. (16) into Eq. (11), we get:

$${\boldsymbol{u}}_c^{n + 1/2} = ({{\boldsymbol{I}}_{6 \times 6}} + \frac{{\Delta t}}{2}{\boldsymbol{A}}){{\boldsymbol{u}}^n} - \Delta t{{\boldsymbol{S}}^n}. $$

Taking Eq. (15) at one time step backward, its inverse relation is derived:

$${{\boldsymbol{u}}^n} = {({{\boldsymbol{I}}_{6 \times 6}} - \frac{{\Delta t}}{2}{\boldsymbol{B}})^{ - 1}}{\boldsymbol{u}}_c^n, $$

Then substituting Eq. (18) into Eq. (17), we get:

$$\begin{aligned} {u}_c^{\boldsymbol{n} + 1/2} &= ({{\boldsymbol{I}}_{6 \times 6}} + \frac{{\Delta t}}{2}{\boldsymbol{A}}){({{\boldsymbol{I}}_{6 \times 6}} - \frac{{\Delta t}}{2}{\boldsymbol{B}})^{ - 1}}{\boldsymbol{u}}_c^n - \Delta t{{\boldsymbol{S}}^n}\\& = ({{\boldsymbol I}_{6 \times 6}} + \frac{{\Delta t}}{2}{\boldsymbol{A}})({{\boldsymbol{I}}_{6 \times 6}} + \frac{{\Delta t}}{2}{\boldsymbol{B}}){\left[ {({{I}_{6 \times 6}} - \frac{{\Delta t}}{2}{\boldsymbol{B}})({{\boldsymbol{I}}_{6 \times 6}} + \frac{{\Delta t}}{2}{\boldsymbol{B}})} \right]^{ - 1}}{\boldsymbol{u}}_c^n - \Delta t{{\boldsymbol{S}}^n} \end{aligned}.$$

According to Eq. (19), the iterative formulas of the electric and magnetic fields at n + 1/2 moments are further derived as:

$$\begin{array}{c} \left[ {\begin{array}{{c}} {{\boldsymbol E}_c^{n + 1/2}}\\ {{\boldsymbol H}_c^{n + 1/2}} \end{array}} \right] = \textrm{(}{{\boldsymbol I}_{6 \times 6}} + \frac{{\Delta t}}{2}\left[ {\begin{array}{{cc}} {{{\mathbf 0}_{3 \times 3}}}&{{{\boldsymbol A}_{12}}}\\ {{{\boldsymbol A}_{21}}}&{{{\mathbf 0}_{3 \times 3}}} \end{array}} \right] + \frac{{\Delta t}}{2}\left[ {\begin{array}{{cc}} {{{\mathbf 0}_{3 \times 3}}}&{{{B}_{12}}}\\ {{{\boldsymbol B}_{21}}}&{{{\mathbf 0}_{3 \times 3}}} \end{array}} \right] + \frac{{\Delta {t^2}}}{4}\left[ {\begin{array}{{cc}} {{{\boldsymbol A}_{12}}{{\boldsymbol B}_{21}}}&{{{\mathbf 0}_{3 \times 3}}}\\ {{{\mathbf 0}_{3 \times 3}}}&{{{\boldsymbol A}_{21}}{{\boldsymbol B}_{12}}} \end{array}} \right])\\ {\textrm{(}{{\boldsymbol I}_{6 \times 6}} - \frac{{\Delta {t^2}}}{4}\left[ {\begin{array}{{cc}} {{{\boldsymbol B}_{12}}{{\boldsymbol B}_{21}}}&{{{\mathbf 0}_{3 \times 3}}}\\ {{{\mathbf 0}_{3 \times 3}}}&{{{\boldsymbol B}_{21}}{{\boldsymbol B}_{12}}} \end{array}} \right])^{ - 1}}\left[ {\begin{array}{{c}} {{\boldsymbol E}_c^n}\\ {{\boldsymbol H}_c^n} \end{array}} \right] - \frac{{\Delta t}}{{{\varepsilon _0}}}\left[ {\begin{array}{{c}} {{\boldsymbol J}_c^n}\\ {{{\mathbf 0}_{1 \times 3}}} \end{array}} \right] \end{array}, $$
where ${{\boldsymbol E}_c} = {[{{E_{cx}},{E_{cy}},{E_{cz}}} ]^T}$, ${{\boldsymbol H}_c} = {[{{H_{cx}},{H_{cy}},{H_{cz}}} ]^T}$, ${{\boldsymbol J}_c} = {[{{J_{cx}},{J_{cy}},{J_{cz}}} ]^T}$.

Therefore, the iterative formula for the electric field Ec at n + 1/2 moments is:

$$\begin{aligned} {\boldsymbol E}_c^{n + 1/2} &= ({{\boldsymbol I}_{3 \times 3}} + \frac{{\Delta {t^2}}}{4}{{\boldsymbol A}_{12}}{{\boldsymbol B}_{21}}){({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol B}_{12}}{{\boldsymbol B}_{21}})^{ - 1}}{\boldsymbol E}_c^n\\& + \frac{{\Delta t}}{2}({{\boldsymbol A}_{12}} + {{\boldsymbol B}_{12}}){({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol B}_{21}}{{\boldsymbol B}_{12}})^{ - 1}}{\boldsymbol H}_c^n - \frac{{\Delta t}}{{{\varepsilon _0}}}{\boldsymbol J}_c^n \end{aligned}, $$

Specifically, the derivation of the electric field Ec at time interval n + 1/2 → n + 1 (Eq. 12) is consistent with that in free space. The literature [32] provides an iterative formula for the electric field Ec at n − 1/2 time step:

$$\begin{array}{c} {\boldsymbol E}_c^{n - 1/2} = ({{\boldsymbol I}_{3 \times 3}} + \frac{{\Delta {t^2}}}{4}{{\boldsymbol A}_{12}}{{\boldsymbol B}_{21}}){({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol B}_{12}}{{\boldsymbol B}_{21}})^{ - 1}}{\boldsymbol E}_c^n\\ - \frac{{\Delta t}}{2}({{\boldsymbol A}_{12}} + {{\boldsymbol B}_{12}}){({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol B}_{21}}{{\boldsymbol B}_{12}})^{ - 1}}{\boldsymbol H}_c^n \end{array}, $$

Subtracting Eq. (21) from Eq. (22), a time step iterative formula for the electric field intensity is obtained:

$${\boldsymbol E}_c^{n + 1/2} - {\boldsymbol E}_c^{n - 1/2} = \Delta t({{\boldsymbol A}_{12}} + {{\boldsymbol B}_{12}}){({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol B}_{21}}{{\boldsymbol B}_{12}})^{ - 1}}{\boldsymbol H}_c^n - \frac{{\Delta t}}{{{\varepsilon _0}}}{\boldsymbol J}_c^n. $$

Similarly, the numerical iterative formula of magnetic field is derived as:

$${\boldsymbol H}_c^{n + 1} - {\boldsymbol H}_c^n = \Delta t({{\boldsymbol A}_{21}} + {{\boldsymbol B}_{21}}){({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol A}_{21}}{{\boldsymbol A}_{12}})^{ - 1}}{\boldsymbol E}_c^{n + 1/2}. $$

Then the auxiliary variables hc and ec are defined as ${{\boldsymbol h}_c} = {({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol B}_{21}}{{\boldsymbol B}_{12}})^{ - 1}}{{\boldsymbol H}_c}$ and ${{\boldsymbol e}_c} = {({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol A}_{21}}{{\boldsymbol A}_{12}})^{ - 1}}{{\boldsymbol E}_c}$, respectively. Substituting them into Eqs. (23) and (24), we get:

$$({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol B}_{21}}{{B}_{12}}){\boldsymbol h}_c^n = {\boldsymbol H}_c^n$$
$${\boldsymbol E}_c^{n + 1/2} = {\boldsymbol E}_c^{n - 1/2} + \Delta t({{\boldsymbol A}_{12}} + {{\boldsymbol B}_{12}}){\boldsymbol h}_c^n - \frac{{\Delta t}}{{{\varepsilon _0}}}{\boldsymbol J}_c^n$$
$$({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol A}_{12}}{{\boldsymbol A}_{21}}){\boldsymbol e}_c^{n + 1/2} = {\boldsymbol E}_c^{n + 1/2}$$
$${\boldsymbol H}_c^{n + 1} = {\boldsymbol H}_c^n + \Delta t({{\boldsymbol A}_{21}} + {{\boldsymbol B}_{21}}){\boldsymbol e}_c^{n + 1/2}. $$

Finally, Eqs. (25)–(28) represent the updating formulae for the frequency-dependent CDI-FDTD method used to simulate general dispersive media. Here, Ec and Hc represent the output electric field and magnetic field, respectively. However, the updating formulas for the current density Jc are not provided in the preceding sections. To address this, this paper presents the updating formulas for Jc for two types of anisotropic dispersive media below:

  • A. Magnetized plasma

When the external magnetic field bias is aligned with the z-axis direction, the electrons present in the plasma will rotate in accordance with the direction of the magnetic field. This rotation gives rise to an anisotropic feature of the plasma. Thus, the susceptibility function in Eq. (2) is represented as:

$$\widetilde \chi (\omega ) = \left[ {\begin{array}{{ccc}} {{\varepsilon_d}(\omega )}&{j{\varepsilon_g}(\omega )}&0\\ { - j{\varepsilon_g}(\omega )}&{{\varepsilon_d}(\omega )}&0\\ 0&0&{{\varepsilon_{zz}}(\omega )} \end{array}} \right], $$
where ${\varepsilon _d}(\omega ) = \frac{{\omega _p^2 + \omega _p^2{v_c}/j\omega }}{{(v_c^2 + \omega _b^2) + 2j\omega {v_c} + {{(j\omega )}^2}}}$, ${\varepsilon _g}(\omega ) = \frac{{ - \omega _p^2{\omega _b}/j\omega }}{{(v_c^2 + \omega _b^2) + 2j\omega {v_c} + {{(j\omega )}^2}}}$, ${\varepsilon _{zz}}(\omega ) = \frac{{ - \omega _p^2}}{{\omega (\omega - j{v_c})}}$. Thus, substituting Eq. (29) into ${{\boldsymbol J}_c}\textrm{(}\omega \textrm{)} = j\omega {\varepsilon _0}\widetilde \chi \textrm{(}\omega \textrm{)}{{\boldsymbol E}_c}\textrm{(}\omega \textrm{)}$, we get:
$${{\boldsymbol J}_c} = \left[ {\begin{array}{{ccc}} {\frac{{j\omega {\varepsilon_0}\omega_p^2 + {\varepsilon_0}\omega_p^2{v_c}}}{{({v_c} + \omega_b^2) + 2j\omega {v_c} + {{(j\omega )}^2}}}}&{\frac{{ - {\varepsilon_0}\omega_p^2{\omega_b}}}{{({v_c} + \omega_b^2) + 2j\omega {v_c} + {{(j\omega )}^2}}}}&0\\ {\frac{{{\varepsilon_0}\omega_p^2{\omega_b}}}{{({v_c} + \omega_b^2) + 2j\omega {v_c} + {{(j\omega )}^2}}}}&{\frac{{j\omega {\varepsilon_0}\omega_p^2 + {\varepsilon_0}\omega_p^2{v_c}}}{{({v_c} + \omega_b^2) + 2j\omega {v_c} + {{(j\omega )}^2}}}}&0\\ 0&0&{\frac{{ - j\omega {\varepsilon_0}\omega_p^2}}{{\omega (\omega - j{v_c})}}} \end{array}} \right]\left[ {\begin{array}{{c}} {{E_{cx}}}\\ {{E_{cy}}}\\ {{E_{cz}}} \end{array}} \right]. $$

Combining the ADE method, the update equations of the current density Jc can be effectively achieved. Firstly, according to Eq. (30), the current density Jc can be defined as follow:

$${{\boldsymbol J}_c} = \left[ {\begin{array}{{c}} {{J_{xx}}(\omega ) + {J_{xy}}(\omega ) + {J_{xz}}(\omega )}\\ {{J_{yx}}(\omega ) + {J_{yy}}(\omega ) + {J_{yz}}(\omega )}\\ {{J_{zx}}(\omega ) + {J_{zy}}(\omega ) + {J_{zz}}(\omega )} \end{array}} \right], $$

It should be noted that the all the non-zero terms of matrix in Eq. (31) can be described by the modified Lorentzian function. Thus, the expression of current density Jp (p = x, y, z) is described by the product of Ec and the modified Lorentzian function [48].

$${{\boldsymbol J}_p}\textrm{(}\omega \textrm{)} = {\varepsilon _0}\frac{{{a_1}j\omega + {a_0}}}{{{b_2}{{(j\omega )}^2} + {b_1}j\omega + {b_0}}}{{\boldsymbol E}_c}\textrm{(}\omega \textrm{)}, $$
where ${{\boldsymbol J}_p} = {[{{J_{px}},{J_{py}},{J_{pz}}} ]^T}$.

Transforming Eq. (32) from frequency domain to time domain expression, we get:

$${b_2}\frac{{{d^2}{{\boldsymbol J}_p}}}{{d{t^2}}} + {b_1}\frac{{d{{\boldsymbol J}_p}}}{{dt}} + {b_0}{{\boldsymbol J}_p} = {\varepsilon _0}{a_1}\frac{{d{{\boldsymbol E}_c}}}{{dt}} + {\varepsilon _0}{a_0}{{\boldsymbol E}_c}. $$

Using the second order central difference strategy to approximate Eq. (33), we get:

$${b_0}{\boldsymbol J}_p^{n - 1} + {b_1}\frac{{{\boldsymbol J}_p^n - {\boldsymbol J}_p^{n - 2}}}{{2\Delta t}}\textrm{ + }{b_2}\frac{{{\boldsymbol J}_p^n - 2{\boldsymbol J}_p^{n - 1} + {\boldsymbol J}_p^{n - 2}}}{{\Delta {t^2}}} = {a_0}\frac{{{\boldsymbol E}_c^{n - 1/2} + {\boldsymbol E}_c^{n - 3/2}}}{2} + {a_1}\frac{{{\boldsymbol E}_c^{n - 1/2} - {\boldsymbol E}_c^{n - 3/2}}}{{\Delta t}}, $$

Thus, the numerical iterative formula for Jp is obtained as:

$${\boldsymbol J}_p^n = {d_1}{\boldsymbol J}_p^{n - 1} + {d_2}{\boldsymbol J}_p^{n - 2} + {d_3}{\boldsymbol E}_c^{n - 1/2} + {d_4}{\boldsymbol E}_c^{n - 3/2}, $$
where ${d_1} = \frac{{4{b_2} - 2{b_0}\Delta {t^2}}}{{{b_1}\Delta t + 2{b_2}}}$, ${d_2} = \frac{{{b_1}\Delta t - 2{b_2}}}{{{b_1}\Delta t + 2{b_2}}}$, ${d_3} = \frac{{{a_0}\Delta {t^2} + 2{a_1}\Delta t}}{{{b_1}\Delta t + 2{b_2}}}$, ${d_4} = \frac{{{a_0}\Delta {t^2} - 2{a_1}\Delta t}}{{{b_1}\Delta t + 2{b_2}}}$.

According to Eqs. (30)–(35), the updating formulae of current density Jxx, Jxy and Jxz are:

$$J_{xx}^n = {f_1}J_{xx}^{n - 1} + {f_2}J_{xx}^{n - 2} + {f_3}E_{cx}^{n - 1/2} + {f_4}E_{cx}^{n - 3/2}, $$
$$J_{xy}^n = {f_1}J_{xy}^{n - 1} + {f_2}J_{xy}^{n - 2} - {f_5}E_{cy}^{n - 1/2} - {f_6}E_{cy}^{n - 3/2}$$
$$J_{xz}^n = 0, $$
where ${f_1} = (4 - 2\Delta {t^2}(\omega _c^2 + {v^2}))/f$, ${f_2} = (2{v_c}\Delta t - 2)/f$, ${f_3} = ({\varepsilon _0}{v_c}\omega _p^2\Delta {t^2} + 2{\varepsilon _0}\omega _p^2\Delta t)/f$, ${f_4} = ({\varepsilon _0}{v_c}\omega _p^2\Delta {t^2} - 2{\varepsilon _0}\omega _p^2\Delta t)/f$, ${f_5} = {f_6} = ({\varepsilon _0}\omega _p^2{\omega _b}\Delta {t^2})/f$, $f = 2{v_c}\Delta t + 2$.

Substituting Eqs. (36a)–(36c) into $J_{cx}^n = J_{xx}^n + J_{xy}^n + J_{xz}^n$, we get

$$\begin{aligned} J_{cx}^n &= {f_1}(J_{xx}^{n - 1} + J_{xy}^{n - 1} + J_{xz}^{n - 1}) + {f_2}(J_{xx}^{n - 2} + J_{xy}^{n - 2} + J_{xz}^{n - 2})\\& + {f_3}E_{cx}^{n - 1/2} + {f_4}E_{cx}^{n - 3/2} - {f_5}E_{cy}^{n - 1/2} - {f_6}E_{cy}^{n - 3/2}\\& = {f_1}J_{cx}^{n - 1} + {f_2}J_{cx}^{n - 2} + {f_3}E_{cx}^{n - 1/2} + {f_4}E_{cx}^{n - 3/2} - {f_5}E_{cy}^{n - 1/2} - {f_6}E_{cy}^{n - 3/2} \end{aligned}.$$

According to the derived updating formulae, the configurations of these related field components are shown in Fig. 1. The fields Ec (Ecx, Ecy, Ecz) and Hc (Hcx, Hcy, Hcz) are spatially staggered in the Yee cell, and the current density vectors (Jxx, Jxy, Jxz) are collocated with the electric field component Ecx. While solving the current density vectors Jcx (Eq. (37)), the electric vectors Ecx and Ecy are required. Since Ecy is not located at the same position as Jcx. Spatial averaging is needed to obtain the value of Ecy at the spatial position of Jcx. This spatial averaging strategy is implemented by averaging the four nearest neighbors to the location of interest (Eq. (38)), as shown in Fig. 2 [49–51].

$$\begin{aligned} {E_{cy}}(i + 1/2,j,k) &= ({E_{cy}}(i,j + 1/2,k) + {E_{cy}}(i,j - 1/2,k)\\& + {E_{cy}}(i + 1,j + 1/2,k) + {E_{cy}}(i + 1,j - 1/2,k))/4 \end{aligned}.$$

 figure: Fig. 1.

Fig. 1. The configurations of related field components in Yee cell.

Download Full Size | PDF

 figure: Fig. 2.

Fig. 2. The illustration of the spatial averaging scheme for Ecy.

Download Full Size | PDF

Similarly, the numerical iteration of Jcy and Jcz at n time step is given by:

$$J_{cy}^n = {f_1}J_{cy}^{n - 1} + {f_2}J_{cy}^{n - 2} + {f_3}E_{cy}^{n - 1/2} + {f_4}E_{cy}^{n - 3/2} + {f_5}E_{cx}^{n - 1/2} + {f_6}E_{cx}^{n - 3/2}, $$
$$J_{cz}^n = {l_1}J_{cz}^{n - 1} + {l_2}J_{cz}^{n - 2} + {l_3}E_{cz}^{n - 1/2} - {l_4}E_{cz}^{n - 3/2}, $$
where ${l_1} = {4 / f}$, ${l_2} = {{({v_c} - 2)} / f}$, ${l_3} = {l_4} = {{2{\varepsilon _0}\omega _p^2} / f}$.

Following the same average strategy as Eq. (38), the value of Ecx at position (i, j + 1/2, k) is approximated by:

$$\begin{aligned} {E_{cx}}(i,j + 1/2,k) &= ({E_{cx}}(i + 1/2,j,k) + {E_{cx}}(i - 1/2,j,k)\\& + {E_{cx}}(i + 1/2,j + 1,k) + {E_{cx}}(i - 1/2,j + 1,k))/4 \end{aligned}.$$
  • B. Magnetic graphene's conductivity

Assuming that a 2-D graphene sheet is placed on the x-y plane of a 3-D Yee grid and biased with a static magnetic field B0 in the z direction, the surface conductivity behavior of the graphene sheet in the microwave frequency range of 10 THz can be described as:

$${\widetilde {\sigma }_s}(\omega ) = \left[ {\begin{array}{{ccc}} {{\sigma_{xx}}(\omega )}&{{\sigma_{xy}}(\omega )}&0\\ {{\sigma_{yx}}(\omega )}&{{\sigma_{yy}}(\omega )}&0\\ 0&0&0 \end{array}} \right] = \left[ {\begin{array}{{ccc}} {{\sigma_d}(\omega )}&{ - {\sigma_g}(\omega )}&0\\ {{\sigma_g}(\omega )}&{{\sigma_d}(\omega )}&0\\ 0&0&0 \end{array}} \right], $$
and
$${\sigma _d}(\omega ) = \frac{{j\omega {\sigma _0} + v{\sigma _0}}}{{(\omega _c^2 + {v^2}) + 2j\omega v + {{(j\omega )}^2}}}, $$
$${\sigma _g}(\omega ) = \frac{{ - {\sigma _0}{\omega _c}}}{{(\omega _c^2 + {v^2}) + 2j\omega v + {{(j\omega )}^2}}}, $$
$${\sigma _0} = \frac{{2{e^2}{k_B}T}}{{\pi {\hbar ^2}}}\ln (2\cosh \frac{{{\mu _c}}}{{2{k_B}T}}), $$
$${\omega _c} = e{B_0}v_F^2/{\mu _c}, $$
where v is the phenomenological scattering rate, ωc is the cyclotron frequency, e is electronic charge, vF is Fermi velocity, kB is Boltzmann constant, T is operating temperature, ħ is reduced Planck constant, B0 is the biased static magnetic field, and µc is the chemical potential. The graphene volume conductivity is defined as $\widetilde \sigma = {{{{\widetilde \sigma }_s}} / {\Delta z}}$. Substituting susceptibility function $\widetilde \chi (\omega ) = \widetilde \sigma /j\omega {\varepsilon _0}$ into ${{\boldsymbol J}_c}\textrm{(}\omega \textrm{)} = j\omega {\varepsilon _0}\widetilde \chi \textrm{(}\omega \textrm{)}{{\boldsymbol E}_c}\textrm{(}\omega \textrm{)}$, the equivalent current density of graphene volume conductivity is obtained:
$$\begin{aligned}{{\boldsymbol J}_c}(\omega ) &= \frac{1}{{\Delta z}}\left[ {\begin{array}{{ccc}} {\frac{{{\sigma_0}j\omega + v{\sigma_0}}}{{(\omega_c^2 + {v^2}) + 2j\omega v + {{(j\omega )}^2}}}}&{\frac{{ - {\omega_c}{\sigma_0}}}{{(\omega_c^2 + {v^2}) + 2j\omega v + {{(j\omega )}^2}}}}&0\\ {\frac{{{\omega_c}{\sigma_0}}}{{(\omega_c^2 + {v^2}) + 2j\omega v + {{(j\omega )}^2}}}}&{\frac{{{\sigma_0}j\omega + v{\sigma_0}}}{{(\omega_c^2 + {v^2}) + 2j\omega v + {{(j\omega )}^2}}}}&0\\ 0&0&0 \end{array}} \right]\left[ {\begin{array}{{c}} {{E_{cx}}}\\ {{E_{cy}}}\\ {{E_{cz}}} \end{array}} \right]\\& = \left[ {\begin{array}{{ccc}} {{J_{xx}}(\omega )}&{{J_{xy}}(\omega )}&0\\ {{J_{yx}}(\omega )}&{{J_{yy}}(\omega )}&0\\ 0&0&0 \end{array}} \right] \end{aligned}.$$

Similarly, the updating formulae of Jcx and Jcy at n time step are given by:

$$J_{cx}^n = {f_1}J_{cx}^{n - 1} + {f_2}J_{cx}^{n - 2} + {f_3}E_{cx}^{n - 1/2} + {f_4}E_{cx}^{n - 3/2} - {f_5}E_{cy}^{n - 1/2} - {f_6}E_{cy}^{n - 3/2}, $$
$$J_{cy}^n = {f_1}J_{cy}^{n - 1} + {f_2}J_{cy}^{n - 2} + {f_3}E_{cy}^{n - 1/2} + {f_4}E_{cy}^{n - 3/2} + {f_5}E_{cx}^{n - 1/2} + {f_6}E_{cx}^{n - 3/2}, $$
where $f = 2v\Delta t + 2$, ${f_1} = (4 - 2\Delta {t^2}(\omega _c^2 + {v^2}))/f$, ${f_2} = (2v\Delta t - 2)/f$, ${f_3} = (v{\sigma _0}\Delta {t^2} + 2{\sigma _0}\Delta t)/(f\Delta z)$, ${f_4} = (v{\sigma _0}\Delta {t^2} - 2{\sigma _0}\Delta t)/(f\Delta z)$, ${f_5} = {f_6} = ({\sigma _0}{\omega _c}\Delta {t^2})/(f\Delta z)$.

Additionally, the treatment of Ecy and Ecx in Eqs. (48) and (49) are constant with the above approach.

2.2 Stability analysis

For the purpose of ascertaining the numerical stability of the anisotropic frequency-dependent ADE-CDI-FDTD method, the numerical iterative expressions for the electric and magnetic field components of graphene with a static magnetic field bias oriented in the z-direction, without any loss of generality, are presented in matrix notation as follows:

$${{\boldsymbol M}_L}\left[ {\begin{array}{{c}} {{\boldsymbol E}_c^{n + 1/2}}\\ {{\boldsymbol H}_c^{n + 1}}\\ {{\boldsymbol J}_c^n} \end{array}} \right] = {{\boldsymbol M}_{R1}}\left[ {\begin{array}{{c}} {{\boldsymbol E}_c^{n - 1/2}}\\ {{\boldsymbol H}_c^n}\\ {{\boldsymbol J}_c^{n - 1}} \end{array}} \right] + {{\boldsymbol M}_{R2}}\left[ {\begin{array}{{c}} {{\boldsymbol E}_c^{n - 3/2}}\\ {{\boldsymbol H}_c^{n - 1}}\\ {{\boldsymbol J}_c^{n - 2}} \end{array}} \right], $$
and
$${{\boldsymbol M}_L} = \left[ {\begin{array}{{ccccccccc}} 1&0&0&0&0&0&a&0&0\\ 0&1&0&0&0&0&0&a&0\\ 0&0&1&0&0&0&0&0&0\\ 0&{ - b{D_z}{F_z}}&{b{D_y}{F_x}}&1&0&0&0&0&0\\ {b{D_z}{F_y}}&0&{ - b{D_x}{F_x}}&0&1&0&0&0&0\\ { - b{D_y}{F_y}}&{b{D_x}{F_z}}&0&0&0&1&0&0&0\\ 0&0&0&0&0&0&1&0&0\\ 0&0&0&0&0&0&0&1&0\\ 0&0&0&0&0&0&0&0&1 \end{array}} \right], $$
$${{\boldsymbol M}_{R1}} = \left[ {\begin{array}{{ccccccccc}} 1&0&0&0&{ - a{D_z}{F_z}}&{a{D_y}{F_x}}&0&0&0\\ 0&1&0&{a{D_z}{F_y}}&0&{ - a{D_x}{F_x}}&0&0&0\\ 0&0&1&{ - a{D_y}{F_y}}&{a{D_x}{F_z}}&0&0&0&0\\ 0&0&0&1&0&0&0&0&0\\ 0&0&0&0&1&0&0&0&0\\ 0&0&0&0&0&1&0&0&0\\ {{f_3}}&{ - {f_5}}&0&0&0&0&{{f_1}}&0&0\\ {{f_5}}&{{f_3}}&0&0&0&0&0&{{f_1}}&0\\ 0&0&0&0&0&0&0&0&0 \end{array}} \right], $$
$${{\boldsymbol M}_{R2}} = \left[ {\begin{array}{{ccccccccc}} 0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0\\ {{f_4}}&{ - {f_6}}&0&0&0&0&{{f_2}}&0&0\\ {{f_6}}&{{f_4}}&0&0&0&0&0&{{f_2}}&0\\ 0&0&0&0&0&0&0&0&0 \end{array}} \right], $$
where a = Δt/ε0, b = Δt/μ0, Dp = ∂/∂p (p = x, y, z) represents a first-order derivative operator with respect to p, ${F_p} = \textrm{ }1 - \frac{{\Delta {t^2}D_p^2}}{{4{\varepsilon _0}{\mu _0}}}$. And the field components of electromagnetic waves can be represented as:
$$\phi _{i,j,k}^n = {\varphi _\phi }{\zeta ^n}\exp [\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over j} ({k_x}i\Delta x + {k_y}j\Delta y + {k_z}k\Delta z)], $$
where $\phi \textrm{ = }E,H,J$, φϕ denotes the amplitude of the field component. ζ denotes the time growth factor, $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over j} = \sqrt { - 1}$, kp (p = x, y, z) denotes the wavenumber in the p-direction. The subscript i, j, and k represent the index of the Yee grid, and Δp (p = x, y, z) represents the grid size.

Then, using a central second-order finite difference to approximate the spatial derivative, we get:

$${D_p}\phi _{i,j,k}^n = \frac{{\phi _{i,j,k}^n(p + \frac{{\Delta p}}{2}) - \phi _{i,j,k}^n(p - \frac{{\Delta p}}{2})}}{{\Delta p}} = {\sigma _p}\phi _{i,j,k}^n, $$
where, ${\sigma _p} = 2\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over j} \sin ({k_p}\Delta p/2)/\Delta p$, (p = x, y, z).

Finally, substituting Eq. (53) into Eq. (48), we get:

$$({{\boldsymbol M}_{R1}} + {{{{\boldsymbol M}_{R2}}} / \zeta } - \zeta {{\boldsymbol M}_L})\left[ {\begin{array}{{c}} {{\boldsymbol E}_c^{n - 1/2}}\\ {{\boldsymbol H}_c^n}\\ {{\boldsymbol J}_c^{n - 1}} \end{array}} \right] = {{\boldsymbol M}}\left[ {\begin{array}{{c}} {{\boldsymbol E}_c^{n - 1/2}}\\ {{\boldsymbol H}_c^n}\\ {{\boldsymbol J}_c^{n - 1}} \end{array}} \right] = 0,$$
where the matrix M is:
$$\scalebox{0.8}{$\displaystyle{\boldsymbol M} = \left[ {\begin{array}{{ccccccccc}} {1 - \zeta }&0&0&0&{ - a{\sigma_z}{F_z}}&{a{\sigma_y}{F_x}}&{ - a\zeta }&0&0\\ 0&{1 - \zeta }&0&{a{\sigma_z}{F_y}}&0&{ - a{\sigma_x}{F_x}}&0&{ - a\zeta }&0\\ 0&0&{1 - \zeta }&{ - a{\sigma_y}{F_y}}&{a{\sigma_x}{F_z}}&0&0&0&0\\ 0&{b\zeta {\sigma_z}F}&{ - b\zeta {\sigma_y}{F_x}}&{1 - \zeta }&0&0&0&0&0\\ { - b\zeta {\sigma_z}{F_y}}&0&{b\zeta {\sigma_x}{F_x}}&0&{1 - \zeta }&0&0&0&0\\ {b\zeta {\sigma_y}{F_y}}&{ - b\zeta {\sigma_x}{F_z}}&0&0&0&{1 - \zeta }&0&0&0\\ {{f_3} + {{{f_4}} / \zeta }}&{ - {f_5} - {{{f_6}} / \zeta }}&0&0&0&0&{{f_1} + {{{f_2}} / {\zeta - \zeta }}}&0&0\\ {{f_5} + {{{f_6}} / \zeta }}&{{f_3} + {{{f_4}} / \zeta }}&0&0&0&0&0&{{f_1} + {{{f_2}} / {\zeta - \zeta }}}&0\\ 0&0&0&0&0&0&0&0&{ - \zeta } \end{array}} \right]$}$$

To guarantee the existence of a non-zero solution to Eq. (56), the determinant of the coefficient matrix M must be 0. Moreover, to ensure numerical stability during the iterative process, the modulus of the growth factor ζ must be less than or equal to 1. In addition to the growth factor ζ, according to Eqs. (25)–(28), (48) and (49), the parameters in the coefficient matrix M include Δt, Δx, Δy, Δz, a, b, ${f_1}$, ${f_2}$, ${f_3}$, ${f_4}$, ${f_5}$, and ${f_6}$. Therefore, a simple and feasible method is to substitute the parameters (Δx = Δy = Δz = 1 nm, $\Delta {t_{FDTD}} = \frac{1}{{{c_0}\sqrt {1/\Delta {x^2} + 1/\Delta {y^2} + 1/\Delta {z^2}} }}$, B0 = 1 T, v = 2.148 × 1011 Hz, T = 300 K, μc = 0.1 eV, vF = 0.96 × 106 m/s) into matrix M to determine the modulus of the amplification factor ζ. The range of values of ${k_p}\Delta p/2$ in ${\sigma _p}$(${\sigma _p} = 2\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over j} \sin ({k_p}\Delta p/2)/\Delta p$) is set as [0, π]. As shown in Fig. 3, the resulting modulus value of the growth factor ζ does not exceed 1 with the ADE-CDI-FDTD method when the CFLN = 5, 10 and 15 (CFLN represents $\Delta {t_{\textrm{CDI}}}/\Delta {t_{\textrm{FDTD}}}$). Consequently, the unconditionally stable property of the frequency-dependent ADE-CDI-FDTD method is verified.

 figure: Fig. 3.

Fig. 3. The modulus value of ζ is obtained for values of ${{{k_p}\Delta p/2}}$ in the ADE-CDI-FDTD method with CFLN (ΔtCDI/ΔtFDTD) = 5, 10, 15.

Download Full Size | PDF

3. Numerical examples and analysis

3.1 Magnetized graphene sheets

First, to ensure the computational soundness of the proposed approach, the current investigation focuses on examining the transmission characteristics of graphene sheets under static magnetic field B0 bias in terahertz and optical bands. Fig. 4 illustrates that the entire computational domain covers 20 × 20 × 180 cells in the x, y and z directions, with a grid cell size of Δx = Δy = Δz = 1 µm. Specifically, a 10-level CPML is positioned at both ends of the z-direction to truncate the computational space, while periodic boundaries are established in the x and y directions. The plane wave propagates in the z-direction in the form of a Gaussian pulse (This expression is $\textrm{pulse(}t\textrm{)} = \exp ( - 4\pi {(t - {t_0})^2}/{\tau ^2})$, where τ = 100Δt/CFLN, and τ = 2t0.), with the graphene sheet placed at the x-y plane. The relevant parameters for the graphene sheet are: B0 = 1 T, v = 2.148 × 1011 Hz, T = 300 K, µc = 0.1 eV, vF = 0.96 × 106 m/s. As depicted in Fig. 5, the transmission coefficient curves of graphene are determined using the analytical method [44], the ADE-FDTD method, the ADE-CDI-FDTD method (CFLN = 5, 10, 15), and the ADE-LADI-FDTD method (CFLN = 5) [25]. The obtained results indicate that the proposed method yields consistent results with those obtained through conventional FDTD and analytical methods, and exhibits superior numerical accuracy compared to the ADE-LADI-FDTD method (CFLN = 5). Concurrently, to quantify the difference between the calculation results of the numerical methods and the analytical solution, the formula of defining relative calculation error is expressed by Eq. (57):

$${E_r}(f) = \frac{{|{{T_{num}}(f) - {T_{Analytical}}(f)} |}}{{\max |{{T_{Analytical}}(f)} |}} \times 100\%$$
where ${T_{num}}$ represents the calculation of transmission coefficients obtained by the ADE-LADI-FDTD method (CFLN = 5) and the ADE-CDI-FDTD method (CFLN = 5, 10, 15), respectively, while ${T_{Analytical}}$ represents the transmission coefficients calculated by the analytical method.

 figure: Fig. 4.

Fig. 4. Model diagram of magnetized graphene sheet.

Download Full Size | PDF

 figure: Fig. 5.

Fig. 5. The transmission coefficient of the magnetized graphene sheet versus frequency.

Download Full Size | PDF

As illustrated in Fig. 6, it is evident that the relative calculation error of the classical LADI-FDTD method (CFLN = 5) is large than that of the ADE-CDI-FDTD method (CFLN = 5). Fig. 7 shows the snapshots of the electric field components Ey recorded at the yoz plane for ADE-FDTD method (CFLN = 1) and ADE-CDI-FDTD (CFLN = 5) method, depicting consistent field distributions, thereby validating the accuracy and reliability of the proposed approach. Simultaneously, to test the late-time stability of the proposed method, a probe field Ex is recorded for 180 ps (the time step of the ADE-CDI-FDTD method is chosen as Δt = 9.6225 fs, corresponding to the condition of CFLN = 5). As shown in Fig. 8 the amplitude of the recorded Ex decay to zero and remain stable, which demonstrate the unconditionally stable property of the ADE-CDI-FDTD method.

 figure: Fig. 6.

Fig. 6. Relative calculation error of the transmission coefficient for ADE-LADI-FDTD and ADE-FDTD methods.

Download Full Size | PDF

 figure: Fig. 7.

Fig. 7. Snapshots of the Ey component in the y-z plane for ADE-FDTD (CFLN = 1) and ADE-CDI-FDTD (CFLN = 5) methods.

Download Full Size | PDF

 figure: Fig. 8.

Fig. 8. Recorded Ex versus time calculated by ADE-CDI-FDTD (CFLN = 5) method.

Download Full Size | PDF

3.2 Magnetized plasma slab

Subsequently, to assess the validation of the ADE-CDI-FDTD algorithm when managing anisotropic dispersive medium using the J-E constitutive relation, a magnetized plasma slab model with a thickness of 9 mm is employed, as depicted in Fig. 9. The computational domain encompasses 20 × 20 × 180 cells in the x, y, and z directions, with a spatial grid discretization of Δx = Δy = Δz = 1 µm and a time step of ΔtFDTD = 0.125 ps. The magnetized plasma slab (with parameters of ωp = (2π) 50 × 109 rad/s, ωb = 3.0 × 1011 rad/s, vc = 2.0 × 1010 [52]) occupies 120 grid cells in the z direction. A Gaussian pulse plane wave is emitted in free space along the z direction. Upon reaching the interface, the magnetized plasma undergoes internal evolution as a result of its anisotropy. The reflection and transmission planes are situated on either side of the plasma slab, and the mean values of the electric field (Er and Et) are subsequently computed.

 figure: Fig. 9.

Fig. 9. Model of magnetized plasma slab with width of 9mm.

Download Full Size | PDF

The reflection and transmission coefficients of both right circularly polarized (RCP) and left circularly polarized (LCP) waves are depicted in Fig. 10 and Fig. 11, and are computed using the analytical method [53], as well as the ADE-FDTD method and the proposed ADE-CDI-FDTD method (CFLN = 1, 3, 5). It can be observed that the results produced by the ADE-CDI-FDTD method are in excellent agreement with both the analytical solution and the results obtained from the conventional FDTD method. These findings further demonstrate the efficacy and accuracy of the proposed methodology.

 figure: Fig. 10.

Fig. 10. (a) RCP reflection coefficient, (b) LCP reflection coefficient for plasma with a biasing magnetic field parallel to the propagation direction.

Download Full Size | PDF

 figure: Fig. 11.

Fig. 11. (a) RCP transmission coefficient, and (b) LCP transmission coefficient for plasma with a biasing magnetic field parallel to the propagation direction.

Download Full Size | PDF

3.3 Simulation of electromagnetic wave propagation in a plasma

Finally, to simulate the propagation characteristics of electromagnetic waves in a lossless plasma medium under the influence of a static magnetic field in the z-direction, the ADE-CDI-FDTD method is utilized and compared with the conventional FDTD method. Fig. 12 depicts the simulation results obtained using a computational domain spanning 40 × 40 × 100 grid cells in the x, y, and z directions, with a grid size of Δx = Δy = Δz = 1 mm. To truncate the computational domain, 10 layers of CPML are incorporated. The plasma medium is cubic in shape and has a side length of 10 mm. Its parameters are given by ωp = 2π × 10 Grad/s, ωb = 10 Grad/s, vc = 10 GHz [25]. An excitation source in the form of a modulated Gaussian pulse is positioned at point A, with its temporal expression provided as:

$$\textrm{pulse}(t) = \cos (2\pi {f_0}t)\exp [ - 4\pi {(t - {t_0})^2}/{\tau ^2}], $$
where f0 = 10 GHz, τ = 240Δt / CFLN, t0 = 6τ.

 figure: Fig. 12.

Fig. 12. Simulation of 3D plasma block model.

Download Full Size | PDF

The validity of the proposed ADE-CDI-FDTD method is demonstrated in Fig. 13, which presents the electric field component Ez recorded at detection point B by both the ADE-FDTD method and the ADE-CDI-FDTD method (CFLN = 1, 3, 5). Notably, when CFLN is set to 5, the time-domain waveform calculated by the ADE-CDI-FDTD method matches that calculated by the conventional FDTD method, providing evidence for the soundness of the proposed method in this study. Moreover, computational efficiency comparisons were conducted between the ADE-CDI-FDTD method and the conventional FDTD method, with CPU times listed in Table 1. Results indicate that when CFLN equals 3, the ADE-CDI-FDTD method exhibits significantly shorter execution times than the traditional FDTD method, underscoring the efficiency of the proposed ADE-CDI-FDTD method in solving problems related to anisotropic dispersive medium.

 figure: Fig. 13.

Fig. 13. The waveform at the probe point B versus of time.

Download Full Size | PDF

Tables Icon

Table 1. THE CPU TIME COMPARISION BASED ON TWO METHODS

View Table

Funding

National Natural Science Foundation of China (2022YFA1404003, 61901001, 62101002, 62201003, U20A20164).

Disclosures

The authors declare no conflicts of interest.

Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

References

1. A. Hartman, D. Romano, G. Antonini, and J. Ekman, “Partial element equivalent circuit models of three-dimensional geometries including anisotropic dielectrics,” IEEE Trans. Electromagn. Compat. 60(3), 696–704 (2018). [CrossRef]  

2. L. Guan, S. Tao, Y. Zhao, and R. Chen, “An anisotropic thin dielectric sheet method for analysis of magnetized graphene,” Antennas Wirel. Propag. Lett. 18(2), 279–282 (2019). [CrossRef]  

3. X. Li, G. Wei, S. Lei, K. Han, T. Qiu, G. Zhang, and Y. Zhou, “Bifunctional Luneburg-Eaton lens fabricated of 3-D-printed anisotropic medium,” Antennas Wirel. Propag. Lett. 21(7), 1462–1466 (2022). [CrossRef]  

4. V. Kamra and A. Dreher, “Efficient analysis of multiple microstrip transmission lines with anisotropic substrates,” IEEE Microw. Wireless Compon. Lett. 28(8), 636–638 (2018). [CrossRef]  

5. J. Wang, J. Li, Y. Chen, F. Han, and Q. H. Liu, “Simulation of 3-D electromagnetic scattering and inverse scattering by arbitrary anisotropic dielectric objects embedded in layered arbitrary anisotropic media,” IEEE Trans. Antennas Propagat. 68(8), 6473–6478 (2020). [CrossRef]  

6. A. Moradi, “Reflection of electrostatic waves on plane interface between magnetized plasma and insulator,” IEEE Trans. Plasma Sci. 48(7), 2687–2690 (2020). [CrossRef]  

7. Y. Yu and J. J. Simpson, “An EJ collocated 3-D FDTD model of electromagnetic wave propagation in magnetized cold plasma,” IEEE Trans. Antennas Propag. 58(2), 469–478 (2009). [CrossRef]  

8. E. Gamliel, “Direct integration 3-D FDTD method for single-species cold magnetized plasma,” IEEE Trans. Antennas Propagat. 65(1), 295–308 (2017). [CrossRef]  

9. S. Liu, S. Liu, and S. Liu, “Analysis for scattering of conductive objects covered with anisotropic magnetized plasma by trapezoidal recursive convolution finite-difference time-domain method,” Int J RF and Microwave Comp Aid Eng 20(4), 465–472 (2010). [CrossRef]  

10. Y. Zhang, Y. Liu, and X. Li, “Total-field/scattered-field formulation for FDTD analysis of plane-wave propagation through cold magnetized plasma sheath,” IEEE Trans. Antennas Propagat. 68(1), 377–387 (2020). [CrossRef]  

11. X. Jia, F. Yang, M. Li, and S. Xu, “Fast Analysis of Metasurfaces Composed of Asymmetric Scatterers Using FDTD With Anisotropic Surface Susceptibility Model,” Antennas Wirel. Propag. Lett. 22(1), 119–123 (2023). [CrossRef]  

12. M. Salmasi, M. E. Potter, and M. Okoniewski, “Implementation of general dispersive anisotropic materials in Lebedev FDTD,” IEEE Trans. Antennas Propagat. 66(12), 7171–7179 (2018). [CrossRef]  

13. A. Pereda, L. A. Vielva, A. Vegas, and A. Prieto, “Analyzing the stability of the FDTD technique by combining the von Neumann method with the Routh-Hurwitz criterion,” IEEE Trans. Microwave Theory Techn. 49(2), 377–381 (2001). [CrossRef]  

14. J. Park and K.-Y. Jung, “Numerical Stability of Modified Lorentz FDTD Unified From Various Dispersion Models,” Opt. Express 29(14), 21639–21654 (2021). [CrossRef]  

15. K. P. Prokopidis and D. C. Zografopoulos, “Investigation of the stability of ADE-FDTD methods for modified Lorentz media,” IEEE Microw. Wireless Compon. Lett. 24(10), 659–661 (2014). [CrossRef]  

16. H. Choi, J.-W. Baek, and K.-Y. Jung, “Comprehensive study on numerical aspects of modified Lorentz model-based dispersive FDTD formulations,” IEEE Trans. Antennas Propagat. 67(12), 7643–7648 (2019). [CrossRef]  

17. H. Choi, J.-W. Baek, and K.-Y. Jung, “Numerical stability and accuracy of CCPR-FDTD for dispersive media,” IEEE Trans. Antennas Propagat.IEEE Trans. Antennas Propagat. 68(11), 7717–7720 (2020). [CrossRef]  

18. K. Hosseini and Z. Atlasbaf, “Unconditionally stable LOD-FDTD method in anisotropic magnetized plasma,” IEEE Microw. Wireless Compon. Lett. 27(3), 212–214 (2017). [CrossRef]  

19. K. K. Hosseini and Z. Atlasbaf, “Efficient three-step LOD-FDTD method in lossy saturated ferrites with arbitrary magnetization,” IEEE Trans. Antennas Propagat. 66(3), 1374–1383 (2018). [CrossRef]  

20. Y.-F. Mao, B. Chen, J.-L. Xia, J. Chen, and J.-Z. Tang, “Application of the leapfrog ADI FDTD method to periodic structures,” Antennas Wirel. Propag. Lett. 12, 599–602 (2013). [CrossRef]  

21. F. Jolani, Y. Yu, and Z. Chen, “Efficient modeling of open structures using nonuniform leapfrog ADI-FDTD,” Antennas Wirel. Propag. Lett. 10, 561–564 (2011). [CrossRef]  

22. B. Zou, S. Liu, L. Zhang, and S. Ren, “Efficient one-step leapfrog ADI-FDTD for far-field scattering calculation of lossy media,” Microw. Opt. Technol. Lett. 62(5), 1876–1881 (2020). [CrossRef]  

23. S.-C. Yang, Z. D. Chen, Y.-q. Yu, and W.-Y. Yin, “The unconditionally stable one-step leapfrog ADI-FDTD method and its comparisons with other FDTD methods,” IEEE Microw. Wireless Compon. Lett. 21(12), 640–642 (2011). [CrossRef]  

24. S.-C. Yang, Z. Chen, Y. Yu, and W.-Y. Yin, “An unconditionally stable one-step arbitrary-order leapfrog ADI-FDTD method and its numerical properties,” IEEE Trans. Antennas Propagat. 60(4), 1995–2003 (2012). [CrossRef]  

25. X.-H. Wang, W.-Y. Yin, and Z. Z. Chen, “One-step leapfrog ADI-FDTD method for simulating electromagnetic wave propagation in general dispersive media,” Opt. Express 21(18), 20565–20576 (2013). [CrossRef]  

26. K. P. Prokopidis and D. C. Zografopoulos, “One-step leapfrog ADI-FDTD method using the complex-conjugate pole-residue pairs dispersion model,” IEEE Microw. Wireless Compon. Lett. 28(12), 1068–1070 (2018). [CrossRef]  

27. X.-H. Wang, J.-Y. Gao, Z. Chen, and F. L. Teixeira, “Unconditionally stable one-step leapfrog ADI-FDTD for dispersive media,” IEEE Trans. Antennas Propagat. 67(4), 2829–2834 (2019). [CrossRef]  

28. G. Xie, Z. Huang, M. Fang, and X. Wu, “A unified 3-D ADI-FDTD algorithm with one-step leapfrog approach for modeling frequency-dependent dispersive media,” Int J Numer Model 33(2), e2666 (2020). [CrossRef]  

29. T. H. Gan and E. L. Tan, “Analysis of the divesrgence properties for the three-dimensional leapfrog ADI-FDTD method,” IEEE Trans. Antennas Propagat. 60(12), 5801–5808 (2012). [CrossRef]  

30. E. L. Tan, “From time-collocated to leapfrog fundamental schemes for ADI and CDI FDTD methods,” Axioms 11(1), 23 (2022). [CrossRef]  

31. E. L. Tan, “A leapfrog scheme for complying-divergence implicit finite-difference time-domain method,” Antennas Wirel. Propag. Lett. 20(5), 853–857 (2021). [CrossRef]  

32. H. Gan and E. L. Tan, “Unconditionally stable fundamental LOD-FDTD method with second-order temporal accuracy and complying divergence,” IEEE Trans. Antennas Propagat. 61(5), 2630–2638 (2013). [CrossRef]  

33. S. Liu, E. L. Tan, and B. Zou, “Incident plane-wave source formulations for leapfrog complying-divergence implicit FDTD method,” IEEE J. Multiscale Multiphys. Comput. Tech. 7, 84–91 (2022). [CrossRef]  

34. S. Liu, E. L. Tan, L. Zhang, and B. Zou, “Unified Formulation of Leapfrog Complying-Divergence Implicit FDTD Method for Lossy Media,” IEEE Trans. Antennas Propag. (to be published).

35. S. Liu, E. L. Tan, B. Zou, and L. Zhang, “Higher Order CPML for Leapfrog Complying-Divergence Implicit FDTD Method and Its Numerical Properties,” IEEE Trans. Microw. Theory Tech. 71(2), 522–535 (2023). [CrossRef]  

36. D. M. Sullivan, “Frequency-dependent FDTD methods using Z transforms,” IEEE Trans. Antennas Propagat. 40(10), 1223–1230 (1992). [CrossRef]  

37. Y. Zhang, N. Feng, J. Zhu, G. Xie, L. Yang, and Z. Huang, “Z-Transform-Based FDTD Implementations of Biaxial Anisotropy for Radar Target Scattering Problems,” Remote Sensing 14(10), 2397 (2022). [CrossRef]  

38. G. Xie, M. Fang, Z. Huang, X. Ren, and X. Wu, “A unified 3-D simulating framework for Debye-type dispersive media and PML technique based on recursive integral method,” Comput. Phys. Commun. 280, 108463 (2022). [CrossRef]  

39. G. Xie, M. Fang, Z. Huang, X. Wu, X. Ren, and N. Feng, “A Numerical Study of Lossy Multipole Debye Dispersive Media Using a Recursive Integral-FDTD Method,” IEEE Trans. Microwave Theory Techn. 71(3), 1009–1018 (2023). [CrossRef]  

40. G. Xie, Z. Song, G. Hou, M. Fang, N. Feng, and Z. Huang, “Efficacious GPR Implementations of Z-Transform-Based Hybrid LOD-FDTD with Subgridding Scheme: Theoretical Formalism and Numerical Study,” Remote Sensing 14(21), 5393 (2022). [CrossRef]  

41. G.-Y. Liu, W.-J. Chen, and J. Quan, “A general ADE-WLP-FDTD method with a new temporal basis for wave propagation,” IEEE Microw. Wireless Compon. Lett. 32(12), 1391–1394 (2022). [CrossRef]  

42. A. A. Al-Jabr, M. A. Alsunaidi, T. Khee, and B. S. Ooi, “A Simple FDTD Algorithm for Simulating EM-Wave Propagation in General Dispersive Anisotropic Material,” IEEE Trans. Antennas Propagat. 61(3), 1321–1326 (2013). [CrossRef]  

43. Y. Zhou, H. Li, L. Li, Y. Cai, K. Zeyde, and X. Han, “Efficient HIE-FDTD method for designing a dual-band anisotropic terahertz absorption structure,” Opt. Express 29(12), 18611–18623 (2021). [CrossRef]  

44. X.-H. Wang, W.-Y. Yin, and Z. Chen, “Matrix exponential FDTD modeling of magnetized graphene sheet,” Antennas Wirel. Propag. Lett. 12, 1129–1132 (2013). [CrossRef]  

45. X.-H. Wang, J.-Y. Gao, and F. L. Teixeira, “Stability-improved ADE-FDTD method for wideband modeling of graphene structures,” Antennas Wirel. Propag. Lett. 18(1), 212–216 (2019). [CrossRef]  

46. M. Feizi, V. Nayyeri, and O. M. Ramahi, “Modeling magnetized graphene in the finite-difference time-domain method using an anisotropic surface boundary condition,” IEEE Trans. Antennas Propagat. 66(1), 233–241 (2018). [CrossRef]  

47. O. Ramadan, “Simplified FDTD Formulations for Magnetostatic Biased Graphene Simulations,” Antennas Wirel. Propag. Lett. 19(12), 2290–2294 (2020). [CrossRef]  

48. K. P. Prokopidis and D. C. Zografopoulos, “An ADI-FDTD formulation with modified Lorentz dispersion for the study of plasmonic structures,” IEEE Photon. Technol. Lett. 26(22), 2267–2270 (2014). [CrossRef]  

49. S. Pokhrel, V. Shankar, and J. J. Simpson, “3-D FDTD modeling of electromagnetic wave propagation in magnetized plasma requiring singular updates to the current density equation,” IEEE Trans. Antennas Propagat. 66(9), 4772–4781 (2018). [CrossRef]  

50. Y. Yu, J. Niu, and J. J. Simpson, “A 3-D global earth-ionosphere FDTD model including an anisotropic magnetized plasma ionosphere,” IEEE Trans. Antennas Propagat. 60(7), 3246–3256 (2012). [CrossRef]  

51. A. Samimi and J. J. Simpson, “An efficient 3-D FDTD model of electromagnetic wave propagation in magnetized plasma,” IEEE Trans. Antennas Propagat. 63(1), 269–279 (2015). [CrossRef]  

52. X. Xi, Z. Li, J. Liu, and J. Zhang, “FDTD simulation for wave propagation in anisotropic dispersive material based on bilinear transform,” IEEE Trans. Antennas Propagat. 63(11), 5134–5138 (2015). [CrossRef]  

53. F. Hunsberger, R. Luebbers, and K. Kunz, “Finite-difference time-domain analysis of gyrotropic media. I. Magnetized plasma,” IEEE Trans. Antennas Propagat. 40(12), 1489–1495 (1992). [CrossRef]  

Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (13)
Fig. 1.
Fig. 1. The configurations of related field components in Yee cell.

View in Article | Download Full Size | PDF

Fig. 2.
Fig. 2. The illustration of the spatial averaging scheme for Ecy.

View in Article | Download Full Size | PDF

Fig. 3.
Fig. 3. The modulus value of ζ is obtained for values of ${{{k_p}\Delta p/2}}$ in the ADE-CDI-FDTD method with CFLN (ΔtCDI/ΔtFDTD) = 5, 10, 15.

View in Article | Download Full Size | PDF

Fig. 4.
Fig. 4. Model diagram of magnetized graphene sheet.

View in Article | Download Full Size | PDF

Fig. 5.
Fig. 5. The transmission coefficient of the magnetized graphene sheet versus frequency.

View in Article | Download Full Size | PDF

Fig. 6.
Fig. 6. Relative calculation error of the transmission coefficient for ADE-LADI-FDTD and ADE-FDTD methods.

View in Article | Download Full Size | PDF

Fig. 7.
Fig. 7. Snapshots of the Ey component in the y-z plane for ADE-FDTD (CFLN = 1) and ADE-CDI-FDTD (CFLN = 5) methods.

View in Article | Download Full Size | PDF

Fig. 8.
Fig. 8. Recorded Ex versus time calculated by ADE-CDI-FDTD (CFLN = 5) method.

View in Article | Download Full Size | PDF

Fig. 9.
Fig. 9. Model of magnetized plasma slab with width of 9mm.

View in Article | Download Full Size | PDF

Fig. 10.
Fig. 10. (a) RCP reflection coefficient, (b) LCP reflection coefficient for plasma with a biasing magnetic field parallel to the propagation direction.

View in Article | Download Full Size | PDF

Fig. 11.
Fig. 11. (a) RCP transmission coefficient, and (b) LCP transmission coefficient for plasma with a biasing magnetic field parallel to the propagation direction.

View in Article | Download Full Size | PDF

Fig. 12.
Fig. 12. Simulation of 3D plasma block model.

View in Article | Download Full Size | PDF

Fig. 13.
Fig. 13. The waveform at the probe point B versus of time.

View in Article | Download Full Size | PDF

Tables (1)

Tables Icon

Table 1. THE CPU TIME COMPARISION BASED ON TWO METHODS

View Table

Equations (63)

Equations on this page are rendered with MathJax. Learn more.

$$\widetilde \varepsilon (\omega ) = {\varepsilon _0} + {\varepsilon _0}\widetilde \chi (\omega ),$$
$$\widetilde \chi (\omega ) = \left[ {\begin{array}{{ccc}} {{\varepsilon_{xx}}(\omega )}&{{\varepsilon_{xy}}(\omega )}&{{\varepsilon_{xz}}(\omega )}\\ {{\varepsilon_{yx}}(\omega )}&{{\varepsilon_{yy}}(\omega )}&{{\varepsilon_{yz}}(\omega )}\\ {{\varepsilon_{zx}}(\omega )}&{{\varepsilon_{zy}}(\omega )}&{{\varepsilon_{zz}}(\omega )} \end{array}} \right],$$
$$j\omega \widetilde \varepsilon \textrm{(}\omega \textrm{)}{\boldsymbol{E}}\textrm{(}\omega \textrm{)} = \nabla \times {\boldsymbol{H}}\textrm{(}\omega \textrm{)},$$
$$j\omega {\mu _0}{\boldsymbol{H}}\textrm{(}\omega \textrm{)} = \nabla \times {\boldsymbol{E}}\textrm{(}\omega \textrm{)}.$$
$$j\omega \widetilde \varepsilon \textrm{(}\omega \textrm{)}{\boldsymbol{E}}\textrm{(}\omega \textrm{)} = j\omega {\varepsilon _0}{\boldsymbol{E}}\textrm{(}\omega \textrm{) + }{\boldsymbol{J}}\textrm{(}\omega \textrm{),}$$
$${J_x}(\omega ) = j\omega {\varepsilon _0}[{{\varepsilon_{xx}}{E_x}(\omega ) + {\varepsilon_{xy}}{E_y}(\omega ) + {\varepsilon_{xz}}{E_z}(\omega )} ],$$
$${J_y}(\omega ) = j\omega {\varepsilon _0}[{{\varepsilon_{yx}}{E_x}(\omega ) + {\varepsilon_{yy}}{E_y}(\omega ) + {\varepsilon_{yz}}{E_z}(\omega )} ],$$
$${J_z}(\omega ) = j\omega {\varepsilon _0}[{{\varepsilon_{zx}}{E_x}(\omega ) + {\varepsilon_{zy}}{E_y}(\omega ) + {\varepsilon_{zz}}{E_z}(\omega )} ].$$
$$\frac{{\partial {\boldsymbol{E}}}}{{\partial t}} + \frac{1}{{{\varepsilon _0}}}{\boldsymbol{J}} = ({{\boldsymbol{A}}_{12}} + {{\boldsymbol{B}}_{12}}){\boldsymbol{H}}, $$
$$\frac{{\partial {\boldsymbol{H}}}}{{\partial t}} = ({{\boldsymbol{A}}_{21}} + {{\boldsymbol{B}}_{21}}){\boldsymbol{E}}, $$
$${{\boldsymbol{A}}_{12}} = \left[ {\begin{array}{{ccc}} 0&0&{\frac{1}{{{\varepsilon _0}}}\frac{\partial }{{\partial y}}} \\ {\frac{1}{{{\varepsilon _0}}}\frac{\partial }{{\partial z}}}&0&0 \\ 0&{\frac{1}{{{\varepsilon _0}}}\frac{\partial }{{\partial x}}}&0 \end{array}} \right],\,{{\boldsymbol{B}}_{12}} = \left[ {\begin{array}{{ccc}} 0&{\frac{{ - 1}}{{{\varepsilon _0}}}\frac{\partial }{{\partial z}}}&0 \\ 0&0&{\frac{{ - 1}}{{{\varepsilon _0}}}\frac{\partial }{{\partial x}}} \\ {\frac{{ - 1}}{{{\varepsilon _0}}}\frac{\partial }{{\partial y}}}&0&0 \end{array}} \right]$$
$${{\boldsymbol{A}}_{21}} = \left[ {\begin{array}{{ccc}} 0&{\frac{1}{{{\mu _0}}}\frac{\partial }{{\partial z}}}&0 \\ 0&0&{\frac{1}{{{\mu _0}}}\frac{\partial }{{\partial x}}} \\ {\frac{1}{{{\mu _0}}}\frac{\partial }{{\partial y}}}&0&0 \end{array}} \right],\,{{\boldsymbol{B}}_{21}} = \left[ {\begin{array}{{ccc}} 0&0&{\frac{{ - 1}}{{{\mu _0}}}\frac{\partial }{{\partial y}}} \\ {\frac{{ - 1}}{{{\mu _0}}}\frac{\partial }{{\partial z}}}&0&0 \\ 0&{\frac{{ - 1}}{{{\mu_0}}}\frac{\partial }{{\partial x}}}&0 \end{array}} \right]$$
$$({{\boldsymbol{I}}_{6 \times 6}} - \frac{{\Delta t}}{2}{\boldsymbol{A}}){{\boldsymbol{u}}^{n + 1/2}} = ({{\boldsymbol{I}}_{6 \times 6}} + \frac{{\Delta t}}{2}{\boldsymbol{A}}){{\boldsymbol{u}}^n} - \Delta t{{\boldsymbol{S}}^n}$$
$$({{\boldsymbol{I}}_{6 \times 6}} - \frac{{\Delta t}}{2}{\boldsymbol{B}}){{\boldsymbol{u}}^{n + 1}} = ({{\boldsymbol{I}}_{6 \times 6}} + \frac{{\Delta t}}{2}{\boldsymbol{B}}){{\boldsymbol{u}}^{n + 1/2}}, $$
$${\boldsymbol{A}} = \left[ {\begin{array}{{cc}} {{\mathbf{0}_{3 \times 3}}}&{{{\boldsymbol{A}}_{12}}} \\ {{{\boldsymbol{A}}_{21}}}&{{\mathbf{0}_{3 \times 3}}} \end{array}} \right],\,{\boldsymbol{B}} = \left[ {\begin{array}{{cc}} {{\mathbf{0}_{3 \times 3}}}&{{{\boldsymbol{B}}_{12}}} \\ {{{\boldsymbol{B}}_{21}}}&{{\mathbf{0}_{3 \times 3}}} \end{array}} \right].$$
$$({{\boldsymbol{I}}_{6 \times 6}} - \frac{{\Delta t}}{2}{\boldsymbol{B}}){{\boldsymbol{u}}^0} = {\boldsymbol{u}}_c^0, $$
$${\boldsymbol{u}}_c^{n + 1} = ({{\boldsymbol{I}}_{6 \times 6}} - \frac{{\Delta t}}{2}{\boldsymbol{B}}){{\boldsymbol{u}}^{n + 1}}. $$
$${\boldsymbol{u}}_c^{n + 1/2} = ({{\boldsymbol{I}}_{6 \times 6}} - \frac{{\Delta t}}{2}{\boldsymbol{A}}){{\boldsymbol{u}}^{n + 1/2}}, $$
$${\boldsymbol{u}}_c^{n + 1/2} = ({{\boldsymbol{I}}_{6 \times 6}} + \frac{{\Delta t}}{2}{\boldsymbol{A}}){{\boldsymbol{u}}^n} - \Delta t{{\boldsymbol{S}}^n}. $$
$${{\boldsymbol{u}}^n} = {({{\boldsymbol{I}}_{6 \times 6}} - \frac{{\Delta t}}{2}{\boldsymbol{B}})^{ - 1}}{\boldsymbol{u}}_c^n, $$
$$\begin{aligned} {u}_c^{\boldsymbol{n} + 1/2} &= ({{\boldsymbol{I}}_{6 \times 6}} + \frac{{\Delta t}}{2}{\boldsymbol{A}}){({{\boldsymbol{I}}_{6 \times 6}} - \frac{{\Delta t}}{2}{\boldsymbol{B}})^{ - 1}}{\boldsymbol{u}}_c^n - \Delta t{{\boldsymbol{S}}^n}\\& = ({{\boldsymbol I}_{6 \times 6}} + \frac{{\Delta t}}{2}{\boldsymbol{A}})({{\boldsymbol{I}}_{6 \times 6}} + \frac{{\Delta t}}{2}{\boldsymbol{B}}){\left[ {({{I}_{6 \times 6}} - \frac{{\Delta t}}{2}{\boldsymbol{B}})({{\boldsymbol{I}}_{6 \times 6}} + \frac{{\Delta t}}{2}{\boldsymbol{B}})} \right]^{ - 1}}{\boldsymbol{u}}_c^n - \Delta t{{\boldsymbol{S}}^n} \end{aligned}.$$
$$\begin{array}{c} \left[ {\begin{array}{{c}} {{\boldsymbol E}_c^{n + 1/2}}\\ {{\boldsymbol H}_c^{n + 1/2}} \end{array}} \right] = \textrm{(}{{\boldsymbol I}_{6 \times 6}} + \frac{{\Delta t}}{2}\left[ {\begin{array}{{cc}} {{{\mathbf 0}_{3 \times 3}}}&{{{\boldsymbol A}_{12}}}\\ {{{\boldsymbol A}_{21}}}&{{{\mathbf 0}_{3 \times 3}}} \end{array}} \right] + \frac{{\Delta t}}{2}\left[ {\begin{array}{{cc}} {{{\mathbf 0}_{3 \times 3}}}&{{{B}_{12}}}\\ {{{\boldsymbol B}_{21}}}&{{{\mathbf 0}_{3 \times 3}}} \end{array}} \right] + \frac{{\Delta {t^2}}}{4}\left[ {\begin{array}{{cc}} {{{\boldsymbol A}_{12}}{{\boldsymbol B}_{21}}}&{{{\mathbf 0}_{3 \times 3}}}\\ {{{\mathbf 0}_{3 \times 3}}}&{{{\boldsymbol A}_{21}}{{\boldsymbol B}_{12}}} \end{array}} \right])\\ {\textrm{(}{{\boldsymbol I}_{6 \times 6}} - \frac{{\Delta {t^2}}}{4}\left[ {\begin{array}{{cc}} {{{\boldsymbol B}_{12}}{{\boldsymbol B}_{21}}}&{{{\mathbf 0}_{3 \times 3}}}\\ {{{\mathbf 0}_{3 \times 3}}}&{{{\boldsymbol B}_{21}}{{\boldsymbol B}_{12}}} \end{array}} \right])^{ - 1}}\left[ {\begin{array}{{c}} {{\boldsymbol E}_c^n}\\ {{\boldsymbol H}_c^n} \end{array}} \right] - \frac{{\Delta t}}{{{\varepsilon _0}}}\left[ {\begin{array}{{c}} {{\boldsymbol J}_c^n}\\ {{{\mathbf 0}_{1 \times 3}}} \end{array}} \right] \end{array}, $$
$$\begin{aligned} {\boldsymbol E}_c^{n + 1/2} &= ({{\boldsymbol I}_{3 \times 3}} + \frac{{\Delta {t^2}}}{4}{{\boldsymbol A}_{12}}{{\boldsymbol B}_{21}}){({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol B}_{12}}{{\boldsymbol B}_{21}})^{ - 1}}{\boldsymbol E}_c^n\\& + \frac{{\Delta t}}{2}({{\boldsymbol A}_{12}} + {{\boldsymbol B}_{12}}){({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol B}_{21}}{{\boldsymbol B}_{12}})^{ - 1}}{\boldsymbol H}_c^n - \frac{{\Delta t}}{{{\varepsilon _0}}}{\boldsymbol J}_c^n \end{aligned}, $$
$$\begin{array}{c} {\boldsymbol E}_c^{n - 1/2} = ({{\boldsymbol I}_{3 \times 3}} + \frac{{\Delta {t^2}}}{4}{{\boldsymbol A}_{12}}{{\boldsymbol B}_{21}}){({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol B}_{12}}{{\boldsymbol B}_{21}})^{ - 1}}{\boldsymbol E}_c^n\\ - \frac{{\Delta t}}{2}({{\boldsymbol A}_{12}} + {{\boldsymbol B}_{12}}){({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol B}_{21}}{{\boldsymbol B}_{12}})^{ - 1}}{\boldsymbol H}_c^n \end{array}, $$
$${\boldsymbol E}_c^{n + 1/2} - {\boldsymbol E}_c^{n - 1/2} = \Delta t({{\boldsymbol A}_{12}} + {{\boldsymbol B}_{12}}){({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol B}_{21}}{{\boldsymbol B}_{12}})^{ - 1}}{\boldsymbol H}_c^n - \frac{{\Delta t}}{{{\varepsilon _0}}}{\boldsymbol J}_c^n. $$
$${\boldsymbol H}_c^{n + 1} - {\boldsymbol H}_c^n = \Delta t({{\boldsymbol A}_{21}} + {{\boldsymbol B}_{21}}){({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol A}_{21}}{{\boldsymbol A}_{12}})^{ - 1}}{\boldsymbol E}_c^{n + 1/2}. $$
$$({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol B}_{21}}{{B}_{12}}){\boldsymbol h}_c^n = {\boldsymbol H}_c^n$$
$${\boldsymbol E}_c^{n + 1/2} = {\boldsymbol E}_c^{n - 1/2} + \Delta t({{\boldsymbol A}_{12}} + {{\boldsymbol B}_{12}}){\boldsymbol h}_c^n - \frac{{\Delta t}}{{{\varepsilon _0}}}{\boldsymbol J}_c^n$$
$$({{\boldsymbol I}_{3 \times 3}} - \frac{{\Delta {t^2}}}{4}{{\boldsymbol A}_{12}}{{\boldsymbol A}_{21}}){\boldsymbol e}_c^{n + 1/2} = {\boldsymbol E}_c^{n + 1/2}$$
$${\boldsymbol H}_c^{n + 1} = {\boldsymbol H}_c^n + \Delta t({{\boldsymbol A}_{21}} + {{\boldsymbol B}_{21}}){\boldsymbol e}_c^{n + 1/2}. $$
$$\widetilde \chi (\omega ) = \left[ {\begin{array}{{ccc}} {{\varepsilon_d}(\omega )}&{j{\varepsilon_g}(\omega )}&0\\ { - j{\varepsilon_g}(\omega )}&{{\varepsilon_d}(\omega )}&0\\ 0&0&{{\varepsilon_{zz}}(\omega )} \end{array}} \right], $$
$${{\boldsymbol J}_c} = \left[ {\begin{array}{{ccc}} {\frac{{j\omega {\varepsilon_0}\omega_p^2 + {\varepsilon_0}\omega_p^2{v_c}}}{{({v_c} + \omega_b^2) + 2j\omega {v_c} + {{(j\omega )}^2}}}}&{\frac{{ - {\varepsilon_0}\omega_p^2{\omega_b}}}{{({v_c} + \omega_b^2) + 2j\omega {v_c} + {{(j\omega )}^2}}}}&0\\ {\frac{{{\varepsilon_0}\omega_p^2{\omega_b}}}{{({v_c} + \omega_b^2) + 2j\omega {v_c} + {{(j\omega )}^2}}}}&{\frac{{j\omega {\varepsilon_0}\omega_p^2 + {\varepsilon_0}\omega_p^2{v_c}}}{{({v_c} + \omega_b^2) + 2j\omega {v_c} + {{(j\omega )}^2}}}}&0\\ 0&0&{\frac{{ - j\omega {\varepsilon_0}\omega_p^2}}{{\omega (\omega - j{v_c})}}} \end{array}} \right]\left[ {\begin{array}{{c}} {{E_{cx}}}\\ {{E_{cy}}}\\ {{E_{cz}}} \end{array}} \right]. $$
$${{\boldsymbol J}_c} = \left[ {\begin{array}{{c}} {{J_{xx}}(\omega ) + {J_{xy}}(\omega ) + {J_{xz}}(\omega )}\\ {{J_{yx}}(\omega ) + {J_{yy}}(\omega ) + {J_{yz}}(\omega )}\\ {{J_{zx}}(\omega ) + {J_{zy}}(\omega ) + {J_{zz}}(\omega )} \end{array}} \right], $$
$${{\boldsymbol J}_p}\textrm{(}\omega \textrm{)} = {\varepsilon _0}\frac{{{a_1}j\omega + {a_0}}}{{{b_2}{{(j\omega )}^2} + {b_1}j\omega + {b_0}}}{{\boldsymbol E}_c}\textrm{(}\omega \textrm{)}, $$
$${b_2}\frac{{{d^2}{{\boldsymbol J}_p}}}{{d{t^2}}} + {b_1}\frac{{d{{\boldsymbol J}_p}}}{{dt}} + {b_0}{{\boldsymbol J}_p} = {\varepsilon _0}{a_1}\frac{{d{{\boldsymbol E}_c}}}{{dt}} + {\varepsilon _0}{a_0}{{\boldsymbol E}_c}. $$
$${b_0}{\boldsymbol J}_p^{n - 1} + {b_1}\frac{{{\boldsymbol J}_p^n - {\boldsymbol J}_p^{n - 2}}}{{2\Delta t}}\textrm{ + }{b_2}\frac{{{\boldsymbol J}_p^n - 2{\boldsymbol J}_p^{n - 1} + {\boldsymbol J}_p^{n - 2}}}{{\Delta {t^2}}} = {a_0}\frac{{{\boldsymbol E}_c^{n - 1/2} + {\boldsymbol E}_c^{n - 3/2}}}{2} + {a_1}\frac{{{\boldsymbol E}_c^{n - 1/2} - {\boldsymbol E}_c^{n - 3/2}}}{{\Delta t}}, $$
$${\boldsymbol J}_p^n = {d_1}{\boldsymbol J}_p^{n - 1} + {d_2}{\boldsymbol J}_p^{n - 2} + {d_3}{\boldsymbol E}_c^{n - 1/2} + {d_4}{\boldsymbol E}_c^{n - 3/2}, $$
$$J_{xx}^n = {f_1}J_{xx}^{n - 1} + {f_2}J_{xx}^{n - 2} + {f_3}E_{cx}^{n - 1/2} + {f_4}E_{cx}^{n - 3/2}, $$
$$J_{xy}^n = {f_1}J_{xy}^{n - 1} + {f_2}J_{xy}^{n - 2} - {f_5}E_{cy}^{n - 1/2} - {f_6}E_{cy}^{n - 3/2}$$
$$J_{xz}^n = 0, $$
$$\begin{aligned} J_{cx}^n &= {f_1}(J_{xx}^{n - 1} + J_{xy}^{n - 1} + J_{xz}^{n - 1}) + {f_2}(J_{xx}^{n - 2} + J_{xy}^{n - 2} + J_{xz}^{n - 2})\\& + {f_3}E_{cx}^{n - 1/2} + {f_4}E_{cx}^{n - 3/2} - {f_5}E_{cy}^{n - 1/2} - {f_6}E_{cy}^{n - 3/2}\\& = {f_1}J_{cx}^{n - 1} + {f_2}J_{cx}^{n - 2} + {f_3}E_{cx}^{n - 1/2} + {f_4}E_{cx}^{n - 3/2} - {f_5}E_{cy}^{n - 1/2} - {f_6}E_{cy}^{n - 3/2} \end{aligned}.$$
$$\begin{aligned} {E_{cy}}(i + 1/2,j,k) &= ({E_{cy}}(i,j + 1/2,k) + {E_{cy}}(i,j - 1/2,k)\\& + {E_{cy}}(i + 1,j + 1/2,k) + {E_{cy}}(i + 1,j - 1/2,k))/4 \end{aligned}.$$
$$J_{cy}^n = {f_1}J_{cy}^{n - 1} + {f_2}J_{cy}^{n - 2} + {f_3}E_{cy}^{n - 1/2} + {f_4}E_{cy}^{n - 3/2} + {f_5}E_{cx}^{n - 1/2} + {f_6}E_{cx}^{n - 3/2}, $$
$$J_{cz}^n = {l_1}J_{cz}^{n - 1} + {l_2}J_{cz}^{n - 2} + {l_3}E_{cz}^{n - 1/2} - {l_4}E_{cz}^{n - 3/2}, $$
$$\begin{aligned} {E_{cx}}(i,j + 1/2,k) &= ({E_{cx}}(i + 1/2,j,k) + {E_{cx}}(i - 1/2,j,k)\\& + {E_{cx}}(i + 1/2,j + 1,k) + {E_{cx}}(i - 1/2,j + 1,k))/4 \end{aligned}.$$
$${\widetilde {\sigma }_s}(\omega ) = \left[ {\begin{array}{{ccc}} {{\sigma_{xx}}(\omega )}&{{\sigma_{xy}}(\omega )}&0\\ {{\sigma_{yx}}(\omega )}&{{\sigma_{yy}}(\omega )}&0\\ 0&0&0 \end{array}} \right] = \left[ {\begin{array}{{ccc}} {{\sigma_d}(\omega )}&{ - {\sigma_g}(\omega )}&0\\ {{\sigma_g}(\omega )}&{{\sigma_d}(\omega )}&0\\ 0&0&0 \end{array}} \right], $$
$${\sigma _d}(\omega ) = \frac{{j\omega {\sigma _0} + v{\sigma _0}}}{{(\omega _c^2 + {v^2}) + 2j\omega v + {{(j\omega )}^2}}}, $$
$${\sigma _g}(\omega ) = \frac{{ - {\sigma _0}{\omega _c}}}{{(\omega _c^2 + {v^2}) + 2j\omega v + {{(j\omega )}^2}}}, $$
$${\sigma _0} = \frac{{2{e^2}{k_B}T}}{{\pi {\hbar ^2}}}\ln (2\cosh \frac{{{\mu _c}}}{{2{k_B}T}}), $$
$${\omega _c} = e{B_0}v_F^2/{\mu _c}, $$
$$\begin{aligned}{{\boldsymbol J}_c}(\omega ) &= \frac{1}{{\Delta z}}\left[ {\begin{array}{{ccc}} {\frac{{{\sigma_0}j\omega + v{\sigma_0}}}{{(\omega_c^2 + {v^2}) + 2j\omega v + {{(j\omega )}^2}}}}&{\frac{{ - {\omega_c}{\sigma_0}}}{{(\omega_c^2 + {v^2}) + 2j\omega v + {{(j\omega )}^2}}}}&0\\ {\frac{{{\omega_c}{\sigma_0}}}{{(\omega_c^2 + {v^2}) + 2j\omega v + {{(j\omega )}^2}}}}&{\frac{{{\sigma_0}j\omega + v{\sigma_0}}}{{(\omega_c^2 + {v^2}) + 2j\omega v + {{(j\omega )}^2}}}}&0\\ 0&0&0 \end{array}} \right]\left[ {\begin{array}{{c}} {{E_{cx}}}\\ {{E_{cy}}}\\ {{E_{cz}}} \end{array}} \right]\\& = \left[ {\begin{array}{{ccc}} {{J_{xx}}(\omega )}&{{J_{xy}}(\omega )}&0\\ {{J_{yx}}(\omega )}&{{J_{yy}}(\omega )}&0\\ 0&0&0 \end{array}} \right] \end{aligned}.$$
$$J_{cx}^n = {f_1}J_{cx}^{n - 1} + {f_2}J_{cx}^{n - 2} + {f_3}E_{cx}^{n - 1/2} + {f_4}E_{cx}^{n - 3/2} - {f_5}E_{cy}^{n - 1/2} - {f_6}E_{cy}^{n - 3/2}, $$
$$J_{cy}^n = {f_1}J_{cy}^{n - 1} + {f_2}J_{cy}^{n - 2} + {f_3}E_{cy}^{n - 1/2} + {f_4}E_{cy}^{n - 3/2} + {f_5}E_{cx}^{n - 1/2} + {f_6}E_{cx}^{n - 3/2}, $$
$${{\boldsymbol M}_L}\left[ {\begin{array}{{c}} {{\boldsymbol E}_c^{n + 1/2}}\\ {{\boldsymbol H}_c^{n + 1}}\\ {{\boldsymbol J}_c^n} \end{array}} \right] = {{\boldsymbol M}_{R1}}\left[ {\begin{array}{{c}} {{\boldsymbol E}_c^{n - 1/2}}\\ {{\boldsymbol H}_c^n}\\ {{\boldsymbol J}_c^{n - 1}} \end{array}} \right] + {{\boldsymbol M}_{R2}}\left[ {\begin{array}{{c}} {{\boldsymbol E}_c^{n - 3/2}}\\ {{\boldsymbol H}_c^{n - 1}}\\ {{\boldsymbol J}_c^{n - 2}} \end{array}} \right], $$
$${{\boldsymbol M}_L} = \left[ {\begin{array}{{ccccccccc}} 1&0&0&0&0&0&a&0&0\\ 0&1&0&0&0&0&0&a&0\\ 0&0&1&0&0&0&0&0&0\\ 0&{ - b{D_z}{F_z}}&{b{D_y}{F_x}}&1&0&0&0&0&0\\ {b{D_z}{F_y}}&0&{ - b{D_x}{F_x}}&0&1&0&0&0&0\\ { - b{D_y}{F_y}}&{b{D_x}{F_z}}&0&0&0&1&0&0&0\\ 0&0&0&0&0&0&1&0&0\\ 0&0&0&0&0&0&0&1&0\\ 0&0&0&0&0&0&0&0&1 \end{array}} \right], $$
$${{\boldsymbol M}_{R1}} = \left[ {\begin{array}{{ccccccccc}} 1&0&0&0&{ - a{D_z}{F_z}}&{a{D_y}{F_x}}&0&0&0\\ 0&1&0&{a{D_z}{F_y}}&0&{ - a{D_x}{F_x}}&0&0&0\\ 0&0&1&{ - a{D_y}{F_y}}&{a{D_x}{F_z}}&0&0&0&0\\ 0&0&0&1&0&0&0&0&0\\ 0&0&0&0&1&0&0&0&0\\ 0&0&0&0&0&1&0&0&0\\ {{f_3}}&{ - {f_5}}&0&0&0&0&{{f_1}}&0&0\\ {{f_5}}&{{f_3}}&0&0&0&0&0&{{f_1}}&0\\ 0&0&0&0&0&0&0&0&0 \end{array}} \right], $$
$${{\boldsymbol M}_{R2}} = \left[ {\begin{array}{{ccccccccc}} 0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0\\ {{f_4}}&{ - {f_6}}&0&0&0&0&{{f_2}}&0&0\\ {{f_6}}&{{f_4}}&0&0&0&0&0&{{f_2}}&0\\ 0&0&0&0&0&0&0&0&0 \end{array}} \right], $$
$$\phi _{i,j,k}^n = {\varphi _\phi }{\zeta ^n}\exp [\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over j} ({k_x}i\Delta x + {k_y}j\Delta y + {k_z}k\Delta z)], $$
$${D_p}\phi _{i,j,k}^n = \frac{{\phi _{i,j,k}^n(p + \frac{{\Delta p}}{2}) - \phi _{i,j,k}^n(p - \frac{{\Delta p}}{2})}}{{\Delta p}} = {\sigma _p}\phi _{i,j,k}^n, $$
$$({{\boldsymbol M}_{R1}} + {{{{\boldsymbol M}_{R2}}} / \zeta } - \zeta {{\boldsymbol M}_L})\left[ {\begin{array}{{c}} {{\boldsymbol E}_c^{n - 1/2}}\\ {{\boldsymbol H}_c^n}\\ {{\boldsymbol J}_c^{n - 1}} \end{array}} \right] = {{\boldsymbol M}}\left[ {\begin{array}{{c}} {{\boldsymbol E}_c^{n - 1/2}}\\ {{\boldsymbol H}_c^n}\\ {{\boldsymbol J}_c^{n - 1}} \end{array}} \right] = 0,$$
$$\scalebox{0.8}{$\displaystyle{\boldsymbol M} = \left[ {\begin{array}{{ccccccccc}} {1 - \zeta }&0&0&0&{ - a{\sigma_z}{F_z}}&{a{\sigma_y}{F_x}}&{ - a\zeta }&0&0\\ 0&{1 - \zeta }&0&{a{\sigma_z}{F_y}}&0&{ - a{\sigma_x}{F_x}}&0&{ - a\zeta }&0\\ 0&0&{1 - \zeta }&{ - a{\sigma_y}{F_y}}&{a{\sigma_x}{F_z}}&0&0&0&0\\ 0&{b\zeta {\sigma_z}F}&{ - b\zeta {\sigma_y}{F_x}}&{1 - \zeta }&0&0&0&0&0\\ { - b\zeta {\sigma_z}{F_y}}&0&{b\zeta {\sigma_x}{F_x}}&0&{1 - \zeta }&0&0&0&0\\ {b\zeta {\sigma_y}{F_y}}&{ - b\zeta {\sigma_x}{F_z}}&0&0&0&{1 - \zeta }&0&0&0\\ {{f_3} + {{{f_4}} / \zeta }}&{ - {f_5} - {{{f_6}} / \zeta }}&0&0&0&0&{{f_1} + {{{f_2}} / {\zeta - \zeta }}}&0&0\\ {{f_5} + {{{f_6}} / \zeta }}&{{f_3} + {{{f_4}} / \zeta }}&0&0&0&0&0&{{f_1} + {{{f_2}} / {\zeta - \zeta }}}&0\\ 0&0&0&0&0&0&0&0&{ - \zeta } \end{array}} \right]$}$$
$${E_r}(f) = \frac{{|{{T_{num}}(f) - {T_{Analytical}}(f)} |}}{{\max |{{T_{Analytical}}(f)} |}} \times 100\%$$
$$\textrm{pulse}(t) = \cos (2\pi {f_0}t)\exp [ - 4\pi {(t - {t_0})^2}/{\tau ^2}], $$
Select as filters
Select Topics Cancel
Top
       
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.
Privacy | Terms of Use