EM modelling of arbitrary shaped dispersive chiral dielectric objects using a 3D leapfrog scheme on unstructured meshes

The standard Yee FDTD algorithm is widely used in computational electromagnetics because of its simplicity and divergence free nature. A generalization of this classical scheme to 3D unstructured co-volume meshes is adopted, based on the use of a Delaunay primal mesh and its high quality Voronoi dual. This circumvents the problem of accuracy losses, which are normally associated with the use of a staircased representation of curved material interfaces in the standard Yee scheme. The procedure has been successfully employed for modelling problems involving both isotropic and anisotropic lossy materials. Here, we consider the novel extension of this approach to allow for the challenging modelling of chiral materials, where the material parameters are frequency dependent. To adequately model the dispersive behaviour, the Z-transform is employed, using second order Padé approximations to maintain the accuracy of the basic scheme. To validate the implementation, the numerical results produced are compared with available analytical solutions. The stability of the chiral algorithm is also studied.


Introduction
Innovative fabrication techniques are currently being employed to create new chiral materials for use in communication technologies. For such materials, the electric permittivity, magnetic permeability and chirality are all frequency dependent parameters. The chirality allows for a negative index of refraction with, simultaneously, a posi- tive relative permittivity and permeability. In addition, chiral materials lead to optical rotatory dispersion (ORD) and circular dichroism (CD). In ORD, the polarization is continuously rotated inside the material, while CD refers to the change of polarization of the incident wave, from linear to elliptical, due to the different absorption coefficients of a right and a left circularly polarized wave. Such materials are being used to create smaller antennas, super lenses [1], polarizers and Radomes or radar cross section (RCS) reducing materials.
The Yee FDTD algorithm is widely used for general computational electromagnetic simulations. The advantages of the algorithm are its explicit nature, simplicity, low computational costs and divergence free nature. A major limitation of the standard implementation is the requirement for a structured Cartesian mesh, which can result in accuracy losses in simulations involving curved material surfaces. This limitation can be alleviated, though not removed, by employing finer meshes, but with a consequent increase in computational costs. To combine geometric flexibility with the efficiency of the Yee FDTD algorithm, a novel hybrid solution method has been developed in which an unstructured mesh procedure is used in the vicinity of material interfaces, while the standard Yee FDTD scheme is used elsewhere [2,3]. The unstructured mesh procedure is essentially an unstruc-tured implementation of the FDTD method, using a primal Delaunay and its Voronoi dual mesh [4]. The approach, which has already been employed for the analysis of scattering and transmission problems involving both isotropic and anisotropic materials, has been shown to require the use of meshes that are up to eight times coarser than those required by standard FDTD methods [5,6].
In this paper, we will consider an extension of this procedure to enable the modelling of problems involving chiral dispersive materials. The approach adopted will employ a Lorentz model for the electric permittivity and magnetic permeability and a Condon model for the chirality [7]. Several methods have been developed for modelling frequency dependent materials within a time domain method, including the auxiliary differential equation method [8], the piecewise recursive convolution method [9] and the Z-transform technique [10]. As these three methods have been found to produce generally comparable results, the Z-transform technique was selected and implemented here. This approach is widely used in signal processing and can be viewed as the discrete version of the Laplace transform. It is readily implemented within the proposed solution algorithm and demonstrates good convergence close to the resonant frequency [11]. Demir et al. [12] appear to have been the first to model chiral materials in 3D with FDTD and using the Z-transform to handle the frequency dependence of the material parameters. In their method, first order approximations were used in the Z-domain and the Z-transform coefficients were derived analytically. Pereda et al. [13] used a different approach, which achieved second order accuracy in the Z-domain by employing a bilinear transformation for calculating Padé approximants [14].
In this paper, we present a hybrid unstructured mesh FDTD approach, in which chiral materials are modelled by employing a generalization of the method introduced by Pereda et al. [13]. The implementation is validated by simulating the transmission of a pulse through a 3D chiral slab in free space and by calculating the RCS of a sphere. The stability of the resulting numerical algorithm, with respect to changes in the mesh and the chirality, is also investigated.

