# Stable 3D FDTD Method for arbitrary Fully Electric and Magnetic Anisotropic Maxwell Equations

## Abstract

We have developed a new fully anisotropic 3D FDTD Maxwell solver for arbitrary electrically and magnetically anisotropic media for piecewise constant electric and magnetic materials that are co-located over the primary computational cells. Two numerical methods were developed that are called non-averaged and averaged methods, respectively. The non-averaged method is first order accurate, while the averaged method is second order accurate for smoothly-varying materials and reduces to first order for discontinuous material distributions. For the standard FDTD field locations with the co-location of the electric and magnetic materials at the primary computational cells, the averaged method require development of the different inversion algorithms of the constitutive relations for the electric and magnetic fields. We provide a mathematically rigorous stability proof followed by extensive numerical testing that includes long-time integration, eigenvalue analysis, tests with extreme, randomly placed material parameters, and various boundary conditions. For accuracy evaluation we have constructed a test case with an explicit analytic solution. Using transformation optics, we have constructed complex, spatially inhomogeneous geometrical object with fully anisotropic materials and a large dynamic range of and , such that a plane wave incident on the object is perfectly reconstructed downstream. In our implementation, the considerable increase in accuracy of the averaged method only increases the computational run time by 20%.

###### keywords:

3D FDTD, Fully-Anisotropic, Stability, Eigenvalue Analysis, Cloaking###### Msc:

[2010] 00-01, 99-00[mycorrespondingauthor]Corresponding author

## 1 Introduction

A rapidly developing field of metamaterials, multiferroic and magnetoelectric materials, applications of transformation optics are at the forefront of the important advancements, new discoveries and practical applications in electomagnetism, (1), (2), (3), (4), (5). Computer simulations of such materials require robust and accurate numerical Maxwell solvers when both electric and magnetic materials are fully anisotropic materials.

The finite-difference time-domain (FDTD) method has a long history of success with isotropic and diagonally-anisotropic materials (6). The key advantages of FDTD are its second-order accuracy, enforcement of continuity of the respective normal and tangential field components across interfaces on Cartesian grids, and its non-dissipative energy preserving nature.

The previous work primarily concentrated on electrically anisotropic materials, (7), (8), (9), (10) and (11). A split-step FDTD method for 3D Maxwellâs equations for the fully anisotropic, but homogeneous media, is presented in Singh *et al.* (12). Our work is inspired by the (10) and (11). In the mathematical proof, we have utilized an idea of ”triplets”, but without introducing new terminology and more importantly, the proof requires changing from the domain of dependence point of view to the domain of influence argument as explained in section 3.2.

In this paper we have developed new fully anisotropic 3D FDTD Maxwell solver for arbitrary electrically and magnetically anisotropic media for piecewise constant electric and magnetic materials that are co-located at the primary computational cells. In particular, we provide a mathematically rigorous proof of the stability of the method. Extensive numerical testing involving extreme and random parameter regimes, and various boundary conditions are utilized to verify the correctness of the implementation and to illustrate the performance of the method. For accuracy evaluation we have constructed a test case with an explicit analytic solution. Using transformation optics, we have constructed complex, spatially inhomogeneous geometrical object with fully anisotropic materials and a large dynamic range of and , such that a plane wave incident on the object is perfectly reconstructed downstream.

The co-location of electric and magnetic materials at the primary computational cells is done for the following reasons: a) since physical fully anisotropic objects do not have electric and magnetic parts shifted by half of the computational cell this approach avoids the necessity for complicated cut-cell algorithms; b) to be suitable for generalizations and applications in finite element context utilizing conventional mesh generators where material discontinuities are allowed across the boundaries between the elements, but not within the individual elements itself. Note, that the co-location of the materials over the primary computational cells breaks the symmetry between the electric and magnetic materials/fields placement. Therefore, two separate algorithms for inverting the constitutive relations were developed, one for the electric field and another one for the magnetic field.

## 2 Fully anisotropic FDTD algorithm

In our implementation of the 3D FDTD method, the relationship between material locations and field components is as follows: for a given computational cell with the origin at the vertex , the -field components are located at the center of their respective low-side faces, while the -field components are located on the low-side of their respective edges, and both the permittivity and permeability matrices are assumed be piecewise constant over the primary computational cell, as shown in figure 1. For example, The and components attributed to computational cell are located at , while and are located at .

If the material distribution is described analytically, either by a smooth or non-smooth function, the material values are replaced by their average values over the primary computational cells.

The standard FDTD discretization of Maxwell’s equations is as follows:

(1) | |||||

(2) | |||||

(3) | |||||

(4) |

