InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Engineering » Electrical and Electronic Engineering » "Solutions and Applications of Scattering, Propagation, Radiation and Emission of Electromagnetic Waves", book edited by Ahmed Kishk, ISBN 978-953-51-0838-2, Published: November 14, 2012 under CC BY 3.0 license. © The Author(s).

Chapter 3

Numerical Modeling of Electromagnetic Wave Propagation Through Bi-Isotropic Materials

By I. Barba, A. Grande, A.C.L. Cabeceira, A. Gómez, J.A. Pereda and J. Represa
DOI: 10.5772/50680

Article top


Transformation of the constitutive equations into the discrete-time domain.
Figure 1. Transformation of the constitutive equations into the discrete-time domain.
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.
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.
Reflection and transmission coefficients: theoretical values and relative errors.
Figure 3. Reflection and transmission coefficients: theoretical values and relative errors.
Haar wavelet functions: left, scale function Φ
						(t). Right, first order wavelet function ψ
						(t). Both cases are particularized for n = 1.
Figure 4. Haar wavelet functions: left, scale function Φ n (t). Right, first order wavelet function ψ n (t). Both cases are particularized for n = 1.
Cross section of a uniform rectangular waveguide fully loaded with an isotropic chiral media (ICM). P.E.C. stands for perfect electric conductor
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
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.
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.
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-.
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-.

Numerical Modeling of Electromagnetic Wave Propagation Through Bi-Isotropic Materials

I. Barba1, A. Grande1, A.C.L. Cabeceira1, A. Gómez2, J.A. Pereda2 and J. Represa1

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:


In which we assume the medium is linear. D,E,B , and 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:


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:


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 (Fe2TeO6 and Cr2O3) 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]:


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 ω 0 [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


and the constitutive relations for chiral media are


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


being ε ( μ ) the permittivity (permeability) at infinite frequency, εs ( μs ) the permittivity (permeability) at zero frequency, ωe,ωh and ωκ the resonance frequencies, δe,δh and δκ the loss factors and τκ a time constant. Introducing (6) into (5) we obtain


where the terms Jhh,Jee,Jeh and Jhe are auxiliary current densities defined as


with σh=s(μμ),  ςh=sκ^/cσe=s(εε) and ςe=ςh . Jhh,Jeh and Jee,Jhe are evaluated at the same point and at the same time step than H and E , respectively.

B Difference model

A discrete time-domain version of (10) is obtained by applying the property sF(s)dF(t)/dt . Then, using central differences for the derivatives, average in time for Jhh,Jee and average in space for Jhe,Jeh , the following expressions are obtained [50]


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(s)dF(t)/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


where σh(Z) is the Z-domain magnetic conductivity. To obtain σh(Z) we adopt Mobius transformation technique [51]. This approach involves the following change of variable


Applying this change to σh(Z) yields


where cm(σh) and dm(σh) are real-valued coefficients. Now considering the Z-transform property ZmF(Z)Fnm equation (17) can be written in difference form, as follows


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


These equations are evaluated at z=(k+(1/2))Δz . Note, that they are coupled to (15). Decoupling them an eliminating Whh,2 we obtain

Hn+12(k+12)=12μ+Δtc0(σh)(2μHn12(k+12)            Δt{2[×E]n(k+12)+Whhn12(k+12)+Jhhn12(k+12)+Jhen(k)+Jhen(k+1)})



where WhhWhh,1 . Using the same procedure to discretize (12), we obtain the following equation, which is evaluated at z=(k+(1/2))Δz


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

En+1(k)=12ε+Δtc0(σe)(2εEn(k)            +Δt{2[×H]n+12(k)Ween(k)Jeen(k)Jehn+12(k12)Jehn+12(k+12)})

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


Equations (25) and (26) are evaluated at z=kΔz .

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

  1. Hn+(1/2) is calculated by means of (21)

  2. Jhhn+(1/2) and Whhn+(1/2) are updated using (22), where Hn+(1/2) is known from step 1.

  3. Jehn+(1/2) and Wehn+(1/2) are obtained by using (23), where Hn+(1/2) has been obtain in 1.

  4. En+1 is evaluated by using (24)

  5. Jeen+1 and Ween+1 are calculated by means of (25), where En+1 is known from step 4.

  6. Jhen+1 and When+1 are updated by using (26), where En+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

R¯¯=(Rco00Rco)      T¯¯=(TcoTcrTcrTco)

being Rco the copolarized reflection coefficient and Tco,Tcr 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 εs=4.2ε0, ε=3.3ε0, μs=1.1μ0, μ=μ0, ωe=ωh=ωκ=2π×7.5GHz , δe=0.05ωe, δh=0.05ωh, δκ=0.07 and τκ=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 Δt=1ps and the space discretization Δz=0.3mm . The thickness of the slab was d=30Δ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 ( E,D,B,H ) in scale and wavelet functions:


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 ( Fniϕ ), and a second one for wavelet coefficients ( Fniψ ).

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 D,B ):


