Time-Domain Studies of General Dispersive Anisotropic Media by the Complex-Conjugate Pole–Residue Pairs Model

: A novel ﬁnite-difference time-domain formulation for the modeling of general anisotropic dispersive media is introduced in this work. The method accounts for fully anisotropic electric or magnetic materials with all elements of the permittivity and permeability tensors being non-zero. In addition, each element shows an arbitrary frequency dispersion described by the complex-conjugate pole–residue pairs model. The efﬁciency of the technique is demonstrated in benchmark numerical examples involving electromagnetic wave propagation through magnetized plasma, nematic liquid crystals and ferrites.


Introduction
The finite-difference time-domain (FDTD) method [1] is widely recognized as one of the dominant tools in computational electromagnetics, particularly in the study of the broadband spectral response of electromagnetic structures. In such large frequency intervals, most materials are dispersive; namely, they exhibit a frequency-dependent variation of their electrical permittivity and magnetic permeability. From the early stages of its development, the FDTD was demonstrated to be capable of accurately simulating the dispersive properties of most common materials, both dielectric and metallic [2].
In several cases involving advanced natural or artificial materials, it is not only their dispersive nature that has to be taken into account but also their anisotropy. This problem has been addressed in various FDTD techniques, which have been proposed for the numerical simulation of certain types of both anisotropic and dispersive media [3][4][5][6][7][8][9].
As concerns the FDTD modeling of dispersive materials, the complex-conjugate poleresidue pairs (CCPR) model, a generalization of the traditional Drude and Lorentz models, has been proven to provide a more accurate description of material dispersion [10,11] with reduced memory requirements and similar conditions for numerical stability [12]. Recently, we have demonstrated how the CCPR can be combined with electric anisotropy for the accurate investigation of nematic liquid crystals (NLC) [13]. The dispersion of the ordinary and extraordinary NLC indices were described by the CCPR model, and the FDTD scheme accounted for both the dispersion and the NLC anisotropy. CCPR modeling has been also combined with implicit FDTD formulations-in particular, the leapfrog alternating-direction implicit FDTD-confirming the robustness of the approach for the accurate description of dispersive materials [14].
In the present work, the CCPR dispersion model is incorporated into an explicit FDTD, which accounts for a generalized form of the permittivity or permeability tensor. All elements of the electric and magnetic tensors are allowed to obtain non-zero values. Moreover, they are allowed to independently exhibit frequency dispersion following the CCPR model with an arbitrary number of terms and are capable of capturing the behavior of dispersive material over broad spectral ranges. The FDTD implementation algorithm is outlined, and the associated generalized update expressions are obtained. The efficiency of the technique is corroborated in benchmark examples studying electromagnetic wave propagation through three types of dispersive and anisotropic materials: (i) magnetized plasma, (ii) NLC and (iii) ferrites. The proposed FDTD scheme provides extensive versatility in the modeling of the arbitrary dispersive characteristics of the tensor elements of anisotropic materials and a robust framework for the study of electromagnetic systems that involve dispersive anisotropic media.

Formulation
The time-harmonic Maxwell's equations in a linear, source-free, electrically and magnetically anisotropic dispersive medium with e jωt time dependence are as follows: where the tilde on the fields denotes the frequency domain, and the frequency-dependent relative permittivity tensor ε(ω) is described by the following general expression: Similarly, the relative permeability tensor µ(ω) is given by In this work, we adopt the CCPR model for the description of the dispersive elements ε αβ (ω), µ αβ (ω). We explicitly present the formulation for electric anisotropic dispersive media, according to which the complex relative permittivity of the tensor elements is given by [10] where ε ∞,αβ is the relative permittivity at infinite frequency, σ e αβ is the static conductivity, α, β denote x, y and z and * represents the complex conjugate. The tensor elements µ αβ (ω) in the case of anisotropic magnetic materials are described by the same model, and the formulation follows the same steps described below, with Equations (5)- (15) referring to the magnetic instead of the electric field.
By the proper selection of the coefficients in Equation (4), the CCPR model incorporates all the standard dispersion models (see Table 1), while also several more advanced models that have been recently proposed for metals, such as the Drude-critical points [15] and modified Lorentz [16]. As a consequence, it can also handle other dispersion models, which can be approximated with those in Table 1, such the Cole-Cole and Havriliak-Negami models that characterize certain biological media and polymers [17,18]. Furthermore, one of the strengths of the proposed formulation lies in the ability of the CCPR model to accurately fit experimental data on material permittivity using the vector fitting technique [19,20] or other related algorithms [21]. Table 1. CCPR parameters for the description of various known dispersion models.

Model Expression Parameters of the CCPR Model: σ, (c p , a p )
Drude are auxiliary variables. The corresponding variablesQ pxy ,Q pxy ,Q pxz andQ pxz are obtained in a similar manner. Next, we convert Equation (6) into the time domain: Since the electric field E x is a real quantity in the time domain, we conclude from Equations (7a) and (7b) that Q pxx = Q * pxx , and as a result, only the variable Q pxx is needed for the simulation. Similar conclusions can be obtained for the other variables.
Subsequently, we convert Equation (5) into the time domain: and using the formula z + z * = 2 (z), with (·) denoting the real part, we obtain We discretize (7a) at time (n + 1/2)∆t and obtain the updated equation of the Q pxx : Similar update expressions can be derived for Q pxy and Q pxz . We also discretize Equation (9) at time (n + 1/2)∆t, and using the update expressions of Q pxx , Q pxy and Q pxz , we obtain where α xx , α xy and α xz are constants defined as and γ n x changes at every time step according to Working in a similar way, two more equations are obtained and a 3 × 3 linear system is formed: The update equation of E x is obtained by the solution of Equation (14) using Cramer's rule: where D is the determinant of the matrix with elements α ij as in the l.h.s. of Equation (14). The update equations of E y and E z are obtained in a similar manner. Finally, the field updating procedure in each iteration of the algorithm is summarized in the following steps:

1.
Update H x , H y and H z ; 2.
Store the current values of E x , E y and E z ; 3.
Update E x , E y and E z ; 4.

Applications and Numerical Results
In this section, we study three benchmark problems involving materials which are both anisotropic and dispersive. Although the problems address different spectral ranges and dispersion relations, they can be successfully handled with the proposed numerical framework with limited modifications of the computer code. In all examined cases, we have opted to present the spectral properties of the transmitted and reflected waves, as calculated by the Fourier transform of the FDTD simulated time-domain signals, in order to (i) provide a direct reference to published results, where available; (ii) highlight some physical properties of the investigated materials, such as bandgaps; and (iii) provide straightforward comparison to the benchmark solutions, which were calculated by frequency-domain methods.

Propagation in Magnetized Plasma
We first study wave propagation in anisotropic magnetized plasma in which the direction of propagation is along the z-axis and the external DC magnetic field is H = H DC z as in [22] (chap. 8). In this case, the nonzero elements of the permittivity tensor, assuming the e jωt time convention, are given by where ω p is the plasma frequency, v is the electron collision rate and ω b is the cyclotron frequency. In view of the tensor elements of Equation (16), the parameters of the CCPR model are derived through the simple algebraic manipulations shown in Table 2.

Tensor Element
Parameters: σ, (c p , a p ) We study the wave propagation of a Gaussian pulse normally incident to a magnetized plasma slab of 9 mm in thickness, assuming that the static biasing magnetic is parallel to the direction of propagation. The plasma parameters are ω p = 2π × 50 × 10 9 rad/s, ω b = 3 × 10 11 rad/s, and v c = 2 × 10 10 Hz. The FDTD cell size is ∆z = 75 µm and the time step is chosen as ∆t = Q∆z/c 0 , where Q = 0.3 is the Courant number and c 0 is the velocity of light in free space. The Gaussian pulse is excited from the left side of the computational domain at the x component of the electric field. The domain is truncated by a 12 cell convolution perfectly matched layer (CPML) [23]. We calculate the reflection (Γ) and transmission (T) coefficients for right circularly polarized (RCP) and left circularized polarized (LCP), defined as where E xt and E yt are the x and y components of the transmitted electric field in the frequency domain, respectively, and E xr and E yr are likewise the reflected components. E xi is the x component of the incident wave in the frequency domain. In Figure 1, we validate the proposed FDTD method using as a reference the semi-analytical solution calculated by the Berreman 4 × 4 method [24]. Excellent agreement is observed in a very broad spectral range (0-100 GHz).