where and are the inverses of the symmetric positive definite matrices and , respectively. This method has been shown to be stable for the diagonally-anisotropic, homogeneous and , (13). For fully anisotropic materials, the only changes to the standard FDTD method are due to the presence of the fully anisotropic material constitutive relations, equations (2) and (4).

### 2.1 Non-Averaged Fully Anisotropic Method

The constitutive relations introduce a first-order error as non co-located field components assigned to a single computational cell () are used, and they are on the order of cell size apart from each other.

For example the update for the computational cell is

(5) |

The update for all other field components is shown in B.

### 2.2 Averaged Fully Anisotropic Method

Here we define two separate update methods, one for the and one for the field components. This is due to the difference in the location of the and the fields that are located at the primary computational cell faces and the computational cell edges, respectively, while both materials () are co-located at the primary computational cells as illustrated in figure 1. While being different in terms of number of surrounding computational cells, both update methods follow a similar strategy: 1) first, average s/s inside the computational cells that share a common face (edge) of interest; 2) invert respective constitutive relations in each of the surrounding computational cells; 3) take the arithmetic average of the resulting / field components. In the following subsections, we describe both algorithms for and components. The update algorithms for all of the and field components are given in C.1 and C.2.

#### E-Update

In order to update , consider two computational cells that share common face where is located, as shown in figure 2.

As mentioned above, computing first averaged values of s in each computational cell, followed by inverting the constitutive relations in each computational cell, and finally, taking the arithmetic average of the resulting values, gives the following update formula for ,

(6) | |||||

For implementation efficiency the common terms are grouped together,

(7) | |||||

#### H-Update

To update the edge-centered field component, consider the four computational cells that share common edge where is located with materials , , , and , respectively.

Since there are four surrounding computational cells for each magnetic field component (vs two for each electric field component), we compute first the averaged values of s in each of the four surrounding computational cells. This is followed by inverting the constitutive relations within each computational cell, and finally, taking the arithmetic average of the resulting values, gives the following update formula for ,

(8) | |||||

The simplification of this expression as well as the update algorithms for all of the field components are given in C.2.

## 3 Stability Analysis

In this section, we prove the stability of the two fully anisotropic FDTD methods described in the previous section. We apply a combined approach of Fourier harmonic ansatz (von Neumann analysis) in time and a matrix stability analysis in space. In (14), under an additional assumption that the curl-curl matrix is diagonalizable, it was shown that requiring that the material matrices are symmetric and positive definite (SPD) is sufficient for the stability under the appropriate CFL restriction. Below, we prove that both methods of the previous section have a diagonalizable curl-curl matrix as well as SPD material matrices. Therefore, similarly to the standard FDTD algorithm, both methods are neutrally stable, e.g. is real for each Fourier mode, and thus there is neither growth nor decay in time, under the CFL restriction.

Consider a fully anisotropic Maxwell’s equations in non-dispersive media (),

(9) | |||||

(10) |

where and are electric and magnetic fluxes, and and are the inverses of the permittivity and permeability matrices, respectively,

(11) | |||||

(12) |

Eliminating from the above equations gives . Assuming a harmonic ansatz in time, , results in the following eigenvalue problem

(13) |

To eliminate exponentially growing solutions should be real, and therefore, the eigenvalues of the right-hand side curl-curl operator should be non-negative. For this it is sufficient to require that the material matrices are SPD matrices, (15).

A similar argument is applied to the fully discrete case below. Consider an FDTD algorithm,

(14) | |||||

(15) |

where and are the discrete curl operators applied to the electric and magnetic fields, respectively. Note that on a uniform grid, for PEC or periodic boundary conditions, the following reciprocity relation holds, (16).

To eliminate , we take the forward time difference of equation (14) and substitute using equation 15. This gives,

(16) |

Finally, assuming periodicity in time and applying Fourier harmonic ansatz, , results in the following eigenvalue problem,

(17) |

As for the continuous case, for the stability requirement on to be real, it is suffices that the right-hand side (discrete curl-curl operator) is diagonalizable and has positive eigenvalues, under the CFL time step restriction that is expressed in terms of the spectral radius of the discrete curl-curl operator in (14); (13); (17),

(18) |

where is the spectral radius of the matrix . To show that SPD material matrices imply that the discrete curl-curl operator is diagonalizable and has non-negative eigenvalues (assuming, the reciprocity relation applies), first multiply both sides of equation (17) by . Then using the fact that the material matrices are SPD, , where the square root matrices themselves are SPD matrices (11); (16), equation (17) can be written as

(19) |

Finally, introducing and , gives

(20) |

