Numerical Modeling of Electromagnetic Wave Propagation Through Bi-Isotropic Materials

1. Introduction

As part of the present wave of interest in electromagnetic propagation in complex materials (left-handed, chiral, nonreciprocal, nonlinear, etc.), a special mention may be made to bi-isotropic materials. After a period of relevance during the nineties [1-4], nowadays they enjoy a renewed interest, because of new fabrication techniques that allow the use of mature technologies (printed or integrated circuit) [5,6], and also the discovery of new properties associated to the bi-isotropy, such us the presence of negative refraction index [6-11]. In order to study electromagnetic wave propagation in such materials, different traditional numerical methods, both in time and frequency domain, have been extended so that they may deal with additional constitutive parameters, as well as their inherently dispersive. This last point is especially important when dealing with time-domain methods, in which the frequency dispersion must be translated to time-domain.

In this work, we want to present a review of different techniques developed during the last years:

• Time domain techniques: Finite Differences in the Time Domain (FDTD), Transmission Line Matrix (TLM) and Multiresolution in Time Domain (MRTD). All of them have been modified to treat with bi-isotropic media, and with several methods to include frequency dispersion.

• In the frequency domain, the extension of the coupled mode method to cover the electromagnetic wave propagation in closed structures containing these new media.

2. Bi-Isotropic materials. Constitutive Relations

In a general bi-isotropic material, the electric and magnetic fields are coupled. There are different ways of formulating its constitutive relations in frequency domain [1]. Among them, we use:

$$\begin{array}{c}\overrightarrow{D}=\epsilon \overrightarrow{E}+\xi \overrightarrow{H}\\ \overrightarrow{B}=\zeta \overrightarrow{E}+\mu \overrightarrow{H}\end{array}$$ E1

In which we assume the medium is linear. $\overrightarrow{D},\overrightarrow{E},\overrightarrow{B}$ , and $\overrightarrow{H}$ are the electromagnetic field vectors. ε and µ are the permittivity and permeability, ζ and ζ would be the cross-coupling parameters: the four parameters may depend on frequency (dispersive behavior). In a bi-anisotropic medium, they have tensor form, and may contain up to 36 scalar material parameters. ζ and ζ are usually expressed in the following way:

$$\begin{array}{c}\zeta =\frac{\chi +j\kappa }{c}\\ \xi =\frac{\chi -j\kappa }{c}\end{array}$$ E2

c being the light speed in vacuum. The parameter is known as “Pasteur” parameter, while χ is the “Tellegen” parameter; thanks to the introduction of c in the equation, both parameters are dimensionless; at the same time, in case of non-dispersive behavior, they are also real [1]. Substituting these values into Eq. 2, it leads to:

$$\begin{array}{c}\overrightarrow{D}=\epsilon \overrightarrow{E}+\frac{\chi -j\kappa }{c}\overrightarrow{H}\\ \overrightarrow{B}=\mu \overrightarrow{H}+\frac{\chi +j\kappa }{c}\overrightarrow{E}\end{array}$$ E3

Using these relations, instead of Eq. 1, the constitutive parameters have a more clear physical meaning: the Pasteur parameter describes the “chirality” of the medium, i.e. if the medium is “non chiral” (that is, it is equivalent to its mirror image), then κ= 0. In other case, the medium is said to be “chiral”. The Tellegen parameter represents the nonreciprocity of the medium. A medium may be nonreciprocal due to other reasons (for example, a magnetized ferrite), but only a cross coupling like the one described in Eq. 3 may be nonreciprocal and isotropic (that is “bi-isotropic”) at the same time [2].

While reciprocal chiral materials have been extensively researched during the last years, and different geometries and fabrication techniques have been designed [5,6], nonreciprocal materials have stayed mainly in the plane of theoretical discussion. In [12,13], Lakhtakia and Weiglhofer affirmed, using a covariance requirement deduced by Post (the so-called “Post constraint” [14]), that a bi-isotropic medium must be reciprocal, so the Tellegen parameter should always be 0. This opened a discussion on the theoretical possibility of the existence of such materials that lasted more than a decade, mainly between Lakhtakia and Weiglhofer [12,13,15-18] and Lindell, Sihvola, Tretyakov et al. [3,4,19-23]. Finally, Tretyakov et al. [24] built a “Tellegen particle” and verified its non-reciprocal behavior. At the same time, it has been found that some natural media (Fe₂TeO₆ and Cr₂O₃) though anisotropic, violate the Post constraint [25], and it may be demonstrated that only in uniform and unbounded media it may be stated that the Tellegen parameter vanishes [26-28], so “the Post constraint as a general dogma should be buried with all due honors” [29]. Other authors, nevertheless, consider the question still open [30].

An alternative formulation of Eq. 3 was proposed by Kong [31]:

$$\begin{array}{c}\overrightarrow{D}=\epsilon \overrightarrow{E}+\frac{1}{c}\left(\chi \overrightarrow{H}-\kappa \frac{d\overrightarrow{H}}{dt}\right)\\ \overrightarrow{B}=\mu \overrightarrow{H}+\frac{1}{c}\left(\chi \overrightarrow{E}-\kappa \frac{d\overrightarrow{E}}{dt}\right)\end{array}$$ E4

Note that this equation is formulated in time domain, instead of frequency domain, and that the chirality parameter κ is not equivalent to the one present in Eq. 3. In addition, this equation does not take into account the intrinsically dispersive behavior of this class of materials [1,6,9-11]. To deal with that dispersion, already in 1937, E.U. Condon [32] proposed a second-order resonant model for optical chiral materials. Though in 1999 Kristensson showed that such model is non-physical [33], it is still useful as an approximation: experimental measurements show that the chirality parameter follows indeed a second order-resonant behavior, that may be modeled as a combination of two orthogonal oscillators [34]. Electric permittivity and magnetic permeability also follow a second-order resonant behavior (Lorentz model [35]) with a resonant frequency virtually identical to ω ₀ [11, 36]

Electromagnetic waves in general dispersive bi-isotropic media show some interesting behavior [1]. First, the existence of the Pasteur parameter leads to:

• Optical/electromagnetic activity (OA): it is the turning of the plane of a linearly polarized electromagnetic wave about the direction of motion as the wave travels through the material. This is caused by the fact that a chiral material transmits right- and left- handed circularly polarized waves with different phase velocities.

• Optical/electromagnetic rotatory dispersion (ORD) is the variation of the optical/electromagnetic activity with frequency. Usually, that dispersion leads to an OA that changes its sign at a resonance frequency (Cotton effect).

• Circular dichroism (CD) is a modification of the nature of field polarization. It refers to the different absorption of left- and right-handed circularly polarized electromagnetic waves. As a consequence, linear polarization of a wave changes into an elliptical polarization while travelling through a chiral material.

Second, the existence of the Tellegen parameter leads to:

• Non-orthogonallity of the electric and magnetic field vectors

• Non-reciprocal scattering parameters

• Isotropic polarization and magnetization induced in it by either a uniform electric or magnetic field [19]

• Broken Time-reversal and space-inversion symmetries, that would lead with different phenomena, still under discussion [30]

Finally, in certain cases, the cross-coupling relation may lead to a negative refraction index [7, 8, 30], which has been confirmed experimentally in chiral media, as mentioned before [9-11]

3. Time Domain Modelling

3.1. Finite Differences in Time Domain (FDTD)

The FDTD method is a time-domain numerical technique for solving Maxwell's equations. Time-domain techniques produce wideband results with a single code execution and provide deeper physical insight into the phenomena under study. In addition, nonlinear or time-varying media are more easily handled in time-domain.

In the classical formulation of the FDTD approach the computational space is divided into cells where the electric and magnetic field components are distributed according to the Yee lattice [37]. Then, the space and time derivatives in Maxwell’s curl equations are approximated by means of second-order central differences, obtaining a system of explicit finite-difference equations that allows us compute the transient electric and the magnetic fields as function of time. The resulting algorithm is a marching-in-time procedure that simulates the actual electromagnetic waves in a finite spatial region.

The original FDTD scheme has been successfully extended to the modeling of different kinds of materials [38]. However, the magnetoelectric coupling that characterizes the bi-isotropic media makes the extension of the FDTD method to incorporate bi-isotropic media a challenging problem. First attempts can be found in [39-41]. In these works extensions of Berenger’s PML for bi-isotropic and bi-anisotropic media are developed assuming that the constitutive parameters are constant with frequency.

The schemes presented in [42-43] to incorporate bi-isotropic media are based on the field decomposition into the right- and the left-handed circularly polarized eigenwaves. These eigenwaves are uncoupled and propagate in an equivalent isotropic media with effective material parameters. Thus, the resulting sets of update equations do not differ in their general forms from the conventional FDTD update equations. Unfortunately, these schemes are only valid in some particular cases, like normal incidence on a bi-isotropic medium or free wave propagation, where the eigenwaves are not coupled. Problems like oblique incidence on a chiral slab cannot be tackled using this technique.

The formulations developed in [44-45] consider the constitutive relations proposed by Kong as governing equations. These approaches may be suitable for problems involving narrowband signals. In [44] the algorithm is formulated for the particular case of Tellegen media. In [45] the scheme is generalized for bi-isotropic media and a closed-form stability criterion, with the same physical meaning than the conventional Courant’s condition, is derived for Tellegen media.

The first formulation that considered the dispersive nature of the constitutive parameters and allowed full transient modeling of general dispersive bi-isotropic media was presented in [46-47]. In this formulation the frequency dependence of the permittivity and the permeability is represented by Lorentzian models while a single-resonance Condon model [1] [32] is adopted to represent the dispersion in the chirality parameter. In [46] the dispersion is handled using the piecewise linear recursive convolution (PLRC) technique, whereas in [47] the auxiliary differential equation (ADE) method is employed.

One important feature of the schemes given in [46-47] is the use of a new FDTD lattice developed to deal with the magnetoelectric coupling [44]. The strategy followed is to overlap two Yee cells making the x-, y-, and z-components to coincide at the same point, respectively, and, therefore, defining x-, y-, and z-nodes staggered in space and time. The advantage of this new mesh is that it simplifies the formulation since interpolation of the values of the fields is not needed. However it presents drawbacks: the standard FDTD boundary conditions are not applicable and it requires implementing a specific interface with Yee’s original lattice.

Other FDTD models for dispersive chiral media have been developed in [48-49]. These schemes use Yee’s conventional mesh and employ the Z-transform technique to convert frequency domain constitutive relations to Z-domain. However, these formulations lose the second-order accuracy in time, since in order to decouple the updating equations backward differences are used to approximate some time derivatives.

In this work we focus on an extension of the FDTD method for modeling frequency-dispersive chiral media presented in [50]. The approach consists of discretizing Maxwell’s curl equations according to Yee’s scheme and employing the Mobius transformation [51] to discretize the constitutive relations. This scheme overcomes the main limitations found in other approaches as it is a full-dispersive Model and is based on a modification of Yee’s original cell. Moreover, the scheme preserves the second-order accuracy and the explicit nature of the original FDTD formulation.

3.2. FDTD Formulation

A Differential model

In the Laplace domain, Maxwell’s curl equations can be written as

$$\begin{array}{l}s\overrightarrow{B}\text{(}\overrightarrow{r},s\text{)}=\nabla \times \overrightarrow{E}\text{(}\overrightarrow{r},s\text{)}\\ s\overrightarrow{D}\text{(}\overrightarrow{r},s\text{)}=\nabla \times \overrightarrow{H}\text{(}\overrightarrow{r},s\text{)}\end{array}$$ E5

and the constitutive relations for chiral media are