Problem formulation for isotropic dielectric materials
We start with the equations for an isotropic material, in order to review the basic idea of the solution algorithm and to highlight the changes required to enable the modelling of frequency dependent chiral materials. If we assume that that material is non-lossy [5], Ampère's Law and Faraday's Law may be expressed, in scattered field form, as Here, ∂ A denotes the closed curve bounding a surface A, dA is an element of surface area, directed normal to the surface and dl is an element of contour length in the direction of the tangent to the curve. In addition, ε is the electric permittivity and μ is the magnetic permeability. The total electric field, E tot , and the total magnetic field, H tot , are each regarded as being formed as the sum of incident and scattered fields, with the subscripts inc and scat being employed to indicate incident and scattered components respectively. The incident field is assumed to be a plane wave, generated in the far field, which has the form E inc = E 0 cos(ωt − k · r), where E 0 is the polarization vector of the electric field, k is the wave vector, r is the position vector, ω is the angular frequency and t is the time. The incident magnetic field may be determined, using Faraday's Law, as wherek is the unit wave vector and η = √ μ o /ε o is the impedance of free space. For numerical solution, Eqs. (1) and (2) are discretized using two mutually orthogonal co-volume meshes. Although the approach is more general, we describe the implementation here on a primal tetrahedral mesh generated, using a standard Delaunay method [15], and its associated dual Voronoi diagram. Each Voronoi face is a perpendicular bisector of the corresponding Delaunay edge and each Delaunay face is perpendicular to the corresponding Voronoi edge. It is assumed that the Delaunay mesh has N D e edges and that the Voronoi mesh has N V e edges. Each Delaunay edge has an associated perpendicular closed loop of Voronoi edges and, similarly, each Voronoi edge is surrounded by a closed loop of Delaunay edges. The unknowns are located at the midpoints of these edges. At the center of the ith Delauany edge, the unknown is the projection, E scat,i , of the scattered electric field onto the direction of the edge. The unknown at the center of the jth Voronoi edge corresponds to the projection, H scat, j , of the scattered magnetic field onto the direction of the edge. Using the leapfrog time stepping scheme, discretization of Eqs. (1) and (2) results in Fig. 1 The ith Delaunay edge, connecting Delaunay vertices p 1 and p 2 , and the corresponding Voronoi face, formed by the Voronoi edges and Here, Δt denotes the time step, the superscript n denotes an With these staggered equations, the magnetic field is updated over the dual graph at the half time step, using Eq. (5), and the electric field is updated over the primal graph at the full time step, using Eq. (4). At interface boundaries, the material parameters are not constant over the area of integration and, in this case, average values of ε, μ are adopted [6].
The extension of this implementation to any primal polyhedral mesh and its orthogonal dual is direct. In particular, the method reverts to the standard Yee FDTD algorithm when hexahedral meshes are employed.

Mesh generation
Consider the problem of generating a body fitted mesh to discretise the computational domain for a general problem of interest. For the solution of wave propagation problems with the FDTD method, a mesh of uniform size is generally desired. The procedure adopted for constructing a suitable mesh can be decomposed into four main stages, making full use of the ability of the algorithm to be implemented on general polyhedral meshes.
In the first stage, an unstructured triangulation of the surface of the object is produced [16] and this triangulation is then placed inside a hexahedral box. The region inside this box is discretized using a regular Cartesian mesh of cubes, of a uniform edge length. Cubes within a prescribed distance of the object, or lying internal to the object, are removed in a second stage, to create a staircase shaped surface that completely encloses the object [15]. A pyramidal element is added to each exposed internal square face. In the third stage, a standard Delaunay mesh generator is used to create a uniform unstructured tetrahedral mesh that completely covers the unmeshed region. The fourth stage consists of stitching together the hexahedral and tetrahedral meshes [17].
The resulting unstructured mesh is optimized to ensure that both the primal and the dual mesh are of the highest possible quality. The approach adopted is to relax, if necessary, the requirement that a dual edge must be a bisector of the corresponding Delaunay edge. At the same time, the corresponding dual mesh vertex is moved to a point which still ensures orthogonality between the two grids and which lies inside the corresponding primal element. Primal elements with a common circumcenter, and hence, a corresponding Voronoi edge length of zero, will automatically be merged prior to the solution process, creating a mesh consisting of general polyhedral cells [4].