Since is symmetric positive semi-definite matrix, it is diagonalizable and has non-negative eigenvalues. In the remaining part of this section we prove that both of our methods have SPD material matrices.

### 3.1 SPD of the Material Matrix for the Non-Averaged Method

For the first order non-averaged method, the global material matrix is block diagonal with each block consisting of physical, , SPD electric permittivity or magnetic permeability matrices with blocks ordered according to the computational cell ordering , where the indices are traversing along , then , and finally along the directions, respectively, as illustrated in equation 26 below,

(21) |

where , , and denote the electric field components and the inverse permittivity matrix associated with the computational cell . For example,

(22) |

The indices label the last computational cell.

### 3.2 SPD of the Material Matrix for the Averaged Method

Consider constitutive relation for electric field, equation (11),

(23) |

where and are column vectors containing all field components of the 3D computational domain written in some fixed ”standard” order. For example, we may chose the standard ordering that corresponds to the computational cell by choosing the electric field components located on the low-side faces of each computational cell as described in section 2, see figure 1 and equation (26).

Now, consider the local vertex labeling for each computational cell as indicated in figure 4. Note, that while we chose to associate with each computational cell the fields located on the lower-side faces, corresponding to vertex 1 in figure 4, we could have chosen any vertex and respective faces to represent the electric field of the computational cell. We enumerate these eight possibilities to represent the global electric field as , and . These vectors will contain all the field components without repetition. When the standard ordering of the global electric field for each computational cell is based on vertex 1 , the global vector in equation (23) is the same as .

Equation (6) (and similar equations for the other field components) may be interpreted either from the ”domain of dependence” point of view, where we focus on the contributions to each field component from the surrounding displacement vectors that share common face, or from the ”domain of influence” argument, where we collect the electric field components that are influenced by a single displacement vector corresponding to a particular vertex within each computational cell. In fact, for any vertex of the computational cell, the corresponding displacement vector will only contribute to the electric field components that are located at the same faces.

Consider an arbitrary computational cell and eight local displacement vectors corresponding to each vertex. The components of these vectors are located on the cell faces intersecting at a particular vertex as shown in figure 4 for vertex 7.

The displacement vector corresponding to vertex 7 contributes to the components of the global electric field that are located at the same faces,

(24) |

Combining these relations for each computational cell into global vectors and , gives

(25) |

where the corresponding material matrix is block diagonal SPD and is ordered according to the computational cell numbering, but the global arrays consisting of the vector components for each vertex , and , are ordered differently from the standard ordering. For example, for computational cell the standard ordering electric field components attributed to and respectively are,

(26) |

Repeating this argument for each vertex labeling of the global electric field we obtain the relations

(27) |

where the material matrix is ordered according to the computational cell numbering which is independent of the vertex chosen to label the global electric field. Since the ordering of the vector components for each and are distinct from the standard ordering, they have to be permuted into a standard order denoted by and , respectively, before they can be summed up to obtain the global field.

This can be accomplished by multiplying vectors and by corresponding orthogonal permutation matrix , distinct for each . Let and , where is a global displacement vector and represents the portion of the global electric field according to the formula (6) that can be written in vector form as

(28) |

Substituting and into equation 27 gives

(29) |

or

(30) |

where . Note, that matrices are SPD matrices for each , since they are orthogonally similar to an SPD matrix . With all of the field components now being in standard order, the electric field update in index form in equation 6 can be written in a vector form as

(31) |

where the global material matrix is an SPD matrix since it is equal to the sum of the SPD matrices. Using equation (8), a similar argument can be applied to show that the magnetic permeability matrix is a SPD matrix as well.

## 4 Numerical Validation of Stability

In this section we validate an implementation of the fully anisotropic methods using two numerical tests. First, we perform long-time integration for a range of random high-contrast material interfaces. Second, we utilize an eigenvalue analysis of the global update matrix on a number of sample domains, showing that their eigenvalue spectra remain on the unit circle. We have not tested the diagonalizability of the global update matrix due to unstable/inconclusive nature of the Jordan form computation in finite precision arithmetic, (18).

### 4.1 Long-Time Integration Test

To perform the long-time integration test, we ran simulations on a 24x24x24 cubic domain with periodic boundary conditions with lattice constant of for 60 million iterations. The simulated geometries ranged from off-centered cubes and spheres to a random distribution, consisting of high-contrast, fully anisotropic materials ( in equation (32)) embedded in a vacuum,

(32) |

In the random distribution each computational cell was randomly assigned as either fully anisotropic in , fully anisotropic in , fully anisotropic in both, or vacuum.