$$\begin{array}{l}\overrightarrow{D}\text{(}\overrightarrow{r},s\text{)}=\epsilon \text{(}s\text{)}\overrightarrow{E}\text{(}\overrightarrow{r},s\text{)}-\frac{1}{c}j\kappa \text{(}s\text{)}\overrightarrow{H}\text{(}\overrightarrow{r},s\text{)}\\ \overrightarrow{B}\text{(}\overrightarrow{r},s\text{)}=\mu \text{(}s\text{)}\overrightarrow{H}\text{(}\overrightarrow{r},s\text{)}+\frac{1}{c}j\kappa \text{(}s\text{)}\overrightarrow{E}\text{(}\overrightarrow{r},s\text{)}\end{array}$$ E6

where the permittivity and permeability are characterized by a Lorentz model [35] and the chirality parameter by a Condon model [32]

$$\epsilon \text{(}s\text{)}={\epsilon }_{\infty }+\frac{\text{(}{\epsilon }_s-{\epsilon }_{\infty }\text{)}{\omega }_e^2}{{\omega }_e^2+2{\delta }_{e}s+s^2}$$ E7

$$\mu \text{(}s\text{)}={\mu }_{\infty }+\frac{\text{(}{\mu }_s-{\mu }_{\infty }\text{)}{\omega }_h^2}{{\omega }_h^2+2{\delta }_{h}s+s^2}$$ E8