where κ’(τ) is the inverse Laplace transform of the chirality parameter in frequency domain κ(ω). In this equation Enix,ϕ , for example, would represent the coefficient for the scale function of the x-component of the electric field 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:


where κ'nmϕ=(nm)Δt(nm+1)Δtκ'(τ)dτ , κ'nmψ=(nm+1/2)Δt(nm+1)Δtκ'(τ)dτ(nm)Δt(nm+1/2)Δtκ'(τ)dτ .If we compare these equations with (27), in order to separate scale and wavelet coefficients of D x and B x , we find:


As we may see, the convolutions sums appear in eq. (30), i.e., they affect only the coefficients of the scale functions of D and B ( Dnix,ϕ and Bnix,ϕ ): 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 D and B , using the previous values of E and H .

  2. Updating of the values of E and 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ε0, μ=μ0, and κ following the Condon model: ω 0=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 valueFDTDMRTD

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 Φ2=εμχ2ε0μ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.

{1} Ezx=Byt, Ezy=Bxt,
{2} Hzx=Dyt,
{3} Dzx=Φ2Hyt+κc2Qeyt2, Dzy=Φ2Hxtκc2Qext2, HyxHxy=Dzt
{4} Bzx=Φ2Eytκc2Qmyt2, Bzy=Φ2Ext+κc2Qmxt2, EyxExy=Bzt

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:


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

Ez Qmz 0
Hz Qez 0
Dz εQmzχcQez Qez
Bz χcQmzμQez Qmz

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 Qm and Qe (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, Δ , in a fixed time step Δ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:


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 .

{1} {2} {3} {4}
Vz... Ez μ0cHz Dz/ε0  cBz
Ix... By/μ0 cDy Φr2Hy Φr2ε0cEy
Iy... Bx/μ0 cDx Φr2Hx Φr2ε0cEx

Table 4.

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

{1} Vzx=μ0Ixt, Vzy=μ0Iyt, Ixx+Iyy=φ2μ0Vzt+κη02Qmzt2
{2} Vzx=μ0Ixt, Vzy=μ0Iyt, Ixx+Iyy=φ2μ0Vzt+κ2Qezt2
{3} Vzx=μ0Ixt+κη02Qeyt2, Vzy=μ0Iytκη02Qext2, Ixx+Iyy=φr2ε0Vzt
{4} Vzx=μ0Ixtκ2Qmyt2, Vzy=μ0Iyt+κ2Qmxt2, Ixx+Iyy=φr2ε0Vzt

Table 5.

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

The corresponding TLM equivalences of the vectors Qm and Qe , introduced in (32), are:


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 yo= 4( ϕr21 ). Then, propagation velocities into each TLM network depend on ϕ , 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 ZTLM  equal ηo/ϕr , where η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 Vsk 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:

Vn{3}=εrVn{1}χVn{2}+κΔt Δt1Vn{2} Vn{4}=μrVn{2}χVn{1}κΔt Δt1Vn{1}Vn{1}=1ϕr2[μrVn{3}+χVn{4}κΔt Δt1(μrVn{2}χVn{1})]Vn{2}=1ϕr2[εrVn{4}+χVn{3}+κΔt Δt1(εrVn{1}χVn{2})]

Here the index n denotes the branch number (n = 15) of the node, and the symbol Δtd 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}:


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


where the symbols and indicate:


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


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 Vz at the node represents the total voltage V T (10), and the currents Ix and Iy 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:


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 εr = 5.0 and μr =1.5, Tellegen parameter χ = 1.2, and Pasteur parameter (Kong model) κ = –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 Δ = 1.00 mm. A monochromatic excitation of frequency 4.00 GHz is applied at the plane x = 500 Δ . The excitation is introduced with linear polarization in the z–direction of the electric field strength Ez=E0sin(ωt) , then:

Vkn{1}=Vk1n{1}+12E0sin(ωt) forn=1, ..., 5
Theoretical values TLM resultsRelative Error (%)
θEH 115.99º116.09º−0.09
θDB 64.01º64.12º−0.17
102 ηEH/o 54.7754.63+0.26
102 ηDB/o 54.7754.92−0.27
10 εr 44.9444.87+0.16
10μr 13.4813.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 ( Δ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 valuesTLM resultsRelative Error (%)

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


where cn are the expansions coefficients, 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[1] - 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


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 H and, B which are related among themselves through the constitutive relations. Thus, only the expansions of two primary fields have to be proposed: one electric ( E or H ) and another magnetic ( E or 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 E and H , whereas in the second one the primary fields are E and 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


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

T[n,m]=η[n,m]kc[n,m]abcos(nπax)cos(mπby),n,m=0,1,2,.... but not n=m=0
  • TM modes



η{n,m}={2   if   n0, m02   if   n=0  or m=0

Moreover, in Eqs. (45) -(52), k c(p) and k c[q] are the cutoff wavenumber of the basis modes; the term h 0 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:


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.


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


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


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 0 have to be determined. Following a similar procedure, the value of h 0 can be obtained from



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


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 0. 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 0 can be obtained from



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 E and B . Therefore, the proposed primary expansions for these fields are


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 0 are needed. When the medium filling the waveguide is an isotropic chiral medium, the coefficient b 0 turns to be null.

Solving the Telegraphist equation eigenvalue problem, the coefficients expansions and the propagation constants for H and 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 H and E obtained using the EB- formulation are enough to obtain the desired solution. However, if we need to characterize a discontinuity, an expansion for H is also necessary. As it was previously mentioned, the obtention of this new expansion for E involves Eqs. (66)-(71) and the constitutive relations (3). Once the 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.


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.


1 - I. V. Lindell, A. H. Sihvola, S. A. Tretyakov, A. J. Viitanen, 1994 Electromagnetic waves in Chiral and Bi-Anisotropic Media. Boston Artech House
2 - A. H. Sihvola, 1994 Electromagnetic Modeling of Bi-Isotropic Media.Progress in Electromagnetic Research, PIER 9 45 86
3 - A. H. Sihvola, I. Lindell, 1995 Material effects on bi-anisotropic electromagnetics.IEICE Trans. Electron. (Tokyo) E78-C: 13831390
4 - S. A. Tretyakov, A. H. Sihvola, A. A. Sochava, C. R. Simovski, 1998 Magnetoelectric interactions in bi-anisotropic media.J. Electromag. Waves Appl. JEMWA 12 481 497
5 - C. M. Soukoulis, M. Wegener, 2011 Past achievements and future challenges in the development of three-dimensional photonic metamaterials. Nature Photonics 5 523 530
6 - I. Barba, A. C. L. Cabeceira, A. J. García-Collado, G. J. Molina-Cuberos, J. Margineda, J. Represa, 2011Quasi-planar Chiral Materials for Microwave Frequencies. In: Kishk A, editor. Electromagnetic Wave Propagation in Complex Matter. Rijeka: Intech97 116
7 - J. B. Pendry, 2004 A Chiral Route to Negative Refraction.Science 306 1353 1355
8 - S. A. Tretyakov, A. H. Sihvola, L. Jylhä, 2005 Backward-wave regime and negative refraction in chiral composites.Photonics and Nanostructures: Fundamentals and Applications, 3 107 115
9 - E. Plum, J. Zhou, J. Dong, V. A. Fedotov, T. Koschny, C. M. Soukoulis, N. I. Zheludev, 2009 Metamaterial with negative index due to chirality, Physical Review B 79 035407 1
10 - J. Zhou, J. Dong, W. Wang, T. Koschny, M. Kafesaki, C. M. Soukoulis, 2009, Negative refractive index due to chirality.Physical Review B 79 121104 1
11 - G. J. Molina-Cuberos, A. J. García-Collado, I. Barba, J. Margineda, 2011 Chiral Metamaterials With Negative Refractive Index Composed by an Eight-Cranks Molecule. IEEE Antennas and Wireless Propagation Letters 10 1488 1490
12 - A. Lakhtakia, 1994 Ten years past the Post. In McCall MW, Dewar G, editors. Complex Mediums V: Light and Complexity. Proc. SPIE (Int. Soc. for Optical Engineering), 5508 85 94
13 - W. S. Weiglhofer, 1994 On a medium constraint arising directly from Maxwell’s equations.J. Phys. A27: L871 L874
14 - E. J. Post, 1962 Formal Structure of Electromagnetics. Amsterdam North-Holland.
15 - A. Lakhtakia, W. S. Weiglhofer, 1995 On a constraint on the electromagnetic constitutive relations of nonhomogeneous linear media.IMA J. of Appl. Math. 54 301 306
16 - A. Lakhtakia, W. S. Weiglhofer, 1996 Lorentz covariance, Occam’s razor, and a constraint on linear constitutive relations. Phys. Lett. A 213 107 111Erratum A222: 459.
17 - W. S. Weiglhofer, A. Lakhtakia, 1998 The Post constraint revisited Arch. Elektron. Uebertrag. 52 276 279
18 - A. Lakhtakia, 2004 On the genesis of the Post constraint in modern electromagnetism. Optik-International Journal of Light and Electron Optics 115 151 158
19 - R. E. Raab, A. H. Sihvola, 1997 On the existence of linear non-reciprocal bi-isotropic (NRBI) media. JPhys. A: Math. Gen. 30 1335 1344
20 - I. V. Lindell, A. H. Sihvola, 2004 Perfect electromagnetic conductor. Journal of Electromagnetic waves and Applications JEMWA 19 861 869
21 - I. V. Lindell, K. H. Wallén, 2004 Differential-form electromagnetics and bi-anisotropic Q-media.J. Electromagn. Waves and Appl. JEMWA 18 957 968
22 - I. V. Lindell, K. H. Wallén, 2004 Generalized Q-media and field decomposition in differential-form approach.J. Electromagn. Waves and Appl. JEMWA 18 1045 1056
23 - I. V. Lindell, A. H. Sihvola, 2005 Transformation method for problems involving perfect electromagnetic conductor (PEMC). IEEE Trans.Antennas Propag. 53 3005 3011
24 - S. A. Tretyakov, S. I. Maslovski, I. S. Nefedov, A. J. Viitanen, P. A. Belov, A. Sanmartin, 2003 Artificial Tellegen particle, Electromagnetics 23 665 680
25 - O. L. de Lange, R. E. Raab, 2001 Post’s constraint for electromagnetic constitutive relations, J. Optics A3L23 L26
26 - A. N. Serdyukov, I. V. Semchenko, S. A. Tretyakov, A. H. Sihvola, 2001 Electromagnetics of Bi-anisotropic Materials: Theory and Applications. Amsterdam Gordon and Breach.
27 - J. R. Canto, C. R. Paiva, A. M. Barbosa, 2011 Modal Analysis of Bi-isotropic H-GuidesProgress In Electromagnetics Research PIER 111 1 24
28 - S. A. Matos, C. R. Paiva, A. M. Barbosa, 2009 A Spacetime Algebra Approach to Moving Bi-isotropic Media. Proc IEEE AP-S/URSI International Symp. APSURSI 09.
29 - F. W. Hehl, Y. N. Obukhov, 2005 Linear media in classical electrodynamics and the Post constraint. Phys. Lett. A 334 249 259
30 - E. O. Kamenetskii, M. Sigalov, R. Shavit, 2009 Tellegen particles and magnetoelectric metamaterialsJournal of Applied Physics 105 013537 1
31 - J. A. Kong, 2008 Electromagnetic Wave Theory.Cambridge, USA EMW Publishing1016.
32 - E. U. Condon, 1937 Theories of optical rotatory power. Rev. Modern. Phys. 9 432 457
33 - G. Kristensson, 1999 Condon’s model on optical rotatory power and causality- a scientific trifle. In: Eds. A. Karlsson and G. Kristensson. Festskrift till Staffan Ström: Lund, Sweden, KF Sigma. 97 119 Available: Accesed 2012 May 30.
34 - A. V. Rogacheva, V. A. Fedotov, A. S. Schwanecke, N. I. Zheludev, 2006 Giant Gyrotropy due to Electromagnetic-Field Coupling in a Bilayered Chiral Structure. Physical Review Letters 97 177401 1
35 - C. A. Balanis, 1989 Advanced Engineering Electromagnetics. New York, USA Wiley
36 - A. J. Bahr, K. R. Clausing, 1994 An approximate model for artificial chiral material. IEEE Trans. Antennas Propag. 42 1592 1599
37 - K. S. Yee, 1966 Numerical solution of initial boundary value problems involving Maxwell’s equation is isotropic media. IEEE Trans. Antennas Propag. 14 302 307
38 - A. Taflove, 1995 Computational Electrodynamics: The Finite-Difference Time-Domain Method. Norwood, MA Artech House.
39 - S. García, I. Perez, R. Martín, B. Olmedo, 1998 Extension of Berenger’s PML for bi-anisotropic media. IEEE Microwave Guided Wave Lett. 8 297 299
40 - F. Teixeira, W. Chew, 1998 General closed-form PML constitutive tensors to match arbitrary bi-anisotropic and dispersive linear media. IEEE Microwave Guided Wave Lett. 8 223 225
41 - S. Garcia, I. Perez, R. Martin, B. Olmedo, 1999 BiPML: A PML to match waves in bi-anisotropic media. Microwave Opt. Technol. Lett. 20 44 48
42 - A. Akyurtlu, D. H. Werner, 2004 BI-FDTD: A novel finite-difference time-domain formulation for modeling wave propagation in bi-isotropic media. IEEE Trans. Antennas Propag. 52 416 425
43 - A. Akyurtlu, D. H. Werner, 2004 A novel dispersive FDTD formulation for modeling transient propagation in chiral Metamaterials IEEE Trans. Antennas Propag. 52 2267 2276
44 - I. Barba, A. Grande, A. C. L. Cabeceira, J. Represa, 2003 Time domain modeling of electromagnetic wave propagation in Tellegen media. Microwave Opt. Technol. Lett. 38 167 168
45 - A. Grande, I. Barba, A. C. L. Cabeceira, J. Represa, K. Kärkkäinen, A. H. Sihvola, 2004 Analyzing the stability of the new FDTD technique for non-dispersive bi-isotropic media using the von Neumann method. VIII ACES Symp. Syracuse. New York.
46 - A. Grande, I. Barba, A. C. L. Cabeceira, J. Represa, P. M. P. So, W. J. R. Hoefer, 2004 FDTD modeling of transient microwave signals in dispersive and lossy bi-isotropic media.IEEE Trans. Microwave. Theory Tech. 52 773 784
47 - A. Grande, I. Barba, A. C. L. Cabeceira, J. Represa, K. Kärkkäinen, A. H. Sihvola, 2005 Two dimensional extension of a novel FDTD technique for modeling dispersive lossy bi-isotropic media using the auxiliary differential equation method. IEEE Microwave Wireless Components Lett. 15 375 377
48 - V. Demir, A. Z. Elsherbeni, E. Arvas, 2005 FDTD formulation for dispersive chiral media using the Z transform method. IEEE Trans. Antennas Propag. 53 3374 3384
49 - V. Nayyerir, M. Soleimani, J. Rashed-Mohassel, M. Dehmollaian, 2011 FDTD modeling of dispersive bianisotropic media using Z-transform method. IEEE Trans. Antennas Propag. 59 2268 2278
50 - J. A. Pereda, A. Grande, O. Gonzalez, A. Vegas, 2006 FDTD modeling of chiral media by using the Mobius transformation technique. IEEE Antennas Wireless Propag. Lett. 5 327 330
51 - J. A. Pereda, A. Vegas, A. Prieto, 2002 FDTD modeling of wave propagation in dispersive media by using the Mobius transformation technique. IEEE Trans. Microwave Theory Tech. 506 1689 1695
52 - N. Bushyager, E. M. Tentzeris, 2006 MRTD (Multi Resolution Time Domain) Method in Electromagnetics (Synthesis Lectures in Computational Electromagnetics). San Francisco, USA: Morgan & Claypool, 116
53 - M. Krumpholz, L. P. B. Katehi, 1996 MRTD: New time-domain schemes based on multiresolution analysis. IEEE Trans. Microwave Theory Tech. 44 555 571
54 - E. M. Tentzeris, A. Cangellaris, L. P. B. Katehi, J. Harvey, 2002 Multiresolution time-domain (MRTD) adaptive schemes using arbitrary resolutions of wavelets, IEEE Trans. Microwave Theory Tech. 50 501 516
55 - M. Fujii, W. J. R. Hoefer, 1998 A three-dimensional Haar-wavelet-based multiresolution analysis similar to the FDTD method-Derivation and application. IEEE Trans. Microwave Theory Tech. 46 2463 2475
56 - S. Grivet-Talocia, 2000 On the accuracy of Haar-based multiresolution time- domain schemeS. IEEE Microwave Guided Wave Lett. 10 397 399
57 - I. Barba, J. Represa, M. Fujii, W. J. R. Hoefer, 2003 Multiresolution model of electromagnetic wave propagation in dispersive materials. International Journal of Numerical Modelling 16 239 247
58 - I. Daubechies, 1992 Ten Lectures on Wavelets.Philadelphia USA SIAM 357
59 - I. Barba, A. Grande, A. C. L. Cabeceira, J. Represa, 2006 A Multiresolution Model of Transient Microwave Signals in Dispersive Chiral Media IEEE Transactions on Antennas & Propagation 54 2808 1812
60 - W. J. R. Hoefer, 1985 The Transmission-Line Matrix method- Theory and Applications IEEE Trans. Microwave Theory Tech. 33 882 893
61 - I. Barba, A. C. L. Cabeceira, J. Represa, M. Panizo, 1996 Modelling dispersive dielectrics for HSCN TLM method Int. J. Numer. Modeling 14 15 30
62 - A. C. L. Cabeceira, I. Barba, A. Grande, J. Represa, 2006 A Time-Domain Modelling for EM wave propagation in Bi-Isotropic media based on the TLM methoD. IEEE Trans. Microwave Theory Tech. 54 2780 2789
63 - S. A. Schelkunoff, 1952 Generalized telegraphist’s equations for waveguides. Bell Syst. Tech. J. 3 784 801
64 - K. Ogusu, 1977 Numerical analysis of rectangular dielectric waveguide and its modifications. IEEE Trans Microwave Theory Tech MTT-25 441 446
65 - Y. Xu, R. G. Bosisio, 1995 An efficient method for study of general bi-anisotropic waveguides. IEEE Trans Microwave Theory Tech MTT-43 873 879
66 - Y. Xu, R. G. Bosisio, 1997 A study on the solutions of chirowaveguides and bianisotropic waveguides with the use of coupled-mode analysis. Microwave Opt Technol Lett 14 308 311
67 - Á. Gómez, A. Vegas, M. A. Solano, 2004 A brief discussion on the different formulations of the coupled mode method in chiral media: Application to the parallel-plate chirowaveguide. Microwave and Optical Technology Letters 42 181 185
68 - Á. Gómez, A. Vegas, M. A. Solano, 2006 On two different formulations of the Coupled Mode Method: Application to 3D rectangular chirowaveguides. AEU-International Journal of Electronics and Communications 60 690 704
69 - Á. Gómez, A. Lakhtakia, A. Vegas, M. A. Solano, 2010 Hybrid technique for analysing metallic waveguides containing isotropic chiral materials. IET Microwaves, Antennas and Propagation 4 305 315
70 - Á. Gómez, I. Barba, A. C. L. Cabeceira, J. Represa, A. Vegas, M. A. Solano, 2008 Application of Schelkunoff’s method for simulating isotropic chiral free propagation: Clarifying some common errors. Journal of Electromagnetic Waves and Applications 22 861 871
71 - M. A. Solano, A. Vegas, A. Prieto, 1994 Numerical analysis of discontinuities in a rectangular waveguide loaded with isotropic or anisotropic obstacles by means of the coupled-mode method and the mode matching method. Int J Numer Model Electron Networks Devices Fields 7 433 452
72 - M. A. Solano, A. Vegas, A. Prieto, 1998 Theoretical and experimental study of two ports structures composed of different size rectangular waveguides partially filled with isotropic dielectrics. Int J Electron 84 521 528
73 - R. F. Harrington, 1993 Field computation by moment methodsIEEE Press Series on Electromagnetic WavesNew JerseyChapter 1
74 - M. A. Solano, A. Vegas, A. Gómez, 2006 Comments on "a comprehensive study of discontinuities in chirowaveguides". IEEE Trans. on Microwave Theory and Techniques 54 1297 1298
75 - Á. Gómez, A. Lakhtakia, J. Margineda, G. J. Molina-Cuberos, M. J. Núñez, Ipiña. J. A. Saiz, A. Vegas, M. A. Solano, 2008 Full-wave hybrid technique for 3-D isotropic-chiral-material discontinuities in rectangular waveguides: Theory and experiment IEEE Trans. on Microwave Theory and Techniques 56 2815 2825


[1] - As it will be shown later these primary fields can be E and  H, D or B and E .