Problem formulation for dispersive chiral materials
The approach outlined above requires some modification before it can be used to model problems involving frequency dependent chiral materials. Incorporation of the Z-transform method enables the development of a time domain formulation that takes into account materials with frequency dependent parameters. It is not possible to employ a pure scattered field formulation and an intermediate total field step is found to be required. In this case, in the frequency domain, Ampère's Law and Faraday's Law may be expressed in the integral form (7) Here, j = √ −1 and D and B represent the electric flux density and the magnetic flux density respectively. The total field formulation removes the terms ε(ω) and μ(ω). The constitutive equations for chiral materials are expressed in the form where c is the speed of light and κ the chiral parameter coupling the electric and magnetic fields. Implementing these relationships within the numerical scheme provides computational challenges, as the electric field projection is stored on the Delaunay mesh while the magnetic field projection is stored on the Voronoi mesh. The Lorentzian model is employed to determine the variation of the permittivity and the permeability with frequency. Here, the subscript ∞ denotes high frequency values, while the subscript refers to the low frequency limit. ω e , ω h refer to the resonance frequency and ξ e , ξ h , to the damping coefficients for the Lorentz model. The value of the chiral parameter is obtained, from the Condon model, as where ω k denotes the resonance frequency, while ξ k corresponds to the damping coefficient for the chiral material. The coupling constant, τ k , defines the magnitude of the chirality. Using these relations, Eqs. (6) and (7) may then be expressed as where

Discrete equations and the Z-transform
To convert Eqs. (13) and (14) from the frequency to the time domain, we use the substitution jω → ∂ ∂t . The staggered leapfrog scheme is then applied to discretize the resulting equations. The bilinear transformation jω = 2 Δt is used, to convert a frequency dependent function to the Zdomain [10]. Equations (15) and (16) can be expressed, with second order accuracy in the Z-domain, by adopting the Padé approximants b In this manner, we obtain the approximation Transferring to the time domain is straightforward, as Z −m acts as a backwards operator in time, e.g. with m = 1, Z −1 J n+1 ee = J n ee . Rearranging Eq. (18) and, for an arbitrary field, F, using the notation Z −m F n = F n−m , we obtain Similar approximations may be written for the terms J n−1/2 κh , J n+1/2 hh and J n+1 κe . Using this notation, discretizing the time dependent form of Eqs. (13) and (14) on our unstructured mesh leads to the update equations is required on the Voronoi edges, but it includes the electric field projections, that are only available on the Delaunay edges. In Fig. 3, we illustrate how we overcome this problem. In the case of a tetrahedron, we have six Delaunay edges and the electric field projection E i , i = 1, . . . , 6 is stored on each edge. Using these projections, the electric field vector E Cell 1 can be reconstructed and linked to the center of the corresponding cell, in this case cell 1. To achieve this, the equation set is as follows Here, E Cell 1 ·ê i corresponds to the projection of the electric field vector E Cell 1 in the direction of the Delaunay edgê e i . The unknown electric field vector at the center of a given cell is approximated by writing Eq. (22) for each of the N Delaunay edges forming the corresponding cell in turn and summing these equations, as depicted in Fig. 3b. In this case, for each cell center, we solve the system with where e i 1 = e i x , e i 2 = e i y and e i 3 = e i z . The resulting field vector is then associated to the corresponding cell center. As a Voronoi edge connects two cells, we need to repeat this procedure for the second cell leading to E Cell 2 . In the next step, the field vectors linked to these cells are averaged and projected to the corresponding Voronoi edge v according to (E Cell1 + E Cell2 )/2 · v, as illustrated in Fig. 3b. With the electric field defined on the Voronoi edges in this manner, we may then compute J n κe,i av . We proceed analogously for the magnetic field, but in this case we have a point from the Delaunay mesh inside a Voronoi cell [6].