$$\kappa \text{'}\text{(}s\text{)}=j\kappa \text{(}s\text{)}=\frac{{\tau }_{\kappa }{\omega }_{\kappa }^{2}s}{{\omega }_{\kappa }^2+2{\delta }_{\kappa }{\omega }_{\kappa }s+s^2}$$ E9

being $${\epsilon }_{\infty }$$ ( $${\mu }_{\infty }$$ ) the permittivity (permeability) at infinite frequency, $${\epsilon }_s$$ ( $${\mu }_s$$ ) the permittivity (permeability) at zero frequency, $${\omega }_e,{\omega }_h$$ and $${\omega }_{\kappa }$$ the resonance frequencies, $${\delta }_e,{\delta }_h$$ and $${\delta }_{\kappa }$$ the loss factors and $${\tau }_{\kappa }$$ a time constant. Introducing (6) into (5) we obtain

$$\begin{array}{l}{\mu }_{\infty }s\overrightarrow{H}\text{(}\overrightarrow{r},s\text{)}=-\nabla \times \overrightarrow{E}\text{(}\overrightarrow{r},s\text{)}-{\overrightarrow{J}}_{hh}\text{(}\overrightarrow{r},s\text{)}-{\overrightarrow{J}}_{he}\text{(}\overrightarrow{r},s\text{)}\\ {\epsilon }_{\infty }s\overrightarrow{E}\text{(}\overrightarrow{r},s\text{)}=\nabla \times \overrightarrow{H}\text{(}\overrightarrow{r},s\text{)}-{\overrightarrow{J}}_{ee}\text{(}\overrightarrow{r},s\text{)}-{\overrightarrow{J}}_{eh}\text{(}\overrightarrow{r},s\text{)}\end{array}$$ E10

where the terms $${\overrightarrow{J}}_{hh},{\overrightarrow{J}}_{ee},{\overrightarrow{J}}_{eh}$$ and $${\overrightarrow{J}}_{he}$$ are auxiliary current densities defined as

$${\overrightarrow{J}}_{hh}\text{(}\overrightarrow{r},s\text{)}={\sigma }_h\text{(}s\text{)}\overrightarrow{H}\text{(}\overrightarrow{r},s\text{)}$$ E11

$${\overrightarrow{J}}_{eh}\text{(}\overrightarrow{r},s\text{)}={\varsigma }_h\text{(}s\text{)}\overrightarrow{H}\text{(}\overrightarrow{r},s\text{)}$$ E12

$${\overrightarrow{J}}_{ee}\text{(}\overrightarrow{r},s\text{)}={\sigma }_e\text{(}s\text{)}\overrightarrow{E}\text{(}\overrightarrow{r},s\text{)}$$ E13

$${\overrightarrow{J}}_{he}\text{(}\overrightarrow{r},s\text{)}={\varsigma }_e\text{(}s\text{)}\overrightarrow{E}\text{(}\overrightarrow{r},s\text{)}$$ E14

with $${\sigma }_h=s\text{(}\mu -{\mu }_{\infty }\text{), }{\varsigma }_h=-s\widehat{\kappa }/c\text{, }{\sigma }_e=s\text{(}\epsilon -{\epsilon }_{\infty }\text{)}$$ and $${\varsigma }_e=-{\varsigma }_h$$ . $${\overrightarrow{J}}_{hh},{\overrightarrow{J}}_{eh}$$ and $${\overrightarrow{J}}_{ee},{\overrightarrow{J}}_{he}$$ are evaluated at the same point and at the same time step than $$\overrightarrow{H}$$ and $$\overrightarrow{E}$$ , respectively.

B Difference model

A discrete time-domain version of (10) is obtained by applying the property $$sF\text{(}s\text{)}\leftrightarrow dF\text{(}t\text{)}/dt$$ . Then, using central differences for the derivatives, average in time for $${\overrightarrow{J}}_{hh},{\overrightarrow{J}}_{ee}$$ and average in space for $${\overrightarrow{J}}_{he},{\overrightarrow{J}}_{eh}$$ , the following expressions are obtained [50]

$$\begin{array}{l}{\overrightarrow{H}}^{n+\frac{1}{2}}\left(k+\frac{1}{2}\right)={\overrightarrow{H}}^{n-\frac{1}{2}}\left(k+\frac{1}{2}\right)-\frac{\Delta t}{{\mu }_{\infty }}\{{[\nabla \times \overrightarrow{E}]}^n\left(k+\frac{1}{2}\right)\\ +\frac{1}{2}\left[{\overrightarrow{J}}_{hh}^{n+\frac{1}{2}}\left(k+\frac{1}{2}\right)+{\overrightarrow{J}}_{hh}^{n-\frac{1}{2}}\left(k+\frac{1}{2}\right)\frac{1}{2}\right]+\frac{1}{2}\left[{\overrightarrow{J}}_{he}^n\left(k\right)+{\overrightarrow{J}}_{he}^n\left(k+1\right)\right]\}\end{array}$$ E15

$$\begin{array}{l}{\overrightarrow{E}}^{n+1}\left(k\right)={\overrightarrow{E}}^n\left(k\right)+\frac{\Delta t}{{\epsilon }_{\infty }}\{{[\nabla \times \overrightarrow{H}]}^{n+\frac{1}{2}}\left(k\right)\\ -\frac{1}{2}\left[{\overrightarrow{J}}_{ee}^{n+1}\left(k\right)+{\overrightarrow{J}}_{ee}^n\left(k\right)\right]-\frac{1}{2}\left[{\overrightarrow{J}}_{eh}^{n+\frac{1}{2}}\left(k-\frac{1}{2}\right)+{\overrightarrow{J}}_{eh}^{n+\frac{1}{2}}\left(k+\frac{1}{2}\right)\right]\}\end{array}$$ E16

Figure 1.

Transformation of the constitutive equations into the discrete-time domain.

We should now derive discrete-time expressions for the constitutive relations (11)-(14). The procedure is outlined in Figure 1. One attempt may consist on first transforming them into the continuous-time domain by again using the property $$sF\text{(}s\text{)}\leftrightarrow dF\text{(}t\text{)}/dt$$ , and then approximating the resulting second-order differential equations by using finite differences. However, there are different ways to perform this discretization and each of them leads to different algorithms with different numerical properties. In this study, the constitutive relations (11)-(14) are discretized following an alternative way, first transforming them into the Z-domain and then obtaining the discrete-time domain version. Thus, for (11) we obtain

$${\overrightarrow{J}}_{hh}\text{(}\overrightarrow{r},Z\text{)}={\sigma }_h\text{(}Z\text{)}\overrightarrow{H}\text{(}\overrightarrow{r},Z\text{)}$$

where $${\sigma }_h\text{(}Z\text{)}$$ is the Z-domain magnetic conductivity. To obtain $${\sigma }_h\text{(}Z\text{)}$$ we adopt Mobius transformation technique [51]. This approach involves the following change of variable

$$s=\frac{2}{\Delta t}\frac{1-Z^{-1}}{1+Z^{-1}}$$

Applying this change to $${\sigma }_h\text{(}Z\text{)}$$ yields

$${\overrightarrow{J}}_{hh}\text{(}\overrightarrow{r},Z\text{)}=\frac{c_0^{({\sigma }_h)}+c_1^{({\sigma }_h)}{Z}^{-1}+c_2^{({\sigma }_h)}{Z}^{-2}}{1+d_1^{({\sigma }_h)}{Z}^{-1}+d_2^{({\sigma }_h)}{Z}^{-2}}\overrightarrow{H}\text{(}\overrightarrow{r},Z\text{)}$$ E17

where $$c_m^{({\sigma }_h)}$$ and $$d_m^{({\sigma }_h)}$$ are real-valued coefficients. Now considering the Z-transform property $$Z^{-m}F\text{(}Z\text{)}\leftrightarrow F^{n-m}$$ equation (17) can be written in difference form, as follows

$${\overrightarrow{J}}_{hh}^{n+\frac{1}{2}}={\displaystyle \sum _{m=0}^2{c}_m^{({\sigma }_h)}{\overrightarrow{H}}_{}^{n+\frac{1}{2}-m}}-{\displaystyle \sum _{m=1}^2{d}_m^{({\sigma }_h)}{\overrightarrow{J}}_{hh}^{n+\frac{1}{2}-m}}$$

To reduce the memory requirements this expression is interpreted as an infinite-impulse response digital filter and is implemented by using the transposed direct form II. This type of implementation is a canonical form, which guarantees the minimum number of delay elements, thus the minimum number of additional storage variables. Applying this approach, the resulting difference equations are

$${\overrightarrow{J}}_{hh}^{n+\frac{1}{2}}={\overrightarrow{W}}_{hh,1}^{n-\frac{1}{2}}+c_0^{({\sigma }_h)}{\overrightarrow{H}}^{n+\frac{1}{2}}$$ E18

$${\overrightarrow{W}}_{hh,1}^{n+\frac{1}{2}}={\overrightarrow{W}}_{hh,2}^{n-\frac{1}{2}}-d_1^{({\sigma }_h)}{\overrightarrow{J}}_{hh}^{n+\frac{1}{2}}+c_1^{({\sigma }_h)}{\overrightarrow{H}}^{n+\frac{1}{2}}$$ E19

$${\overrightarrow{W}}_{hh,2}^{n+\frac{1}{2}}=c_2^{({\sigma }_h)}{\overrightarrow{H}}^{n+\frac{1}{2}}-d_2^{({\sigma }_h)}{\overrightarrow{J}}_{hh}^{n+\frac{1}{2}}$$ E20

These equations are evaluated at $$z=(k+(1/2))\Delta z$$ . Note, that they are coupled to (15). Decoupling them an eliminating $${\overrightarrow{W}}_{hh,2}$$ we obtain

$$\begin{array}{l}{\overrightarrow{H}}^{n+\frac{1}{2}}\left(k+\frac{1}{2}\right)=\frac{1}{2{\mu }_{\infty }+\Delta t{c}_0^{({\sigma }_h)}}(2{\mu }_{\infty }{\overrightarrow{H}}^{n-\frac{1}{2}}\left(k+\frac{1}{2}\right)\\ -\Delta t\left\{2{[\nabla \times \overrightarrow{E}]}^n\left(k+\frac{1}{2}\right)+{\overrightarrow{W}}_{hh}^{n-\frac{1}{2}}\left(k+\frac{1}{2}\right)+{\overrightarrow{J}}_{hh}^{n-\frac{1}{2}}\left(k+\frac{1}{2}\right)+{\overrightarrow{J}}_{he}^n\left(k\right)+{\overrightarrow{J}}_{he}^n\left(k+1\right)\right\})\end{array}$$ E21

and

$$\begin{array}{l}{\overrightarrow{J}}_{hh}^{n+\frac{1}{2}}={\overrightarrow{W}}_{hh}^{n-\frac{1}{2}}+c_0^{({\sigma }_h)}{\overrightarrow{H}}^{n+\frac{1}{2}}\\ {\overrightarrow{W}}_{hh}^{n+\frac{1}{2}}=c_2^{({\sigma }_h)}{\overrightarrow{H}}^{n-\frac{1}{2}}-d_2^{({\sigma }_h)}{\overrightarrow{J}}_{hh}^{n-\frac{1}{2}}-d_1^{({\sigma }_h)}{\overrightarrow{J}}_{hh}^{n+\frac{1}{2}}+c_1^{({\sigma }_h)}{\overrightarrow{H}}^{n+\frac{1}{2}}\end{array}$$ E22

where $${\overrightarrow{W}}_{hh}\equiv {\overrightarrow{W}}_{hh,1}$$ . Using the same procedure to discretize (12), we obtain the following equation, which is evaluated at $$z=(k+(1/2))\Delta z$$

$$\begin{array}{l}{\overrightarrow{J}}_{eh}^{n+\frac{1}{2}}={\overrightarrow{W}}_{eh}^{n-\frac{1}{2}}+c_0^{({\xi }_h)}{\overrightarrow{H}}^{n+\frac{1}{2}}\\ {\overrightarrow{W}}_{eh}^{n+\frac{1}{2}}=c_2^{({\xi }_h)}{\overrightarrow{H}}^{n-\frac{1}{2}}-d_2^{({\xi }_h)}{\overrightarrow{J}}_{eh}^{n-\frac{1}{2}}-d_1^{({\xi }_h)}{\overrightarrow{J}}_{eh}^{n+\frac{1}{2}}+c_1^{({\xi }_h)}{\overrightarrow{H}}^{n+\frac{1}{2}}\end{array}$$ E23

After that (13) is discretized in the same way than (11) y (12). The resulting equations are coupled to (16). Decoupling them we get to the next expression

$$\begin{array}{l}{\overrightarrow{E}}^{n+1}\left(k\right)=\frac{1}{2{\epsilon }_{\infty }+\Delta t{c}_0^{({\sigma }_e)}}(2{\epsilon }_{\infty }{\overrightarrow{E}}^n\left(k\right)\\ +\Delta t\left\{2{[\nabla \times \overrightarrow{H}]}^{n+\frac{1}{2}}\left(k\right)-{\overrightarrow{W}}_{ee}^n\left(k\right)-{\overrightarrow{J}}_{ee}^n\left(k\right)-{\overrightarrow{J}}_{eh}^{n+\frac{1}{2}}\left(k-\frac{1}{2}\right)-{\overrightarrow{J}}_{eh}^{n+\frac{1}{2}}\left(k+\frac{1}{2}\right)\right\})\end{array}$$ E24

and

$$\begin{array}{l}{\overrightarrow{J}}_{ee}^{n+1}={\overrightarrow{W}}_{ee}^n+c_0^{({\sigma }_e)}{\overrightarrow{E}}^{n+1}\\ {\overrightarrow{W}}_{ee}^{n+1}=c_2^{({\sigma }_e)}{\overrightarrow{E}}^n-d_2^{({\sigma }_e)}{\overrightarrow{J}}_{ee}^n-d_1^{({\sigma }_e)}{\overrightarrow{J}}_{ee}^{n+1}+c_1^{({\sigma }_e)}{\overrightarrow{E}}^{n+1}\end{array}$$ E25

Finally, (14) is discretized using the Mobius transformation, as was previously done for (11)–(13), and implemented as

$$\begin{array}{l}{\overrightarrow{J}}_{he}^{n+1}={\overrightarrow{W}}_{he}^n+c_0^{({\xi }_e)}{\overrightarrow{E}}^{n+1}\\ {\overrightarrow{W}}_{he}^{n+1}=c_2^{({\xi }_e)}{\overrightarrow{E}}^n-d_2^{({\xi }_e)}{\overrightarrow{J}}_{he}^n-d_1^{({\xi }_e)}{\overrightarrow{J}}_{he}^{n+1}+c_1^{({\xi }_e)}{\overrightarrow{E}}^{n+1}\end{array}$$ E26

Equations (25) and (26) are evaluated at $$z=k\Delta z$$ .

Thus, the resulting computational procedure comprises the following calculations on each time step:

1. $${\overrightarrow{H}}^{n+(1/2)}$$ is calculated by means of (21)

2. $${\overrightarrow{J}}_{hh}^{n+(1/2)}$$ and $${\overrightarrow{W}}_{hh}^{n+(1/2)}$$ are updated using (22), where $${\overrightarrow{H}}^{n+(1/2)}$$ is known from step 1.

3. $${\overrightarrow{J}}_{eh}^{n+(1/2)}$$ and $${\overrightarrow{W}}_{eh}^{n+(1/2)}$$ are obtained by using (23), where $${\overrightarrow{H}}^{n+(1/2)}$$ has been obtain in 1.

4. $${\overrightarrow{E}}^{n+1}$$ is evaluated by using (24)

5. $${\overrightarrow{J}}_{ee}^{n+1}$$ and $${\overrightarrow{W}}_{ee}^{n+1}$$ are calculated by means of (25), where $${\overrightarrow{E}}^{n+1}$$ is known from step 4.

6. $${\overrightarrow{J}}_{he}^{n+1}$$ and $${\overrightarrow{W}}_{he}^{n+1}$$ are updated by using (26), where $${\overrightarrow{E}}^{n+1}$$ has been obtained in 4.

This algorithm preserves the second-order accuracy and the explicit nature of the conventional FDTD formulation with only 4 additional back-stored variables per field component and cell.

3.3. Results

Wave propagation in a chiral slab

FDTD schemes emulate the progression of the fields as they actually evolve in space and time. This feature of time-domain simulators will allow us to visualize the noteworthy properties of wave propagation in chiral media.

Figure 2.

Evolution of the polarization of a wave traveling in a chiral slab: (a) lossless case: optical activity, (b) lossy case: optical activity and circular dichroism.

Thus, in a first simulation we have considered a time-harmonic electric field that propagates in the z positive direction and impinges on a slab of chiral material. Figure 2.a shows a snapshot of the electric field in the entire simulation domain. It can be seen how, due to the optical activity, the plane of polarization rotates as the wave travels through the chiral slab.

When losses are present in a chiral medium circular dichroism appears. The circular dichroism modifies the nature of the field polarization. To visualize this phenomenon we performed the same simulation than in the previous case but now considering a lossy chiral slab. Figure 2.b shows how the polarization rotates as the wave propagates through the chiral slab and how the linearly polarized wave degenerates progressively into an elliptically polarized wave.

Reflection and transmission from a chiral slab

Consider the propagation of an electromagnetic wave that travels in air and impinges on a slab of chiral medium. In a chiral medium, the reflected and transmitted waves can be obtained through the following reflection and transmission dyadics

$$\overline{\overline{R}}=\left(\begin{array}{cc}{R}_{co}& 0\\ 0& R_{co}\end{array}\right)\overline{\overline{T}}=\left(\begin{array}{cc}{T}_{co}& T_{cr}\\ T_{cr}& T_{co}\end{array}\right)$$

being $$R_{co}$$ the copolarized reflection coefficient and $$T_{co},T_{cr}$$ the co- and crosspolarized transmission coefficients, respectively. The crosspolarized reflection coefficient equals zero.

The dispersive behavior of the constitutive parameters were assumed to follow the models given in (7)-(9) with the following values $${\epsilon }_s=4.2{\epsilon }_0,$$ $${\epsilon }_{\infty }=3.3{\epsilon }_0,$$ $${\mu }_s=1.1{\mu }_0,$$ $${\mu }_{\infty }={\mu }_0,$$ $${\omega }_e={\omega }_h={\omega }_{\kappa }=2\pi \times 7.5GHz$$ , $${\delta }_e=0.05{\omega }_e,$$ $${\delta }_h=0.05{\omega }_h,$$ $${\delta }_{\kappa }=0.07$$ and $${\tau }_{\kappa }=1ps.$$

Figure 3.

Reflection and transmission coefficients: theoretical values and relative errors.

A band limited Gaussian pulse was injected at a source point in the air region. The time step was $$\Delta t=1\text{ps}$$ and the space discretization $$\Delta z=0.3\text{mm}$$ . The thickness of the slab was $$d=30\Delta z.$$ Figure 3 depicts the theoretical reflection and transmission coefficients [1] and the relative error of the coefficients computed by the present FDTD formulation. These results demonstrate the accuracy and the full transient capability of this FDTD approach.

4. Multiresolution in Time Domain (MRTD)

Multiresolution in Time Domain (MRTD) technique is a new approach to reduce the simulation time of FDTD schemes (actually, it may be considered a generalization of FDTD [52]), keeping the same degree of accuracy. In 1996, Krumpholz showed that Yee’s FDTD scheme could be derived by applying the method of moments for the discretization of Maxwell’s equations, using pulse basis functions for the expansion of the unknown fields [53]. In a MRTD scheme the expansion is completed by means of a twofold expansion in scaling and wavelet functions with respect to time/space. Scaling functions takes into account smoothly varying fields. In regions or time periods characterized by strong field variations or field singularities, higher resolution is enhanced by incorporating wavelet functions in the field expansions. The major advantage of the use of MRTD in the time domain is the capability to develop time and space adaptive grids. The scheme is, then, a generalization of Yee’s FDTD and can extend this technique’s capabilities by improving computational efficiency and substantially reducing computer resources.

The MRTD approach has been successfully applied to finite-difference time-domain (FDTD) simulations of passive structures [53-56]. In those cases, multiresolution expansion has been usually applied to the space dependence of the electromagnetic fields. In [57], the expansion was applied to the time dependence of the fields. In this case, such expansion has the purpose of improving the resolution in the calculation of convolutions in dispersive media. As we have seen previously, the chirality parameter, present in eq. (3) is intrinsically dispersive (eqs. (7)-(9)), so we expand here a similar technique suitable for chiral media.

Figure 4.

Haar wavelet functions: left, scale function Φ _n (t). Right, first order wavelet function ψ _n (t). Both cases are particularized for n = 1.

To discretize the constitutive equations (2), we must develop all the components of the field vectors ( $\overrightarrow{E},\overrightarrow{D},\overrightarrow{B},\overrightarrow{H}$ ) in scale and wavelet functions:

$$F={\displaystyle \sum _i{\varphi }_i\left(z\right){\displaystyle \sum _n\left({}_{n}F{}_i^{\varphi }{\varphi }_n\left(t\right)+{}_{n}F{}_i^{\psi }{\psi }_n\left(t\right)\right)}}$$ E27

F may be any component of any of the field vectors above mentioned. We have included two levels of resolution in the time dependence of the function, while the space discretization is represented only in one dimension, assuming dependence on z. The extension to three dimensions is straightforward: if the scale (Φ) and wavelet (ψ) functions are Haar functions [58], as shown in Figure 4, the space discretization is equivalent to the FDTD discretization explained in the previous section. That means that the discretization of Maxwell curl equations will lead to two parallel Yee’s networks, one for scale coefficients ( $${}_{n}F{}_i^{\varphi }$$ ), and a second one for wavelet coefficients ( $${}_{n}F{}_i^{\psi }$$ ).

The results for both networks are combined in the constitutive relations (Eq. (3)). If we discretize them in time domain, using a recursive convolution method (in a similar way to FDTD [46]), we obtain the following expressions (x component of $\overrightarrow{D},\overrightarrow{B}$ ):

$$\begin{array}{c}{D}^x\left(z,t\right)={\displaystyle \sum _i\left\{\begin{array}{c}{\displaystyle \sum _n\epsilon \left({}_{n}E{}_i^{x,\varphi }{\varphi }_n\left(t\right)+{}_{n}E{}_i^{x,\psi }{\psi }_n\left(t\right)\right)}\\ -{\displaystyle \sum _m\left\{\frac{1}{c}{\displaystyle {\int }_0^t\left({}_{m\text{'}}H{}_i^{x,\varphi }{\varphi }_m\left(t-\tau \right)+{}_{m}H{}_i^{x,\psi }{\psi }_m\left(t-\tau \right)\right)\kappa \text{'}\left(\tau \right)d\tau }\right\}}\end{array}\right\}{\varphi }_i\left(z\right)}\\ B^x\left(z,t\right)={\displaystyle \sum _i\left\{\begin{array}{c}{\displaystyle \sum _n\mu \left({}_{n}H{}_i^{x,\varphi }{\varphi }_n\left(t\right)+{}_{n}H{}_i^{x,\psi }{\psi }_n\left(t\right)\right)}\\ +{\displaystyle \sum _m\left\{\frac{1}{c}{\displaystyle {\int }_0^t\left({}_{m}E{}_i^{x,\varphi }{\varphi }_m\left(t-\tau \right)+{}_{m}E{}_i^{x,\psi }{\psi }_m\left(t-\tau \right)\right)\kappa \text{'}\left(\tau \right)d\tau }\right\}}\end{array}\right\}{\varphi }_i\left(z\right)}\end{array}$$ E28

where κ’(τ) is the inverse Laplace transform of the chirality parameter in frequency domain κ(ω). In this equation $${}_{n}E{}_i^{x,\varphi }$$ , for example, would represent the coefficient for the scale function of the x-component of the electric field $$\overrightarrow{E}$$ , when t = n t and z = i Δ z. Note that the equation includes only the dispersion in the chirality response: dielectric and magnetic response could be included as in [57]. If we are working with Haar functions, this equation may be simplified [59], obtaining:

$$\begin{array}{c}{D}^x\left(z,t\right)={\displaystyle \sum _{i,n}\left\{\begin{array}{c}\epsilon \left({}_{n}E{}_i^{x,\varphi }{\varphi }_n\left(t\right)+{}_{n}E{}_i^{x,\psi }{\psi }_n\left(t\right)\right)\\ -\frac{1}{c}{\varphi }_n\left(t\right){\displaystyle \sum _m\left\{{}_{n-m}H{}_i^{x,\varphi }\kappa {\text{'}}_m^{\varphi }+{}_{n-m}H{}_i^{x,\psi }\kappa {\text{'}}_m^{\psi }\right\}}\end{array}\right\}{\varphi }_i\left(z\right)}\\ B^x\left(z,t\right)={\displaystyle \sum _{i,n}\left\{\begin{array}{c}\mu \left({}_{n}H{}_i^{x,\varphi }{\varphi }_n\left(t\right)+{}_{n}H{}_i^{x,\psi }{\psi }_n\left(t\right)\right)\\ +\frac{1}{c}{\varphi }_n\left(t\right){\displaystyle \sum _m\left\{{}_{n-m}E{}_i^{x,\varphi }\kappa {\text{'}}_m^{\varphi }+{}_{n-m}E{}_i^{x,\psi }\kappa {\text{'}}_m^{\psi }\right\}}\end{array}\right\}}{\varphi }_i\left(z\right)\end{array}$$ E29

where $$\kappa {\text{'}}_{n-m}^{\varphi }={\displaystyle {\int }_{\left(n-m\right)\Delta t}^{\left(n-m+1\right)\Delta t}\kappa \text{'}\left(\tau \right)d\tau }$$ , $$\kappa {\text{'}}_{n-m}^{\psi }={\displaystyle {\int }_{\left(n-m+1/2\right)\Delta t}^{\left(n-m+1\right)\Delta t}\kappa \text{'}\left(\tau \right)d\tau }-{\displaystyle {\int }_{\left(n-m\right)\Delta t}^{\left(n-m+1/2\right)\Delta t}\kappa \text{'}\left(\tau \right)d\tau }$$ .If we compare these equations with (27), in order to separate scale and wavelet coefficients of D ^x and B ^x , we find:

$$\begin{array}{c}{}_{n}D{}_i^{x,\varphi }=\epsilon {}_{n}E{}_i^{x,\varphi }-\frac{1}{c_0}{\displaystyle \sum _m\left\{{}_{n-m}E{}_i^{x,\varphi }\kappa {\text{'}}_m^{\varphi }+{}_{n-m}E{}_i^{x,\psi }\kappa {\text{'}}_m^{\psi }\right\}}\\ {}_{n}B{}_i^{x,\varphi }=\mu {}_{n}H{}_i^{x,\varphi }+\frac{1}{c_0}{\displaystyle \sum _m\left\{{}_{n-m}E{}_i^{x,\varphi }\kappa {\text{'}}_m^{\varphi }+{}_{n-m}E{}_i^{x,\psi }\kappa {\text{'}}_m^{\psi }\right\}}\end{array}$$ E30

$$\begin{array}{c}{}_{n}D{}_i^{x,\psi }=\epsilon \cdot {}_{n}E{}_i^{x,\psi }\\ {}_{n}B{}_i^{x,\psi }=\mu \cdot {}_{n}H{}_i^{x,\psi }\end{array}$$ E31

As we may see, the convolutions sums appear in eq. (30), i.e., they affect only the coefficients of the scale functions of $$\overrightarrow{D}$$ and $\overrightarrow{B}$ ( $${}_{n}D{}_i^{x,\varphi }$$ and $${}_{n}B{}_i^{x,\varphi }$$ ): if the wavelet functions give account of the faster variations of the fields, it is logic that they do not depend on the long term values stored in the integral is reasonable.

So, we may follow a procedure similar to the one described in [57] for isotropic dispersive materials, i.e. the algorithm would be composed by two different steps:

1. Propagation of the field, following the discretized Maxwell curl equations. We obtain updated values of vectors $\overrightarrow{D}$ and $\overrightarrow{B}$ , using the previous values of $\overrightarrow{E}$ and $\overrightarrow{H}$ .

2. Updating of the values of $\overrightarrow{E}$ and $\overrightarrow{H}$ , using the constitutive relations of the medium; we should update first the wavelet (ψ) coefficients, following Eq. (31), because their value is necessary to solve Eq. (30) and obtain the new values for the scale (ϕ) coefficients.

4.1. Results

To validate numerically the accuracy of our approximation, we have simulated the propagation of a plane wave linearly polarized through a chiral medium, with ε =2ε₀, μ=μ₀, and κ following the Condon model: ω ₀=1 GHz, τ=3 ps and a damping factor δ=0 (medium without losses). To simplify the results, we have not taken into account the dielectric and magnetic dispersions (they do not affect the electromagnetic activity). The discretization is Δ t = 1 ps and Δ x = 0.375 mm. We have introduced a harmonic source, with a frequency of 7 GHz, and after 8000 time steps, we have calculated the rotation of the polarization plane in different points: the result is compared with the one obtained using a FDTD scheme [46]. The results are shown in Table 1: [59]

	Theoretical value	FDTD		MRTD
Value	Error	Value	Error
1000Δz	11.55º	11.54º	0.09%	11.54º	0.09%
2000Δz	23.10º	23.09º	0.04%	23.08º	0.09%
3000Δz	34.64º	34.56º	0.2%	34.65º	0.03%
4000Δz	46.21º	46.11º	0.2%	46.04º	0.15%

Table 1.

Rotation of the polarization plane during the propagation of a plane wave through a chiral medium [59]

5. Transmission Line Matrix (TLM)

The Transmission Line Matrix (TLM) method is a numerical technique that solves the wave equation of a propagation phenomenon in time domain. As in the previous cases, the basic TLM network [60] must be modified to include the cross terms emerging in the Maxwell equations, thus in the wave equation, when dealing with bi-isotropic media. We consider the propagation of waves in the xy–plane: although it is a bi-dimensional problem, the medium presents cross coupling on all the field components on every plane perpendicular to propagation direction, so the number of unknown field components cannot be reduced less than twelve.

After a simple manipulation, Maxwell equations can be arranged in four sets of equations as shown in the Table 2, related by the constitutive relations described in equation (4).

where Φ is defined as $${\Phi }^2=\epsilon \mu -{\chi }^2{\epsilon }_0{\mu }_0$$ , and must have a non–zero value. This mathematical condition, verified by all bi-isotropic media [1] means that any wave solution has associated a finite value of the phase velocity.

$$\left\{1\right\}$$	$$\frac{\partial E_z}{\partial x}=\frac{\partial B_y}{\partial t},$$	$$\frac{\partial E_z}{\partial y}=-\frac{\partial B_x}{\partial t},$$
$$\left\{2\right\}$$	$$\frac{\partial H_z}{\partial x}=-\frac{\partial D_y}{\partial t},$$
$$\left\{3\right\}$$	$$\frac{\partial D_z}{\partial x}={\Phi }^2\frac{\partial H_y}{\partial t}+\frac{\kappa }{c}\frac{{\partial }^2{Q}_{ey}}{\partial t^2},$$	$$\frac{\partial D_z}{\partial y}=-{\Phi }^2\frac{\partial H_x}{\partial t}-\frac{\kappa }{c}\frac{{\partial }^2{Q}_{ex}}{\partial t^2},$$	$$\frac{\partial H_y}{\partial x}-\frac{\partial H_x}{\partial y}=\frac{\partial D_z}{\partial t}$$
$$\left\{4\right\}$$	$$\frac{\partial B_z}{\partial x}=-{\Phi }^2\frac{\partial E_y}{\partial t}-\frac{\kappa }{c}\frac{{\partial }^2{Q}_{my}}{\partial t^2},$$	$$\frac{\partial B_z}{\partial y}={\Phi }^2\frac{\partial E_x}{\partial t}+\frac{\kappa }{c}\frac{{\partial }^2{Q}_{mx}}{\partial t^2},$$	$$\frac{\partial E_y}{\partial x}-\frac{\partial E_x}{\partial y}=-\frac{\partial B_z}{\partial t}$$

Table 2.

The Four Sets of Electromagnetic Field Vectors Equations [62]

To obtain a simplified representation, two next vectors, with no physical meaning, have been introduced:

$$\begin{array}{cc}{\overrightarrow{Q}}_m=\mu \overrightarrow{H}+\frac{\chi }{c}\overrightarrow{E}+\overrightarrow{B},& \begin{array}{cc}& \end{array}{\overrightarrow{Q}}_e=\epsilon \overrightarrow{E}+\frac{\chi }{c}\overrightarrow{H}+\overrightarrow{D}\end{array}$$ E32

For every set of resulting equations, we may derive an independent wave equation describing the propagation, and having the form:

$$\frac{{\partial }^{2}A}{\partial x^2}+\frac{{\partial }^{2}A}{\partial y^2}-{\varphi }^2\frac{{\partial }^{2}A}{\partial t^2}+\frac{\kappa }{c}\frac{{\partial }^3}{\partial t^3}\left[R+\frac{\partial S}{\partial t}\right]=0$$ E33

A	R	S
$$E_z$$	$$Q_{mz}$$	0
$$H_z$$	$$-Q_{ez}$$	0
$$D_z$$	$$\epsilon Q_{mz}-\frac{\chi }{c}{Q}_{ez}$$	$$Q_{ez}$$
$$B_z$$	$$\frac{\chi }{c}{Q}_{mz}-\mu Q_{ez}$$	$$Q_{mz}$$

Table 3.

Symbols $$A$$ , $$R$$ and $$S$$ in wave equations [62]

The symbol A represents whatever z-component of the electromagnetic field vectors, and each symbol R or S equals the z-component of a linear combination of $${\overrightarrow{Q}}_m$$ and $${\overrightarrow{Q}}_e$$ (both introduced in (32)) depending on the component A, as detailed in Table 3.

5.1. TLM Modeling

The classical TLM method works with interconnected transmission–lines along the coordinate axes establishing a network. The points at which the transmission–lines intersect are referred to as “nodes”. Then, space and time are discretized, and voltage and current pulses scatter from point to point in space, $$\Delta \ell$$ , in a fixed time step $$\Delta t$$ , [60]. The starting point to obtain any TLM model is to discretize the electromagnetic field equations (either the Maxwell equations in Table 2 or the wave equations (33)), and to compare it with the equations modeled by the TLM algorithm, that is, the equations relating the voltages and currents in the transmission–line network. For a parallel connection of the transmission–line in the xy–plane, the voltages and currents are related by the differential equations:

$$\begin{array}{ccccc}\frac{\partial V_z}{\partial x}=-\text{L}\frac{\partial I_x}{\partial t}& & \frac{\partial V_z}{\partial y}=-\text{L}\frac{\partial I_y}{\partial t}& & \frac{\partial I_x}{\partial x}+\frac{\partial I_y}{\partial y}=-2\text{C}\frac{\partial V_z}{\partial t}\end{array}$$ E34

where L and C represent the inductance and the capacitance per unit length of the transmission–lines, respectively. In these equations, the voltage V _z , and the currents I _x and I _y symbolize the standard electric quantities in the transmission–line network [60].

A look at Table 2 allows us to implement the model of every set of equations with a classical TLM mesh of shunt-connected nodes. Furthermore, the twelve field vectors components are required in the new model due to the existing magneto–electric coupling, even in the 2D case. In this way, a single set is not enough, but four meshes are needed to simulate each of the twelve field components. Every TLM network will be labeled {1}, {2}, {3}, or {4}, following their equivalent set of Maxwell equations in Table 2. This leads to a specific analogy between each TLM network electric quantity and its respective electromagnetic field vectors component, as shown in Table 4 where ϕ r = ϕ / c .

	$$\left\{1\right\}$$	$$\left\{2\right\}$$	$$\left\{3\right\}$$	$$\left\{4\right\}$$
$$V_z\equiv \mathrm{...}$$	$E_z$	$$-{\mu }_{0}c{H}_z$$	$$D_z/{\epsilon }_0$$	$$-c{B}_z$$
$$-I_x\equiv \mathrm{...}$$	$$B_y/{\mu }_0$$	$$c{D}_y$$	$${\Phi }_r^2{H}_y$$	$${\Phi }_r^2{\epsilon }_{0}c{E}_y$$
$$I_y\equiv \mathrm{...}$$	$$B_x/{\mu }_0$$	$$c{D}_x$$	$${\Phi }_r^2{H}_x$$	$${\Phi }_r^2{\epsilon }_{0}c{E}_x$$

Table 4.

Analogies between TLM meshes quantities and field vectors [62].

$$\left\{1\right\}$$	$$\frac{\partial V_z}{\partial x}=-{\mu }_0\frac{\partial I_x}{\partial t},$$	$$\frac{\partial V_z}{\partial y}=-{\mu }_0\frac{\partial I_y}{\partial t},$$	$$\frac{\partial I_x}{\partial x}+\frac{\partial I_y}{\partial y}=-\frac{{\phi }^2}{{\mu }_0}\frac{\partial V_z}{\partial t}+\frac{\kappa }{{\eta }_0}\frac{{\partial }^2{Q}_{mz}}{\partial t^2}$$
$$\left\{2\right\}$$	$$\frac{\partial V_z}{\partial x}=-{\mu }_0\frac{\partial I_x}{\partial t},$$	$$\frac{\partial V_z}{\partial y}=-{\mu }_0\frac{\partial I_y}{\partial t},$$	$$\frac{\partial I_x}{\partial x}+\frac{\partial I_y}{\partial y}=-\frac{{\phi }^2}{{\mu }_0}\frac{\partial V_z}{\partial t}+\kappa \frac{{\partial }^2{Q}_{ez}}{\partial t^2}$$
$$\left\{3\right\}$$	$$\frac{\partial V_z}{\partial x}=-{\mu }_0\frac{\partial I_x}{\partial t}+\kappa {\eta }_0\frac{{\partial }^2{Q}_{ey}}{\partial t^2},$$	$$\frac{\partial V_z}{\partial y}=-{\mu }_0\frac{\partial I_y}{\partial t}-\kappa {\eta }_0\frac{{\partial }^2{Q}_{ex}}{\partial t^2},$$	$$\frac{\partial I_x}{\partial x}+\frac{\partial I_y}{\partial y}=-{\phi }_r^2{\epsilon }_0\frac{\partial V_z}{\partial t}$$
$$\left\{4\right\}$$	$$\frac{\partial V_z}{\partial x}=-{\mu }_0\frac{\partial I_x}{\partial t}-\kappa \frac{{\partial }^2{Q}_{my}}{\partial t^2},$$	$$\frac{\partial V_z}{\partial y}=-{\mu }_0\frac{\partial I_y}{\partial t}+\kappa \frac{{\partial }^2{Q}_{mx}}{\partial t^2},$$	$$\frac{\partial I_x}{\partial x}+\frac{\partial I_y}{\partial y}=-{\phi }_r^2{\epsilon }_0\frac{\partial V_z}{\partial t}$$

Table 5.

The differential equations for the four TLM meshes [62].

The corresponding TLM equivalences of the vectors $${\overrightarrow{Q}}_m$$ and $${\overrightarrow{Q}}_e$$ , introduced in (32), are:

$$\begin{array}{ccc}\begin{array}{l}{Q}_{mx}=\frac{{\mu }_0}{{\varphi }_r^2}\left[{\mu }_r{I}_y\left\{3\right\}+\chi I_y\left\{4\right\}+{\varphi }_r^2{I}_y\left\{1\right\}\right]\\ Q_{my}=-\frac{{\mu }_0}{{\varphi }_r^2}\left[{\mu }_r{I}_x\left\{3\right\}+\chi I_x\left\{4\right\}+{\varphi }_r^2{I}_x\left\{1\right\}\right]\\ Q_{mz}=-\frac{1}{c}\left[{\mu }_r{V}_z\left\{2\right\}+\chi V_z\left\{1\right\}+V_z\left\{4\right\}\right]\end{array}\}& & \begin{array}{l}{Q}_{ex}=\frac{1}{{\varphi }_r^{2}c}\left[{\epsilon }_r{I}_y\left\{4\right\}+\chi I_y\left\{3\right\}+{\varphi }_r^2{I}_y\left\{2\right\}\right]\\ Q_{ey}=-\frac{1}{{\varphi }_r^{2}c}\left[{\epsilon }_r{I}_x\left\{4\right\}+\chi I_x\left\{3\right\}+{\varphi }_r^2{I}_x\left\{2\right\}\right]\\ Q_{ez}=\frac{}{}{\epsilon }_0\left[{\epsilon }_r{V}_z\left\{1\right\}-\chi V_z\left\{2\right\}+V_z\left\{3\right\}\right]\end{array}\}\end{array}$$ E35

Again, these vectors have no special physical meaning. Please note that whatever the relationship between field components and TLM variables is, these networks have the same wave propagation characteristics of the classic TLM network, well studied in the literature, [60]. Therefore, the numerical dispersion and consequently the velocity error of the bi-isotropic TLM network are well known.

To take into account the propagation velocity (the same in the four meshes following (3)), a permittivity stub is shunt connected at all the nodes. The normalized admittance of these stubs is $$y_o=$$ 4( $${\varphi }_{\text{r}}^{\text{2}}-1$$ ). Then, propagation velocities into each TLM network depend on $$\varphi$$ , different from the velocity in the modeled bi-isotropic medium. The same transmission–lines characteristics are imposed for all the meshes, and consequently the TLM meshes impedances $$Z_{TLM}$$ equal $${\eta }_o/{\varphi }_r$$ , where $${\eta }_o$$ denotes the vacuum intrinsic impedance.

Additionally, a new circuit element is needed at every node in order to incorporate the terms in Table 2 with no equivalence in Table 4. The technique developed in [60] to model the waves propagation in complex media, based on the voltage source connection, is used for this purpose. It is important to remark that, since these voltages sources $$V_s^k$$ will introduce the adequate signal correction at any time iteration, we obtain a TLM algorithm that models the correct propagation velocity.

The simulation begins with the excitation of one or several field vectors components, depending on the problem (symmetry, boundaries, output, …). Regarding the equivalences in Table 4, voltage impulses are initially introduced into the system according to the desired field excitation. The scattering and incidence processes of these voltage pulses follow the classical TLM method procedures [60], but, for a bi-isotropic 2D–TLM model, the algorithm is not yet completed, since at every time iteration k, the voltages at each node must be transferred from a pair of networks to the other pair:

$$\begin{array}{c}{V}_n\left\{3\right\}={\epsilon }_r{V}_n\left\{1\right\}-\chi V_n\left\{2\right\}+\frac{\kappa }{\Delta t}{\Delta }_t^1{V}_n\left\{2\right\}\\ V_n\left\{4\right\}={\mu }_r{V}_n\left\{2\right\}-\chi V_n\left\{1\right\}-\frac{\kappa }{\Delta t}{\Delta }_t^1{V}_n\left\{1\right\}\\ V_n\left\{1\right\}=\frac{1}{{\varphi }_r^2}\left[{\mu }_r{V}_n\left\{3\right\}+\chi V_n\left\{4\right\}-\frac{\kappa }{\Delta t}{\Delta }_t^1\left({\mu }_r{V}_n\left\{2\right\}-\chi V_n\left\{1\right\}\right)\right]\\ V_n\left\{2\right\}=\frac{1}{{\varphi }_r^2}\left[{\epsilon }_r{V}_n\left\{4\right\}+\chi V_n\left\{3\right\}+\frac{\kappa }{\Delta t}{\Delta }_t^1\left({\epsilon }_r{V}_n\left\{1\right\}-\chi V_n\left\{2\right\}\right)\right]\end{array}$$ E36

Here the index n denotes the branch number (n = 15) of the node, and the symbol $${\Delta }_t^d$$ represents the d–order time variation of the pulses. The values of the new elements to be connected at the TLM node equivalent circuit are obtained from the discretized form of wave equations (33). The TLM method deals with the total voltage at the node, so it is more useful to write the wave equations in terms of such voltages. For example, the wave equation for the z- component of the electric field become for the mesh {1}:

$$\begin{array}{l}\frac{{\partial }^2{}_{k}V{}_{\text{T}}\left\{1\right\}}{\partial x^2}+\frac{{\partial }^2{}_{k}V{}_{\text{T}}\left\{1\right\}}{\partial y^2}\\ -{\varphi }_r^2{\epsilon }_0{\mu }_0\frac{{\partial }^2{}_{k}V{}_{\text{T}}\left\{1\right\}}{\partial t^2}+\frac{\kappa }{c}\frac{{\partial }^3}{\partial t^3}\left[-\frac{\mu }{{\mu }_0}{}_{k}V{}_{\text{T}}\left\{2\right\}+\chi {}_{k}V{}_{\text{T}}\left\{1\right\}-{}_{k}V{}_{\text{T}}\left\{4\right\}\right]=0\end{array}$$ E37

The voltage of the four sources $$V_s^k$$ at every time iteration k can be obtained directly from these TLM wave equations. The resulting expressions for each TLM mesh are:

$$\begin{array}{c}{V}_s^{k+1}\left\{1\right\}=V_s^{k-1}\left\{1\right\}-\frac{1}{{\varphi }_r^2}\frac{\kappa }{\Delta t}{\left[\gamma \right]}^k\\ V_s^{k+1}\left\{2\right\}=V_s^{k-1}\left\{2\right\}+\frac{1}{{\varphi }_r^2}\frac{\kappa }{\Delta t}{\left[\rho \right]}^k\\ V_s^{k+1}\left\{3\right\}=V_s^{k-1}\left\{3\right\}-\frac{1}{{\varphi }_r^2}\frac{\kappa }{\Delta t}{\left[{\epsilon }_r\gamma +\chi \rho -\frac{\kappa }{\Delta t}{\Delta }_t^1\rho \right]}^k\\ V_s^{k+1}\left\{4\right\}=V_s^{k-1}\left\{4\right\}+\frac{1}{{\varphi }_r^2}\frac{\kappa }{\Delta t}{\left[{\mu }_r\rho +\chi \gamma +\frac{\kappa }{\Delta t}{\Delta }_t^1\gamma \right]}^k\end{array}$$ E38

where the symbols and indicate:

$$\begin{array}{ccc}\gamma ={\Delta }_t^3\left[{\mu }_r{}_{k}V{}_{\text{T}}\left\{2\right\}-\chi {}_{k}V{}_{\text{T}}\left\{1\right\}+{}_{k}V{}_{\text{T}}\left\{4\right\}\right]& & \rho ={\Delta }_t^3\left[{\epsilon }_r{}_{k}V{}_{\text{T}}\left\{1\right\}-\chi {}_{k}V{}_{\text{T}}\left\{2\right\}+{}_{k}V{}_{\text{T}}\left\{3\right\}\right]\end{array}$$ E39

The sources are introduced into the algorithm by adding them to the voltage impulses, before the scattering process, as follows:

$${}_{k}V{}_{\text{T}}\left\{m\right\}=\frac{2}{4+y_0}\left(\sum _{n=1}^4{}_{k}V{}_n^i\left\{m\right\}+y_0{}_{k}V{}_5^i\left\{m\right\}\right)+V_s^k\left\{m\right\}$$ E40

Finally, the information to get out any field vectors component, at any point of the network and at any time k Δt, may be found in Table 4. The voltage $$V_z$$ at the node represents the total voltage V _T (10), and the currents $$I_x$$ and $$I_y$$ are proportional to the voltage pulses difference. For example, the electric field strength components can be expressed in terms of the incident voltage pulses upon the node:

$$\begin{array}{l}\begin{array}{ccc}{}_{k}E{}_x=\frac{{}_{k}V{}_3^i\left\{4\right\}-{}_{k}V{}_1^i\left\{4\right\}}{\Delta \ell \sqrt{2}{\varphi }_r^2}& {}_{k}E{}_y=\frac{{}_{k}V{}_2^i\left\{4\right\}-{}_{k}V{}_4^i\left\{4\right\}}{\Delta \ell \sqrt{2}{\varphi }_r^2}& \end{array}\\ {}_{k}E{}_z=\frac{{\displaystyle \sum _{n=1}^4{}_{k}V{}_n^i\left\{1\right\}}+y_o{}_{k}V{}_5^i\left\{1\right\}}{\Delta \ell 2{\varphi }_r^2}+\frac{V_s^k\left\{1\right\}}{\Delta \ell }\end{array}$$ E41

Finally, all these steps (pulse scattering, node connection, signal transfer between meshes (36), voltage sources updating (38), voltage impulses summation (40), output of the fields values (41)…) are repeated as many time steps as required.

5.2. Results

To validate the TLM model, we have modeled again the plane wave propagation in bi-isotropic media and compared the results with known theoretical ones [1]. In order to simulate a plane wave propagating along the x–direction, the mesh boundaries are modeled with absorbing conditions [62] for limits in x–direction to avoid not desired reflections, and magnetic walls for the y–direction. It is to be noted that in the four meshes (Table 4), the total voltage at the node is either equivalent to an electric component ({1} and {3}) or to a reversed–sign magnetic component ({2} and {4}) of the electromagnetic field, thus the same boundary conditions in the y–direction have been modeled at all the meshes.

The propagation of a plane wave propagation in a bi-isotropic medium with relative permittivity and permeability $${\epsilon }_r$$ = 5.0 and $${\mu }_r$$ =1.5, Tellegen parameter $$\chi$$ = 1.2, and Pasteur parameter (Kong model) $$\kappa$$ = –0.01 ps is now computed. In the simulation, the computational space domain is n _x = 5,000 nodes, and n _y = 10 nodes along the x– and y–axes, respectively, with a space step $$\Delta \ell$$ = 1.00 mm. A monochromatic excitation of frequency 4.00 GHz is applied at the plane x = 500 $$\Delta \ell$$ . The excitation is introduced with linear polarization in the z–direction of the electric field strength $$E_z=E_0\sin \left(\omega t\right)$$ , then:

$${}_{k}V{}_n\left\{1\right\}={}_{k-1}V{}_n\left\{1\right\}+\frac{1}{2}{E}_0\sin \left(\omega t\right)forn=\text{1, }\mathrm{...}\text{, 5}$$ E42

	Theoretical values	TLM results	Relative Error (%)
θEH	115.99º	116.09º	−0.09
θDB	64.01º	64.12º	−0.17
10^2 ηEH/o	54.77	54.63	+0.26
10^2 ηDB/o	54.77	54.92	−0.27
10 εr	44.94	44.87	+0.16
10μr	13.48	13.46	+0.15

Table 6.

Angles, impedances and effective parameters [62]

In Table 6, the angle tilt between electric and magnetic field vectors θ_EH and θ_DB, the intrinsic impedances η_EH and η_DB, and the effective relative permittivity and permeability are shown after 10,000 time iterations ( $$\Delta t$$ = 2.36 ps) of the algorithm. In Table 7, the rotation of the angle of the polarization plane at different distances from the source is presented. The respective theoretical values [62] and relative errors are included also in these two tables of results: they are in good agreement with the theoretical values, with relative errors less than 0.3 %.

Distance to the source (Δℓ)	Theoretical values	TLM results	Relative Error (%)
500	90.60º	90.74º	−0.15
1000	91.21º	91.36º	−0.17
1500	91.81º	91.97º	−0.18
2000	92.41º	92.60º	−0.20
2500	93.02º	93.26º	−0.26

Table 7.

Rotation of the polarization plane [62]

6. Frequency Domain Modelling: The Coupled Mode Method (CMM)

The Coupled Mode Method is one of the available numerical methods in the frequency domain used for analyzing the electromagnetic wave propagation inside waveguides that contain any isotropic, anisotropic or complex medium [63-69]. Originally, this method was formulated by Schelkunoff [63] for uniform closed structures involving heterogeneous and non-isotropic media. Afterwards, it was extended to cover other structures such us open dielectric guides [64], waveguides containing bi-isotropic media [65-69], and free propagation in unbounded bi-isotropic media [70].

Moreover, the application of this numerical method is not only restricted to the analysis of uniform structures. In fact, one of the main advantages of the coupled mode method is that it can be combined with the Mode Matching Method (MMM) in order to analyze discontinuities between two waveguides that contain isotropic, anisotropic or complex media. Thus, joining the coupled mode method and the mode matching method a hybrid numerical technique capable of analyzing 3D structures implemented in a waveguide that presents discontinuities in the energy propagation direction is obtained [71-72]. Due to the mode matching method properties, there is no limitation in the number of discontinuities that we can analyze. Therefore 3D periodic structures can also be examined with this numerical tool.

6.1. Theory

In few words, the coupled mode method is a method of moments [63] which consists on expanding the electromagnetic field components inside a uniform waveguide partially filled with any of the former media in terms of a set of base functions previously defined. Thus, the general expression of any electromagnetic field component is expressed as

$$G(x,y,z)={\displaystyle \sum _{n=1}^{\infty }{c}_n{g}_n(x,y,z)}={\displaystyle \sum _{n=1}^{\infty }{c}_n{\tilde{g}}_n(x,y)}{e}^{-\Gamma z}$$ E43

where $$c_n$$ are the expansions coefficients, $${\tilde{g}}_n$$ the selected base functions and the propagation constant of the proper modes, i.e., the modes of the partially filled waveguide. In many cases these base functions are the electric and magnetic fields of the different modes supported by the basis structure, i.e., the empty waveguide. Therefore, they are called basis modes. When the guiding structure is a rectangular waveguide these basis modes correspond to the TE and TM modes of the empty rectangular waveguide [63].

The obtention of the expansions coefficients and the propagation constants of the proper modes constitutes the mathematical procedure of the coupled mode method. Mainly, it consists on selecting two electromagnetic field vectors as primary expansions

As it will be shown later these primary fields can be $\overrightarrow{E}$ and $\overrightarrow{H},\overrightarrow{D}$ or $\overrightarrow{B}$ and $\overrightarrow{E}$ .

and then substituting them into the x- and y- components of the Maxwell Equations.

Next, the other two field vectors are expressed in terms of the primary ones through the constitutive relations (3). Thus, the obtained equations only depend of the primary fields expansions. In order to solve the problem, the Galerkin method is applied [73] and a set of differential equations, named the Generalized Telegraphist Equations, is obtained [63].

Finally, after some analytical work involving the z-components of the Maxwell Equations [63,69] and the constitutive relations, we reach to an eigenvalue problem whose solution provides the expansion coefficients and the propagation constants of the proper modes

$$\left\{\left[A\right]-\Gamma \left[U\right]\right\}\left[C\right]=0.$$ E44

In Eq. (44) [A] is a coupling matrix between the basis modes taken into account in the procedure, [U] is an unitary matrix, is a vector related to the propagation constants of the electromagnetic field and [C] is a matrix containing the coefficients of the electromagnetic field expansions. In following sections this procedure is analyzed in more detail.

6.1.1. A. The different formulations of the coupled mode method:

As it is know, the electromagnetic field is described by four field vectors $\overrightarrow{H}$ and, $\overrightarrow{B}$ which are related among themselves through the constitutive relations. Thus, only the expansions of two primary fields have to be proposed: one electric ( $\overrightarrow{E}$ or $\overrightarrow{H}$ ) and another magnetic ( $\overrightarrow{E}$ or $\overrightarrow{B}$ ). The other two field vectors are obtained by substituting the primary fields into the constitutive relations. Hence, for the same problem, the coupled mode method can be formulated in different ways giving rise to different formulations that produce similar solutions.

So far, in the literature two alternative formulations of the coupled mode method can be found: the EH- and EB-formulations [67-69]. In the first one the primary fields are $\overrightarrow{E}$ and $\overrightarrow{H}$ , whereas in the second one the primary fields are $\overrightarrow{E}$ and $\overrightarrow{B}$ . As we will see later, the main difference between the EH- and EB-formulations of the coupled mode method resides in the satisfaction or not of all boundary conditions of the electromagnetic field components over the perfect electric conductor that define the waveguide. In particular, the issue is related to the normal component of the magnetic field over the walls in contact with the chiral material, as they must be no null-valued [69].

EH-formulation of the coupled mode method.

In the EH-formulation, the primary expansions, particularized for a uniform rectangular waveguide [63], are

$$E_x\left(x,y\right)=\stackrel{nte}{{\displaystyle \sum }_{}^{}_{q=1}}{V}_{\left[q\right]}\frac{\partial T_{\left[q\right]}\left(x,y\right)}{\partial y}+\stackrel{ntm}{{\displaystyle \sum }_{}^{}_{p=1}}{V}_{\left(p\right)}\frac{\partial T_{\left(p\right)}\left(x,y\right)}{\partial x}$$ E45

$$E_y\left(x,y\right)=-\stackrel{nte}{{\displaystyle \sum }_{}^{}_{q=1}}{V}_{\left[q\right]}\frac{\partial T_{\left[q\right]}\left(x,y\right)}{\partial x}+\stackrel{ntm}{{\displaystyle \sum }_{}^{}_{p=1}}{V}_{\left(p\right)}\frac{\partial T_{\left(p\right)}\left(x,y\right)}{\partial y}$$ E46

$$E_z\left(x,y\right)=\stackrel{ntm}{{\displaystyle \sum }_{}^{}_{p=1}}{k}_{c\left(p\right)}{V}_{\left(p\right)}^z{T}_{\left(p\right)}\left(x,y\right)$$ E47

$$H_x\left(x,y\right)=\stackrel{nte}{{\displaystyle \sum }_{}^{}_{q=1}}{I}_{\left[q\right]}\frac{\partial T_{\left[q\right]}\left(x,y\right)}{\partial x}-\stackrel{ntm}{{\displaystyle \sum }_{}^{}_{p=1}}{I}_{\left(p\right)}\frac{\partial T_{\left(p\right)}\left(x,y\right)}{\partial y}$$ E48

$$H_y\left(x,y\right)=\stackrel{nte}{{\displaystyle \sum }_{}^{}_{q=1}}{I}_{\left[q\right]}\frac{\partial T_{\left[q\right]}\left(x,y\right)}{\partial y}+\stackrel{ntm}{{\displaystyle \sum }_{}^{}_{p=1}}{I}_{\left(p\right)}\frac{\partial T_{\left(p\right)}\left(x,y\right)}{\partial x}$$ E49

$$H_z\left(x,y\right)=\stackrel{nte}{{\displaystyle \sum }_{}^{}_{q=1}}{k}_{c_{\left[q\right]}}{I}_{\left[q\right]}^z{T}_{\left[q\right]}\left(x,y\right)+h_0$$ E50

Here, V’s and I’s represent the electric and magnetic field’s expansion coefficients, and the T-functions are the T-potentials that give the electromagnetic field components of the basis modes TE and TM in the empty waveguide [63,68]. For rectangular waveguides they are

• TE modes

$$\begin{array}{l}{T}_{\left[n,m\right]}=\frac{{\eta }_{\left[n,m\right]}}{k_{c_{\left[n,m\right]}}\sqrt{ab}}\cos \left(\frac{n\pi }{a}x\right)\cos \left(\frac{m\pi }{b}y\right),\\ n,m=0,1,2,\mathrm{....}butnotn=m=0\end{array}$$ E51

• TM modes

$$\begin{array}{l}{T}_{\left(n,m\right)}=\frac{{\eta }_{\left(n,m\right)}}{k_{c_{\left(n,m\right)}}\sqrt{ab}}\sin \left(\frac{n\pi }{a}x\right)\sin \left(\frac{m\pi }{b}y\right)\\ n,m=1,2,3,\mathrm{....}\end{array}$$ E52

with

$${\eta }_{\left\{n,m\right\}}=\{\begin{array}{c}2ifn\ne 0,m\ne 0\\ \sqrt{2}ifn=0orm=0\end{array}$$ E53

Moreover, in Eqs. (45) -(52), k _c(p) and k _c[q] are the cutoff wavenumber of the basis modes; the term h ₀ is a constant that must be taken into account in some media in order to obtain the correct solution of the electromagnetic field. Also, nte and ntm correspond, respectively, to the number of basis modes TE and TM taken into account in the expansions. In all expressions, parentheses refer to TM-modes and brackets to TE-modes.

Following the procedure of the coupled mode method previously mentioned at the beginning of this section and detailed in [63,68], the Generalized Telegraphist Equations for the EH-formulation turn out as:

$$\frac{d{V}_{\left(n\right)}}{dz}=j\omega {{\displaystyle }}_s\left(B_x\frac{\partial T_{\left(n\right)}}{\partial y}-B_y\frac{\partial T_{\left(n\right)}}{\partial x}\right)ds+k_{c_{\left(n\right)}}{V}_{\left(n\right)}^z$$ E54

$$\frac{d{V}_{\left[m\right]}}{dz}=-j\omega {{\displaystyle }}_s\left(B_x\frac{\partial T_{\left[m\right]}}{\partial x}+B_y\frac{\partial T_{\left[m\right]}}{\partial y}\right)ds$$ E55

$$\frac{d{I}_{\left(n\right)}}{dz}=-j\omega {{\displaystyle }}_s\left(D_x\frac{\partial T_{\left(n\right)}}{\partial x}+D_y\frac{\partial T_{\left(n\right)}}{\partial y}\right)ds$$ E56

$$\frac{d{I}_{\left[m\right]}}{dz}=-j\omega {{\displaystyle }}_s\left(D_x\frac{\partial T_{\left[m\right]}}{\partial y}-D_y\frac{\partial T_{\left[m\right]}}{\partial x}\right)ds+k_{c_{\left[m\right]}}{I}_{\left[m\right]}^z$$ E57

Inspecting the set of equations (54)-(57) it can be seen that it is constituted by 4 equations and 6 unknowns. Therefore, in order to find their solution, two additional expressions relating the longitudinal coefficients with the transversal ones are needed. These relations are provided by the longitudinal components of the Maxwell equations and the constitutive relations; they can be obtained following two alternative ways that give two different formulations of the EH-formulation.

EHD-formulation

The first presented way to obtain the additional relations between the longitudinal and transversal coefficients it the one proposed by Ogusu [64]. Although it is not the original procedure, it provides the more intuitive way to relate these coefficients. Therefore, it is firstly introduced.

In this formulation the additional relations are obtained by rewriting the constitutive relations (3) as

$$E_z=\frac{1}{{\epsilon }_0}\frac{{\mu }_r}{\left[{\mu }_r{\epsilon }_r-\left({\chi }^2+{\kappa }^2\right)\right]}{D}_z-c\frac{\left(\chi -j\kappa \right)}{\left[{\mu }_r{\epsilon }_r-\left({\chi }^2+{\kappa }^2\right)\right]}{B}_z$$ E58

$$H_z=c\frac{\left(\chi +j\kappa \right)}{\left[\left({\chi }^2+{\kappa }^2\right)-{\mu }_r{\epsilon }_r\right]}{D}_z-\frac{1}{{\mu }_0}\frac{{\epsilon }_r}{\left[\left({\chi }^2+{\kappa }^2\right)-{\mu }_r{\epsilon }_r\right]}{B}_z$$ E59

Then, after introducing the longitudinal components of the Maxwell equations and the expansions (45)-(50) into (58)-(59), and applying the Galerkin mode [73], we get the desired expressions

$$\begin{array}{c}{k}_{c_{\left(n\right)}}{V}_{\left(n\right)}^z=\frac{c}{j\omega }\stackrel{nte}{{\displaystyle \sum }_{}^{}_{j=1}}{V}_{\left[j\right]}{k}_{c_{\left[j\right]}}^2{k}_{c_{\left(n\right)}}^2{{\displaystyle }}_S\frac{\left(\chi -j\kappa \right)}{\left[{\mu }_r{\epsilon }_r-\left({\chi }^2+{\kappa }^2\right)\right]}{T}_{\left[j\right]}{T}_{\left(n\right)}ds\\ -\frac{1}{j\omega {\epsilon }_0}\stackrel{ntm}{{\displaystyle \sum }_{}^{}_{i=1}}{I}_{\left(i\right)}{k}_{c_{\left(i\right)}}^2{k}_{c_{\left(n\right)}}^2{{\displaystyle }}_S\frac{{\mu }_r}{\left[{\mu }_r{\epsilon }_r-\left({\chi }^2+{\kappa }^2\right)\right]}{T}_{\left(i\right)}{T}_{\left(n\right)}ds\end{array}$$ E60

$$\begin{array}{c}{k}_{c_{\left[m\right]}}{I}_{\left[m\right]}^z=\frac{1}{j\omega {\mu }_0}\stackrel{nte}{{\displaystyle \sum }_{}^{}_{j=1}}{V}_{\left[j\right]}{k}_{c_{\left[j\right]}}^2{k}_{c_{\left[m\right]}}^2{{\displaystyle }}_S\frac{{\epsilon }_r}{\left[\left({\chi }^2+{\kappa }^2\right)-{\mu }_r{\epsilon }_r\right]}{T}_{\left[j\right]}{T}_{\left[m\right]}ds\\ -\frac{c}{j\omega }\stackrel{ntm}{{\displaystyle \sum }_{}^{}_{i=1}}{I}_{\left(i\right)}{k}_{c_{\left(i\right)}}^2{k}_{c_{\left[m\right]}}^2{{\displaystyle }}_S\frac{\left(\chi +j\kappa \right)}{\left[\left({\chi }^2+{\kappa }^2\right)-{\mu }_r{\epsilon }_r\right]}{T}_{\left(i\right)}{T}_{\left[m\right]}ds\end{array}$$ E61

Note that in (60)-(61) the longitudinal coefficients are directly related with the transversal ones, i.e., no matrix inversion is required. Therefore we call it direct EH-formulation of the coupled mode method or, simply, EHD-formulation.

Once the longitudinal and transversal coefficients are known, only the coefficient h ₀ have to be determined. Following a similar procedure, the value of h ₀ can be obtained from

$$\begin{array}{l}{h}_0{{\displaystyle }}_{S}ds=\\ \frac{1}{j\omega {\mu }_0}\stackrel{nte}{{\displaystyle \sum }_{}^{}_{j=1}}{V}_{\left[j\right]}{k}_{c_{\left[j\right]}}^2{{\displaystyle }}_S\frac{{\epsilon }_r}{\left[\left({\chi }^2+{\kappa }^2\right)-{\epsilon }_r\right]}{T}_{\left[j\right]}ds-\frac{c}{j\omega }\stackrel{ntm}{{\displaystyle \sum }_{}^{}_{i=1}}{I}_{\left(i\right)}{k}_{c_{\left(i\right)}}^2{{\displaystyle }}_S\frac{\left(\chi +j\kappa \right)}{\left[\left({\chi }^2+{\kappa }^2\right)-{\epsilon }_r\right]}{T}_{\left(i\right)}ds\end{array}$$ E62

EHI-formulation

The alternative process to obtain the relations between the longitudinal and transversal coefficients is the original one formulated by Schelkunoff [63]. Now, the longitudinal components of the Maxwell equations and the expansions (45)-(50) are directly introduced in the constitutive relations (3). After applying the Galerkin method [73], the desired equations are obtained

$$\begin{array}{c}{I}_{\left(n\right)}=-j\omega {\epsilon }_0\stackrel{ntm}{{\displaystyle \sum }_{}^{}_{i=1}}{k}_{c_{\left(i\right)}}{V}_{\left(i\right)}^z{{\displaystyle }}_S{\epsilon }_r{T}_{\left(i\right)}{T}_{\left(n\right)}ds\\ -j\frac{\omega }{c}\stackrel{nte}{{\displaystyle \sum }_{}^{}_{j=1}}{k}_{c_{\left[j\right]}}{I}_{\left[j\right]}^z{{\displaystyle }}_S\left(\chi -j\kappa \right)T_{\left[j\right]}{T}_{\left(n\right)}ds\\ -j\frac{\omega }{c}{h}_0{{\displaystyle }}_S\left(\chi -j\kappa \right)T_{\left(n\right)}ds\end{array}$$ E63

$$\begin{array}{c}{V}_{\left[m\right]}=-j\omega {\mu }_0\stackrel{nte}{{\displaystyle \sum }_{}^{}_{j=1}}{k}_{c_{\left[j\right]}}{I}_{\left[j\right]}^z{{\displaystyle }}_S{\mu }_r{T}_{\left[j\right]}{T}_{\left[m\right]}ds-j\omega {\mu }_0{h}_0{{\displaystyle }}_S{\mu }_r{T}_{\left[m\right]}ds\\ -j\frac{\omega }{c_0}\stackrel{ntm}{{\displaystyle \sum }_{}^{}_{i=1}}{k}_{c_{\left(i\right)}}{V}_{\left(i\right)}^z{{\displaystyle }}_S\left(\chi +j\kappa \right)T_{\left(i\right)}{T}_{\left[m\right]}ds\end{array}$$ E64

Inspecting (62) and (63) two main differences with the EHD-formulation can be seen. The first one is that in order to isolate the longitudinal coefficients, matrix inversions are required. Therefore, as the longitudinal coefficients cannot be directly expressed in terms of the transversal ones, we call this formulation “Indirect EH-formulation” of the coupled mode method or, simply, EHI-formulation. The second difference deals with the obtention of the constant h ₀. As it influences over the expansions coefficients, it must be taken into account in the eigenvalue problem. Thus a third relation must be included. Following a similar procedure, the value of h ₀ can be obtained from

$$h_0{{\displaystyle }}_S{\mu }_{r}ds=-\stackrel{nte}{{\displaystyle \sum }_{}^{}_{j=1}}{k}_{c_{\left[m\right]}}{I}_{\left[m\right]}^z{{\displaystyle }}_S{\mu }_r{T}_{\left[m\right]}ds-\sqrt{\frac{{\epsilon }_0}{{\mu }_0}}\stackrel{ntm}{{\displaystyle \sum }_{}^{}_{i=1}}{k}_{c\left(i\right)}{V}_{\left(i\right)}^z{{\displaystyle }}_S\left(\chi +j\kappa \right)T_{\left(i\right)}ds$$ E65

EB-formulation

We have seen that, in the EH formulation of the coupled mode method, both electric and magnetic fields inside the partially filled waveguide are expressed in terms of the electric and magnetic fields of the empty waveguide. Therefore they satisfy the same boundary conditions over the perfect electric conductor. As long as the media contained in the waveguide allow the proper modes and the basis modes to satisfy the same boundary conditions, the EH-formulation provides correct results. This happens when the media are dielectric, isotropic or anisotropic. The problem appears when the medium produces a discontinuity in any component of the electromagnetic field that cannot be reproduced with the basis modes of the waveguide. This occurs, for example, when a bi-isotropic medium fills the waveguide. If we focus our attention on the normal component of the magnetic field over the perfect electric conductor, we see that it must be not null valued, whereas the normal component of the magnetic induction does. However, when the problem is analyzed with the coupled mode method, it can be seen that, if the magnetic field satisfies the same boundary conditions than the basis modes, the normal component of the magnetic field will be zero over the PEC, which is incorrect. In order to avoid this inconvenience we can choose as primary fields those that satisfy the same boundary conditions over the perfect electric conductors than the basis modes fields. These are $\overrightarrow{E}$ and $\overrightarrow{B}$ . Therefore, the proposed primary expansions for these fields are

$$E_x=\stackrel{nte}{{\displaystyle \sum }_{}^{}_{m=1}}{V}_{\left[m\right]}\frac{\partial T_{\left[m\right]}}{\partial y}+\stackrel{ntm}{{\displaystyle \sum }_{}^{}_{n=1}}{V}_{\left(n\right)}\frac{\partial T_{\left(n\right)}}{\partial x}$$ E66

$$E_y=-\stackrel{nte}{{\displaystyle \sum }_{}^{}_{m=1}}{V}_{\left[m\right]}\frac{\partial T_{\left[m\right]}}{\partial x}+\stackrel{ntm}{{\displaystyle \sum }_{}^{}_{n=1}}{V}_{\left(n\right)}\frac{\partial T_{\left(n\right)}}{\partial y}$$ E67

$$E_z=\stackrel{ntm}{{\displaystyle \sum }_{}^{}_{n=1}}{k}_{c_{\left(n\right)}}{V}_{\left(n\right)}^z{T}_{\left(n\right)}$$ E68

$$B_x={\mu }_0\stackrel{nte}{{\displaystyle \sum }_{}^{}_{m=1}}{b}_{\left[m\right]}\frac{\partial T_{\left[m\right]}}{\partial x}-{\mu }_0\stackrel{ntm}{{\displaystyle \sum }_{}^{}_{n=1}}{b}_{\left(n\right)}\frac{\partial T_{\left(n\right)}}{\partial y}$$ E69

$$B_y={\mu }_0\stackrel{nte}{{\displaystyle \sum }_{}^{}_{m=1}}{b}_{\left[m\right]}\frac{\partial T_{\left[m\right]}}{\partial y}+{\mu }_0\stackrel{ntm}{{\displaystyle \sum }_{}^{}_{n=1}}{b}_{\left(n\right)}\frac{\partial T_{\left(n\right)}}{\partial x}$$ E70

$$B_z={\mu }_0\stackrel{nte}{{\displaystyle \sum }_{}^{}_{m=1}}{k}_{c_{\left[m\right]}}{b}_{\left[m\right]}^z{T}_{\left[m\right]}+{\mu }_0{b}_0$$ E71

Following a similar procedure than the presented for the conventional formulation of the coupled mode method, an eigenvalue problem of Eq. (54)-(57) type is obtained. As it happened in the previous formulations, in order to solve this set of equations, two new equivalent relations between longitudinal and transversal coefficients and another third one to obtain b ₀ are needed. When the medium filling the waveguide is an isotropic chiral medium, the coefficient b ₀ turns to be null.

Solving the Telegraphist equation eigenvalue problem, the coefficients expansions and the propagation constants for $\overrightarrow{H}$ and $\overrightarrow{H}$ can be obtained in a way similar to the EH-formulation.

Magnetic Field

When only the knowledge of the electric field or the constant propagations of any proper modes is required, the expansions of $\overrightarrow{H}$ and $\overrightarrow{E}$ obtained using the EB- formulation are enough to obtain the desired solution. However, if we need to characterize a discontinuity, an expansion for $\overrightarrow{H}$ is also necessary. As it was previously mentioned, the obtention of this new expansion for $\overrightarrow{E}$ involves Eqs. (66)-(71) and the constitutive relations (3). Once the $\overrightarrow{B}$ expansions are known, the mode matching method can be easily extended [71-72, 75] and any discontinuity can be analyzed.

6.2. Results

In order to check the behavior of the different formulations of the coupled mode method previously presented, the analysis of electromagnetic wave propagation inside a uniform rectangular waveguide filled with an isotropic chiral medium is proposed (Figure 5). It is important to note that there is no analytical solution available of this problem. Hence, a priori one could not establish which formulation provides the most accurate results. However, according to the conclusions extracted from the analysis of parallel-plate waveguide containing an isotropic chiral media [67], and the fact that the EB-formulation was developed to ensure the satisfaction of all boundary conditions of the electromagnetic field over the perfect electric conductor, we can infer that the EB- is the optimal formulation for analyzing this problem. Thus, along this section this statement is verified with numerical results. For that purpose we focus our attention in two properties of the electromagnetic field: the propagation constant values of the proper modes and their magnetic field components profile. In the first case it can be seen how the convergence behavior of each formulation is. In the second one it is shown how the magnetic field obtained by the EB-formulation satisfies the boundary conditions of the normal component of the magnetic field whereas that extracted from the EH-formulation does not.

Figure 6 shows the variation of the phase constant of the fundamental proper mode normalized to the vacuum wavenumber versus the number of basis TE and TM modes (convergence diagram) for the structure of Figure 5 provided by the mentioned formulations of the coupled mode method. Here, the basis modes TE and TM have been taken by increasing their cutoff frequency and the simulation data is displayed in the figure caption. An inspection of this figure evidences that the final values of the propagation constant provided by the EB- and EHI-formulations are very similar but slightly different than the EHD ones. Moreover, the convergences of the EB- and EHI- are faster than the EHD- one as they need a smaller number of basis modes TE and TM to reach to the final value. Therefore, with these evidences, which can be extended to the rest of proper modes, we can infer that the EHI-formulation gives more accurate results than the EHD-.

Figure 5.

Cross section of a uniform rectangular waveguide fully loaded with an isotropic chiral media (ICM). P.E.C. stands for perfect electric conductor

Figure 6.

Convergence diagram of the fundamental mode for a uniform waveguide completely filled with an isotropic chiral media. Dimensions in mm: a × b = 22.86 × 10.16;ε_r = 3, µ_r = 1, χ = 0 and κ= 1.5.

As it was previously mentioned, since the EHI-formulation implies a matrix inversion that introduces a numerical error, one could think that the results provided by the EHD-formulation should be more accurate than those from the EHI. However, this statement is false and the EHI- results are closer to the EB- than the EHD- ones. The explanation of this disagreement is related to the fact that in the EHD-formulation the continuous E _z component is written as the product of discontinuous functions. This situation, which is right analytically, is not correct when is obtained numerically, as it is impossible to reproduce a discontinuous function in terms of a finite number of sines and cosines. In the other hand, in the EHI-formulation, although matrix inversions are required, a discontinuous function is expressed in terms of the product of discontinuous and continuous functions, which is correct both analytically and numerically. For this reason, the EHI-formulation results agree with the EB-formulation ones better than those obtained from the EHD-formulation.

Figure 7.

Module of the magnetic field components of the fundamental mode, normalized to their maximum value, as a function of the x- and y- coordinates provided by (a) EHI-formulation and (b) EB-formulation. Dimensions in mm: a × b = 22.86 × 10.16. ε_r = 3, µ_r = 1, χ= 0 and κ= 0.75. Number of basis modes 120 for the EB-formulation and 150 for the EH-.

Finally, in Figure 7 the module of the three components of the magnetic field normalized to their maximum values, as a function of the x- and y- coordinates are presented. As it is clear that the EHI-formulation is more appropriate that the EHD when analyzing this problem, only results for EHI- and EB-formulations are shown. On inspecting Figure 7 it is easy to prove that the magnetic field obtained with the EB-formulation does satisfy the boundary conditions of the normal component of the magnetic field whereas that extracted from the EH-formulation does not. It can be seen in Figure 7 how the normal component of the magnetic field over the perfect electric conductor walls is null in the EHI-formulation meanwhile it is not null in the EB- one.

7. Conclusions

A review of different numerical techniques for the modeling of electromagnetic wave propagation through bi-isotropic materials has been presented. In the time domain, different FDTD models have been mentioned, with emphasis in one based on Mobius transformation, as well as one MRTD approach, suitable to improve recursive convolution techniques and one TLM model. In the frequency domain, two different formulations of the coupled mode method, the EH- and the EB-formulations, have been presented. It has been seen that the EB-formulation is the optimal for analyzing electromagnetic wave propagation inside a rectangular waveguide that contains an isotropic chiral media because the solution that it gives verify all the boundary conditions of the electromagnetic field over the perfect electric walls that constitute the waveguide.

Acknowledgments

Partial support for this work was provided by the Spanish MICINN through the TEC2010-21496-C03-01, subprogram “Ramón y Cajal” (RYC-2010-06922) and CONSOLIDER CSD2008-00066 projects.