We excited the domain with a broad-band, 2fs Gaussian pulse, exciting all possible modes of the empty cavity, from the lowest mode dictated by the lattice constant to the highest mode dictated by Nyquist. The computed solutions remained bounded with no evidence of exponentially growing oscillations up to the 60 millionth iteration, corresponding to approximately 2 million periods of the lowest vacuum mode.

### 4.2 Eigenvalue Analysis

Here, we demonstrate that the update matrix has eigenvalues lying on the unit circle for an arbitrary distribution of anisotropic materials.

Both update methods being linear operators may be represented as

(33) |

where is the update matrix, is a vector containing the and fields at each spatial location in the computational domain, and the superscript denotes the iteration number in time. We numerically construct the matrix by applying a single iteration of the field update to each of the canonical basis elements of the vector , e.g. . The corresponding eigenvalue equation for is

(34) |

where is an eigenvector of and is its corresponding eigenvalue. For a domain the vector contains 6 field components from and resulting in a matrix with elements, producing eigenvalues. Computing the eigenvalues of the large matrix puts a practical upper limit on the size of the computational domains on which this stability analysis can be performed.

The update matrix, , encodes information regarding the entire computational domain configuration, including material properties and layout, CFL, and boundary conditions; thus the eigenvalue stability analysis is only relevant to the computational domains tested. We numerically construct update matrices and compute their eigenvalues for computational domains with a range of materials and material distributions. In the following, we present the results for two representative material layouts: an isolated asymmetrically located sphere and a random distribution of high-contrast fully anisotropic materials. Using high-contrast material interfaces in test domains is particularly important due to their tendency to introduce instability (11). Smooth material distributions, such as those used as metamaterial cloaks, were also tested with the same, stable results. A sphere is chosen because, when discretized, it presents a complex staircased material interface with a wide range of nearest neighbor material arrangements. Similarly, a domain composed of a random distribution of small, high-index structures presents a large material interface surface, maximizing the possibility of seeding exponentially growing solutions. For the analysis, grids are used, resulting in update matrices, and the CFL used was . The eigenvalues corresponding to each of the update matrices are calculated and plotted along with the unit circle in figure 9. In each case we note that each of the eigenvalues of lie on the unit circle within machine precision of the eigenvalue solver.

## 5 Accuracy and Convergence

Having proved the stability and validated the implementation numerically of both fully anisotropic methods under consideration, we now turn our attention to their accuracy properties. Due to the limited availability of analytical test solutions for fully anisotropic materials, we construct a test inspired by the method of manufactured solutions, the key idea of which to construct an exact solution, without being concerned about its physical realism (19).
Using the transformation optics approach, we construct a complex, spatially inhomogeneous geometrical object using fully anisotropic materials such that a plane-wave incident on the object is perfectly reconstructed downstream. Such an object was described in the field of transformation optics by Pendry *et al.* (20) where the objects are referred to as electromagnetic cloaks.

As the test was designed for numerical investigation of the accuracy properties, we are not concerned with any physical cloaking properties of the objects other than their ability to reconstruct the plane-wave downstream. Such cloaking objects have a number of desirable properties: they are composed of fully anisotropic materials, contain materials with a large dynamic range of and , and exhibit a simple plane wave solution against which we can compute relative error.

To characterize the field error due to discretization, we illuminate a cloak from one side with a plane wave, and quantify the field error by comparing the computed solution to the analytically expected plane-wave over a rectangular region on the downstream side of the cloak, as shown in figure 10.

Using the transformation optics approach, we define two versions of a cloak: one with a smoothly-varying material distribution function with a continuous first derivative, and one for the discontinuous material distributions.

Details of the spatial transformations used to generate the two cloaks are shown in A. The smooth cloak is used to characterize the accuracy of smoothly-varying material distributions, and the non-smooth cloak is used to characterize the accuracy of non-smoothly-varying or high-contrast material distributions. Both of these cloaks are discretized onto the computational grid. Matrix components and along a 1D cut through the computational grid are plotted in figure 11.

Each cloak is simulated in a 3D computational domain of size . The domain is discretized at various resolutions, quantified as points per wavelength (ppw), with periodic boundary conditions employed in and , and uniaxial perfectly matched layer (UPML) boundary in , the propagation direction. A polynomial of order 3 () and 10 cells of UPML () were used in computing the grading of the UPML layers. The and are fixed to 1 and respectively, where is the cell size in the UPML region. A spherical cloaking region with an approximate radius of influence of is located at the center of the domain. Each cell in the spherical region is assigned a particular tensor value of and dictated by either a smooth or non-smooth function defining the cloak. A plane wave (), propagating in the negative z-direction, is generated by total-field scatter-field source, located at the edge of the domain.