Terahertz Wave Propagation through a Nematic Liquid Crystal Cell
For NLC materials, the relative permittivity tensor is given by ε o + ∆ε cos 2 θ cos 2 φ ∆ε cos 2 θ sin φ cos φ ∆ε sin θ cos θ cos φ ∆ε cos 2 θ sin φ cos φ ε o + ∆ε cos 2 θ sin 2 φ ∆ε sin θ cos θ sin φ ∆ε sin θ cos θ cos φ ∆ε sin θ cos θ sin φ where θ and φ are the tilt and twist angles of the spatially varying LC molecular orientation, respectively, and ε o and ε e are the frequency-dependent ordinary and extraordinary dielectric permittivities, with ∆ε = ε e − ε o . The tilt angle is defined as the angle between the nematic director, a unit vector parallel to the NLC optical axis, and the x-y plane, whereas the twist angle is measured between the x axis and the projection of the director on the x-y plane. The ordinary and extraordinary dielectric permittivities are described as the sum of CCPR terms according to Equation (4). We demonstrate the efficiency of the proposed FDTD formulation for the case of THz wave propagation through an LC cell. The NLC material is 5CB, whose material dispersion has been described via the modified Lorentz model [25], and the equivalent CCPR parameters are obtained using Table 1. The real and imaginary parts of the extraordinary and ordinary relative permittivity (ε o,e = ε − jε ) in the considered spectral range are shown in Figure 2a. A layer of 0.5 mm is embedded in a substrate material with a refractive index of 1.5. An x polarized plane wave impinges on the slab, and the transmittance after the slab is calculated. The FDTD cell size is chosen to be equal to 0.5 µm, and the computational domain is terminated by a 12-cell CPML. The FDTD time step is set to ∆t = 0.3∆z/c 0 and the simulation is run for 60,000 time steps. A direct comparison of the transmittance of the structure obtained by the proposed FDTD scheme and the reference matrix technique is provided in Figure 2b for indicative values of the θ and φ angles. Excellent agreement is observed.

Propagation in Ferrites
As a last example, we investigate the case of wave propagation through a ferrite slab. The nonzero elements of the permeability tensor are given by [22] µ xx (ω) = µ yy (ω) = µ(ω) (19a) where it is assumed that the magnetic bias is in the z-direction. Loss is accounted for via the complex resonant frequency; namely, ω 0 → ω 0 + jωa. The CCPR parameters for this case are given in Table 3 and they are described by modified Lorentz relations.
A computational domain with 5000 cells with a cell size of 15 µm is assumed, which is backed with a 12-cell CPML. A ferrite slab of thickness 15 mm is surrounded by air, and an amplitude-modulated Gaussian pulse with a frequency range from 500 MHz to 14 GHz is excited and impinges upon the slab. The gyromagnetic response frequency is ω 0 = 2π × 4 × 10 9 rad/s, while the saturation magnetization frequency is ω M = 2π × 5.6 × 10 9 rad/s with a = 0.05 and ε r = 10 [26]. The transmission and reflection coefficients for RCP and LCP waves are calculated using the proposed formulation and compared to the finite element method (FEM) solution implemented in the software COMSOL Multiphysics. A direct comparison of the results is provided in Figure 3, once again demonstrating the efficiency of the proposed FDTD algorithm.

Conclusions
An FDTD formulation is introduced for the study of anisotropic media with generalized dispersive properties, described by the CCPR function. Both electric and magnetic anisotropy are considered with independently dispersive elements of the electric permittivity and magnetic permeability tensors, respectively. The validity of the proposed algorithm is demonstrated in three case examples involving magnetized plasma, liquid crystals and ferrites by comparing the FDTD results with those obtained with benchmark semi-analytic or numerical solutions.
Author Contributions: Conceptualization, K.P.P. and D.C.Z.; methodology, K.P.P.; software, K.P.P.; validation, K.P.P. and D.C.Z.; writing-original draft preparation, K.P.P.; writing-review and editing, D.C.Z.; supervision, D.C.Z. All authors have read and agreed to the published version of the manuscript.