Updating scheme
With the total field formulation employed as an intermediate step, to obtain the desired scattered fields, the complete solution algorithm may be described by the steps:

Interface boundary conditions
To satisfy the boundary conditions at interfaces between different materials, a weighted averaging is applied to the Padé approximants. This is similar to the process that has already been successfully employed for simulations involving isotropic and anisotropic media [5]. However, for the Padé approximants, which are required in expressions such as Eq. (18), we only average the coefficients in the numerator.
To illustrate this process, consider two elements, separated by a boundary face, as illustrated in Fig. 4, and consider the calculation of the permeability for the corresponding Voronoi edge. Applying the bilinear transformation to μ ∞ − μ from Eq. (11) leads to (μ ∞ − P(z)), where The magnetic flux is determined as where the subscript i(= 1, 2) refers to the medium 1 or 2. The average magnetic flux at the interface is calculated as where w i = l i /(l 1 + l 2 ) for i = 1, 2.

Rotation of the plane of polarization
To verify the implementation of the scheme, the transmission coefficients for an EM pulse incident upon a chiral slab are computed, along with an evaluation of the variation of the angle of the plane of polarisation of the incident field. The dielectric slab has dimensions 0.3 m, 0.3 m and 0.1 m in the x, y, and z directions respectively. The mesh employed is such that 40 cells are used to discretise a length of 0.1 m. The electric permittivity and magnetic permeability of the slab are described by a Lorentz model and the chirality by a Condon model. The slab is located in free space, as depicted in Fig. 5.
In total the mesh with Δz = 0.0025 m consists of 3, 512, 025 cells corresponding to 10, 676, 124 Delaunay and 10, 417, 154 Voronoi edges. The PML is represented in terms of 1, 220, 300 hexahedral cells. We only consider a real chirality, which induces a rotation of the plane of polarization of the incident field (ORD). In this case, the variation of the angle of polarization, Φ, can be computed as   where κ is the real value of the chirality, ω is the angular frequency of the incident wave and L is the thickness of the slab [18]. To compute this angle, we use a Fourier transform, which allows us to evaluate the variation of the amplitude without time dependence. A point in the Delaunay mesh is selected, such that, with respect to the incident wave, it lies behind the chiral slab. Then, for a wave propagating in z direction, where FT denotes the Fourier transform. Figure 6 shows the numerically computed rotation compared to the theoretical result, for different values of the chirality. For values of κ less than 0.0407, the numerically computed angles of rotation of the plane of polarization remain almost constant over time and are in good agreement with the theoretical values. For higher chiralities, however, the numerical solution shows very strong fluctuations indicating an instability of the scheme.

Stability investigation
Initially, to investigated the effect of the time step size on the results, we performed calculations on the same mesh but with different time steps. From the results of these calculations, presented in Fig. 7a, it is apparent that reducing the time step has little effect on the stability of the scheme.
Next, we investigated the effect of spatial step size on the numerical stability. The slab thickness, L, is expressed as L = nΔz, where n is the number of cells through the slab, in the z direction, and Δz is the edge length of a given cell. A cell normalized rotation of the plane of polarization may then be defined as In the next step, we fix the chirality and create different meshes by just increasing the spatial step Δz and compute the relative error of Φ cell .This leads to an increase in Φ cell for coarser meshes, as represented in Fig. 7b. Finally, for several real chiralities we employ two meshes with Δz = 0.0025 m for the finer mesh and Δz = 0.003 m for the coarser mesh. The results produced after 20 cycles are displayed in Fig. 8a, b. It can be clearly seen that, for the same chirality, the oscillations decrease when the mesh size reduces.
We conclude that the chirality induced instabilities are directly linked to the spatial discretisation and not to the time step. For a stable simulation, the results suggest that the mesh used should ensure that Φ Cell ≤ 0.4 • for each cell.

Reflection/transmission on a chiral slab
Consider the problem of determining the transmission and reflection coefficients for a chiral slab, of infinite length and width, with a finite thickness L, set between two isotropic media. The first medium extends from z ∈] − ∞, −L[, the bi-isotropic chiral slab forms the second medium that extends from z ∈ [−L, 0[ and the third medium extends from z ∈ [0, +∞[. For our case, medium 1 and 3 will be taken to be free space. This problem has an analytical solution [18]. The linearly polarized narrowband pulse propagating in the z direction, is sent towards the slab. The values adopted for the material parameters are ε s = 1.8ε 0 , ε ∞ = 1.8ε 0 , μ s = 1.1μ 0 , μ ∞ = 1.0μ 0 , ω e = ω h = ω k = 2π × 3.5 GHz, ξ e = 0.14, ξ h = 0.12, ξ k = 0.1 and τ k = 1 ps. To compute the frequency dependent reflection, we chose a point p1 in the mesh located in free space just ahead of the slab. For the transmission coefficient, a point p2, slightly behind the slab, is selected. The co-polarized reflection coefficient, R co , the co-polarized transmission coefficient, T co , and the cross polarized transmission coefficient, T cr , are then determined as The results for the transmission and reflection coefficients are displayed in Figs. 10 and 11 respectively. For a simple isotropic material, T cr is equal to zero. Overall, good agreement is obtained between the numerical and the analytical solution. The differences in the results arise because the analytical model assumes a slab with infinite extension in the directions perpendicular to the propagation. In addition, the numerical dispersion inherent in the scattered field formulation leads to a shift in the peaks for the reflection coefficient.

Radar cross section (RCS)
This case involves the computation of the radar cross section (RCS) of a chiral sphere. The RCS corresponds to the amount of scattered power from a target towards a radar and may be readily determined from the numerical solution [5,6]. This problem is chosen as an analytical solution is available [19]. The scattered wave is decomposed into co-polarized and cross-polarized components. We consider a λ/40 sphere and employ a mesh consisting of 4,742,854 Delaunay cells, 1,727,445 Delaunay vertices, 7,604,002 Delaunay edges and 10,544,148 Voronoi edges. The expression λ/n expresses the fact that our mesh consists of a spacing of n points per wavelength. We will use three different wavelengths, 0.8 m, 1.0 m and 1.2 m, corresponding to spheres of electrical length 2.5 m, 2.0 m and 1.66 m respectively. A PML region consisting of 10 layers of cells was found to be sufficient for this case. As this mesh is too dense to display, the form of a coarser mesh, suitable for the analysis of a λ/15 sphere, is shown in Fig. 12. The bi-static RCS of a λ/40 sphere, with varying free space wavelength λ 0 and coupling constants between the electric and magnetic fields, is computed. The wavelength λ inside a dielectric shrinks with respect to free space according to λ = λ 0 / ε r μ r − κ 2 . The material parameters for the electric permittivity and permeability remain the same, because we are interested in the effect of chirality on the RCS. For the simulations, the values ε s = 1.8 ε 0 , ε ∞ = 1.6 ε 0 , μ s = 1.1 μ 0 , μ ∞ = 1.0 μ 0 , ω e = ω h = ω k = 2π × 1 GHz, ξ e = 0.12, ξ h = 0.12 were employed for the material parameters. A time step size Δt = 6.18 ps was used and the wavelength λ 0 of the incident wave took the values 0.8 m, 1.0 m and 1.2 m. The numerical solution was advanced through 20 cycles and the results produced were compared to the analytical solution, using the program developed by Demir et al. [19].
As the RCS distributions will be symmetrical in this case, only the results for angles lying between 180 • and 360 • are shown in Figs. 13, 14, 15 and 16. All the results demonstrate a very good agreement between the numerical and analytical solution for different values of the coupling constants τ k and the wavelength λ 0 . In particular, the very weak signal at high angles, from 300 • to 360 • , is well captured.

Practical example
The final example demonstrates the use of the method in a predictive mode and involves the calculation of the transmission coefficient of a dielectric radome. The radome is in the form of half of an ellipsoidal shell, with a lateral radius of 0.5 m, a length of 1 m and a thickness of 0.05 m. A triangulation of the radome surface is shown in Fig. 17a. The shell is made of a chiral material, with properties characterised by the Lorentz and Condon models of Eqs. (11) and (12), with parameters values ε s = 1.8 ε 0 , ε ∞ = 1.6 ε 0 , μ s = 1.1 μ 0 , μ ∞ = 1.0 μ 0 , ω e = ω h = ω k = 2π × 1 GHz, ξ e = 0.12 and ξ h = 0.12.
A narrowband pulse, as defined in Eq. (32), is used to illuminate the radome, with a central frequency at 1.2 GHz and a width of 1 ns. Relative to an (x, y, z) cartesian coordinate system, the pulse is polarized in the y direction and propagates in the x direction. The mesh employed is designed to have 40 degrees of freedom per wavelength on the shell and a PML represented by 10 layers of hexahedral elements. The shortest distance between the inner boundary of the PML and the radome surface is discretised with 8 cells. This mesh, which is illustrated in Fig. 17b, contains 876,314 cells, 1,673,283 Delaunay edges and 2,076,114 Voronoi edges. The simulation is stopped after 20 periods, a time of 16.6 ns, when all the wave reverberations become insignificant. This required 40 min of CPU time on a single Intel Nehalem processor. Figure 17b also shows computed contours of the scattered field component E scat,y at a time corresponding to 8 ns of propagation of the narrowband pulse.
The transmission coefficient spectrum is computed using a Fourier transform of the time domain solution, as indicated in Eq. (33). For this example, the reference points p1 and p2 are chosen to lie outside and inside the radome respectively. Figure 18 shows the variation of the computed co-polarized reflection and transmission coefficients, R co and T co respectively, and the cross-polarized transmission coefficient, T cr . The results show that the reflection is higher below 3 GHz. At this frequency, the cross-polarized transmission factor shows its maximum and above this frequency the transmission is higher.

Conclusion
In this paper, we have demonstrated how a co-volume solution method, originally developed for applications involving Fig. 17 a View of the surface discretisation of a portion of the radome shell and b contours of E scat,y , and a view of the mesh employed, after propagation of the narrowband pulse for 8 ns Fig. 18 Calculated variation of the co-polarized reflection coefficient, R co , the co-polarized transmission coefficient, T co , and the cross polarized transmission, T cr , for the radome example PECs, isotropic dielectric and anisotropic dielectric materials [5,6], can be further extended to enable the modelling of isotropic frequency dependent materials and chirality. These additional features have been incorporated with the use of the Z-transform method [10,13]. We have investigated the effect of the chirality on the stability of the resulting scheme and demonstrated that stability is enhanced by the use of a finer mesh or by reducing the real part of the chirality. Increasing the damping also has a stabilising effect. As the chirality is obtained by fitting the Condon model to experimental data, varying Re(κ) or I m(κ) is not possible, so that, normally, a stable solution can only be obtained by using finer meshes. With an appropriate mesh, we have successfully simulated the transmission and reflection of an electromagnetic pulse through a chiral slab and computed the RCS of a sphere. The modelling of chiral materials requires the use of a fine mesh, with Φ Cell < 0.4 • for stability.