To ensure the results were not sensitive to the domain boundaries and/or PML thickness, several configurations were considered when determining a suitable computational domain. The particular simulation configuration presented here has a number of advantages that minimizes the error introduced by PML boundaries. At the source end we use a uni-directional source which limits the reflections from the source-side PML layer. At the transmitted end, there is near-normal incidence of the rays impinging on the PML layer. Also, periodic boundary conditions are used on all faces normal to the propagation direction which introduce no error into the domain.

Electric fields are recorded from each computational cell in a sampling box (white region in figure 10) on the transmitted side of the cloaked region. The sampling box is chosen to encloses half of a wavelength. We then compare the electric fields in the sampling region to the analytical solution in two stages. First, knowing that the isotropic FDTD method converges with second-order accuracy to the analytical plane wave solution, we record a numerical simulation of a plane wave in a vacuum using the isotropic FDTD algorithm. Second, we simulate a cloak with the fully anisotropic FDTD algorithm and compute the relative error to the corresponding isotropic FDTD plane-wave simulation on the same grid.

We compute the scaled norm of the relative error using,

(35) |

where is the number of points in the sampling box, indexes each point in the box, is the reference vacuum electric field at point , is the electric field of the test run containing the cloak at point . Thus, we can determine the accuracy of the fully anisotropic FDTD to the analytical plane wave up to second-order, the order of accuracy of the isotropic FDTD to the analytical plane wave solution.

We compute the scaled norm of the relative error over the sampling region and plot the relative error versus the grid resolution for both the non-averaged and averaged methods in figure 12. The non-averaged method exhibits first-order accuracy for each cloak, as is expected due to the first-order error in material assignment on the discrete grid. With the averaged method, we have demonstrated second-order accuracy for the smoothly-varying cloak. This is consistent to what was seen in the fully anisotropic case for continuously varying material values(11). In our implementation, the considerable increase in accuracy and higher-order accuracy of the averaged method only increases the computational run time by 20%.

## 6 Summary

We have developed a new FDTD Maxwell solver for fully anisotropic electric and magnetic materials assuming that the material is piecewise constant and are co-located over the primary computational cells. We have provided a mathematically rigorous proof of the stability of the two algorithms under consideration, and illustrated it by extensive numerical tests that include long-time integration, an eigenvalue analysis for extreme and random material parameter values, and various boundary conditions; accuracy evaluation for a test case designed to have an explicit analytic solution.

We found that the non-averaged method is first order accurate, while the averaged method is second order accurate for materials with continuous first derivative and it reduces to first order for the discontinuous material distributions. In our implementation, the considerable increase in accuracy of the averaged method only increases the computational run time by 20%.

## 7 Acknowledgments

This material is based on upon work supported by the Air Force Office of Scientific Research under award numbers FA9550-16-1-0088 and FA9550-16-1-0199.

## References

## Appendix A Electromagnetic Cloak

To construct the cloaking objects used in section 5, we follow the general principles were outlined in (20). The cloaks are created by specifying a spatial transformation to 3D uniform grids that pull ray-paths out a desired region (located at the center of the domain), calculating the Jacobian matrix for the transformation, and combining the Jacobian into the material parameters ( and ).

For the smooth cloak, the smooth spatial transformation is defined by

(36) |

where is the uniform spatial coordinate and is the transformed coordinate.

For the non-smooth cloak, the non-smooth piecewise-linear transformation outlined in (21),

(37) |

This is used over the canonical piecewise-linear cloak (20) due to the infinite material parameters needed for the singularity of an ideal cloak by mapping the finite volume to another finite volume rather than a single point (22).

We then obtain the material parameters for the cloak from the above transformations using the formulation outlined in (21) as follows:

(38) | |||||

(39) |

where is the Jacobian matrix, which we compute numerically.

The parameters for the smooth cloak in the accuracy test runs were , and in equation (36). The parameters for the non-smooth cloak used in the accuracy test runs were , and in equation 37. Note we employ cloaks that are non-ideal in that they allow a percentage of light into the cloaked region. This is required to keep the materials finite (22), though under plane wave illumination, the downstream solution is the desired plane wave.

## Appendix B Non-Averaged Constitutive Relations Update

Each material and field component has a subscript that indicates the position of FDTD computational cell in the overall rectangular domain. The relative positions of each material and field component are shown in figure 1. Furthermore, each material also has a subscript denoting the specific matrix element. The and field component update are defined as follows.

### b.1 E-field update

:

:

:

### b.2 H-field update

:

:

:

## Appendix C Averaged Constitutive Relations Update

In this section we explicitly describe field and material averaging methods for fully anisotropic and .

### c.1 E-field update

:

: