Stable 3D FDTD Method for arbitrary Fully Electric and Magnetic Anisotropic Maxwell Equations
J. Nehls C. Dineen J. Liu C. Poole M. Brio J. V. Moloney
math.NA 65M06




  
  



ACMS,COS]J. Nehls mycorrespondingauthor
[mycorrespondingauthor]Corresponding author
nehls@email.arizona.edu
ACMS,DoM]C. Dineen
NC]J. Liu
WI]C. Poole
ACMS,DoM]M. Brio
ACMS,COS,DoM]J.V. Moloney

[ACMS]Arizona Center for Mathematical Sciences, University of Arizona, Tucson, Arizona 85721
[DoM]Department of Mathematics, University of Arizona, Tucson, Arizona 85721 
[COS]College of Optical Sciences, University of Arizona, Tucson, Arizona 85721 
[NC]Department of Mathematical Sciences, Delaware State University, Dover, DE 19901 
[WI]Department of Physics, Lawrence University, Appleton, WI 54911













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%.





3D FDTD Fully-AnisotropicStability Eigenvalue Analysis Cloaking
[2010] 00-0199-00








§ 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, <cit.>, <cit.>, <cit.>, <cit.>, <cit.>. 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 <cit.>. 
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, <cit.>, <cit.>, <cit.>,   <cit.>  and <cit.>. A split-step FDTD method for 3D Maxwell’s equations for the fully anisotropic, but homogeneous media, is presented in Singh et al. <cit.>. Our work is inspired by the  <cit.>  and <cit.>. In the mathematical proof, we have utilized their idea of "triplets", but  without introducing new terminology and more importantly, the proof switches from the domain of dependence point of view to the domain of influence argument as explained in section <ref>.  We have generalizes the results of these two papers on FDTD that deal with anisotropic permittivity only and semi-discrete stability considerations to fully electric and magnetic anisotropic materials and provided a rigorous mathematical proof of stability in a fully discrete case. We have also shown that stability of a particular algorithm may depend on the material placement. In particular, "stable" algorithm of <cit.> that was determined to be unstable and modified to be stable in  <cit.>, are all the same algorithm in our material placement when applied to the electric field, and is proven to be stable in fully discrete case. The magnetic field treatment requires a different algorithm to be provable stable. We have also added a missing component in the proof in  <cit.> that the global material matrix is an SPD matrix by including necessary permutations applied to global electric field and electric displacement vectors, as well as respective material matrices. 

 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. 

The paper is organized as follows: in section <ref>, we present the fully anisotropic FDTD algorithm. 
In section <ref>, we analyze the stability of the anisotropic FDTD algorithm, followed by the 
numerical stability tests in section <ref>. The accuracy study is presented in section <ref>. 





 



§ 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  (i,j,k), the D/E-field components are located at the center of their respective low-side faces, while the B/H-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 <ref>. For example, The D_x(i,j,k) and E_x(i,j,k)
 components attributed to computational cell (i,j,k) are located at (x_i, y_j+1/2, z_k+1/2), while  B_x(i,j,k) and H_x(i,j,k)   are located at (x_i+1/2, y_j, z_k). 
  


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:


    B^n+1/2   =   B^n-1/2 - Δ t  ∇× (E^n), 
    H^n+1/2   =   ζ B^n+1/2,  
    D^n+1   =   D^n +  Δ t  ∇× (H^n+1/2), 
    E^n+1   =   ξ D^n+1,
 
where ζ and ξ are the inverses of the   3× 3 symmetric positive definite matrices μ and ϵ, respectively. This method has been shown to be stable for the diagonally-anisotropic, homogeneous  ϵ and μ, <cit.>.
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 (<ref>) and (<ref>).   
 


 §.§ 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 (i,j,k) are used, and they are on the order of cell size apart from each other.

For example the E_z update for the computational cell (i, j, k) is 


    E_z (i,j,k) =  [ [ ξ_zx ξ_zy ξ_zz ]]_ijk[ [ D_x; D_y; D_z ]]_ijk.


The  update for all other field components is shown  in <ref>.




 §.§ Averaged Fully Anisotropic Method

Here we define two separate update methods, one for the E and one for the H field components.  
This is due to the difference in the location of the E and the H fields that are located  at the primary computational cell faces and the computational cell edges, respectively, while both  materials (ϵ and μ)  are co-located at the primary computational cells as illustrated in figure <ref>.
While being different in terms of number of surrounding computational cells, both update methods follow a similar strategy: 1) first, average Ds/Bs 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 E/H field components. 
In the following subsections, we describe both algorithms for E_z and H_z components. The update algorithms for all of the E and H field components are given  in <ref> and <ref>.



  §.§.§ E-Update

In order to update E_z(i,j,k), consider two computational cells that share common face where E_z(i,j,k) is located, as shown in figure <ref>.




As mentioned above, computing first averaged values of Ds in each computational cell, followed by inverting the constitutive relations in each computational cell, and finally, taking the arithmetic average of the resulting E_z values, gives the following update formula for E_z(i,j,k), 




    E_z(i,j, k)   =   1/8[ [ ξ_zx ξ_zy ξ_zz ]]_i,j,k(     [ [ D_x(i,j,k); D_y(i,j,k); D_z(i,j,k) ]]   +
      [ [ D_x(i+1,j,k);   D_y(i,j,k);   D_z(i,j,k) ]] +
      [ [   D_x(i,j,k); D_y(i,j+1,k);   D_z(i,j,k) ]]  +
      [ [ D_x(i+1,j,k); D_y(i,j+1,k);  D_z (i,j,k) ]]   ) 
       +   1/8[ [ ξ_zx ξ_zy ξ_zz ]]_ij,k-1(     [ [ D_x(i,j,k-1); D_y(i,j,k-1);   D_z(i,j,k) ]]   +
      [ [ D_x(i+1, j,k-1);    D_y(i,j,k-1);      D_z(i,j,k) ]] +
      [ [   D_x(i,j,k-1); D_y(i,j+1,k-1);     D_z(i,j,k) ]]  +
      [ [ D_x(i+1,j,k-1); D_y(i,j+1,k-1);     D_z(i,j,k) ]]   ).


For implementation efficiency the common terms are grouped together,  

    E_z(i,j,k)   =   1/4( ξ_zx(i,j,k) (D_x(i,j,k) + D_x(i+1,j,k)) +        ξ_zx(i,j,k-1) (D_x(i,j,k-1) + D_x(i+1,k,k-1)))
       +   1/4( ξ_zy(i,j,k) (D_y(i,j,k) + D_y(i,j+1,k)) +    ξ_zy(i,j,k-1) (D_y(i,j,k-1) + D_y(i,j+1,k-1)))
       +   1/2(ξ_zz(i,j,k) + ξ_zz(i,j,k-1)) D_z(i,j,k).






  §.§.§ H-Update

To update the edge-centered H_z(i,j,k) field component, consider the four computational cells that share common edge where H_z(i,j,k) is located with materials ζ_(i,j,k), ζ_(i-1,j,k),  ζ_(i,j-1,k), and  ζ_(i-1,j-1,k), 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 Bs 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 H_z values,  gives the following update formula for H_z(i,j,k), 





    H_z(i,j, k)   =   1/8[ [ ζ_zx ζ_zy ζ_zz ]]_i,j,k(     [ [ B_x(i,j,k); B_y(i,j,k); B_z(i,j,k) ]]   +
      [ [ B_x(i, j,k+1);  B_y(i,j,k+1);    B_z(i,j,k) ]] ) 
       +   1/8[ [ ζ_zx ζ_zy ζ_zz ]]_i-1,j,k(     [ [ B_x(i-1,j,k);   B_y(i,j,k);   B_z(i,j,k) ]]  +
                          [ [ B_x(i-1,j,k+1);   B_y(i,j,k+1);     B_z(i,j,k) ]] ) 
       +   1/8[ [ ζ_zx ζ_zy ζ_zz ]]_i,j-1,k(     [ [    B_x(i,jk); B_y(i,j-1,k);    B_z(i,jk) ]]   +
      [ [ B_x(i, j,k +1); B_y(i,j-1,k+1);      B_z(i,jk) ]] ) 
       +   1/8[ [ ζ_zx ζ_zy ζ_zz ]]_i-1,j-1,k( [ [ B_x(i-1,j,k); B_y(i,j-1,k);   B_z(i,j,k) ]]  +
      [ [ B_x(i-1,j,k+1); B_y(i,j-1,k+1);      B_z(i,jk) ]]  ).


The simplification of this expression as well as the update algorithms for all of the H field components are given   in <ref>.




§ 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 <cit.>, 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 (J=0), 


    D_t    =   ∇× (M_ζ B),  
     
    B_t    =    - ∇×  (M_ξ D),


where B and D are electric and magnetic fluxes, and M_ζ and M_ξ  are the inverses of the permittivity and permeability matrices, respectively,


    E    =    M_ξ D,   
     
    H    =    M_ζ B.


Eliminating B from the above equations gives 
D_tt = -∇× (M_ζ∇× (M_ξ D)). Assuming a harmonic ansatz in time,
D(x,y,z,t) = e^ i ω t D̃(x,y,z), results in the following eigenvalue problem 


    ω^2 D̃ = ∇× (M_ζ∇× (M_ξD̃)).


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, <cit.>.  

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


    D^n-D^n-1/Δ t   =    C_h M_ζ B^n -1/2, 
    B^n+1/2-B^n-1/2/Δ t   =    -  C_e M_ξ D^n,


where C_h and C_e 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, C_h=C_e ^T <cit.>.    

To eliminate B, we take the forward time difference of equation (<ref>) and substitute  (B^n+1/2-B^n-1/2 )/Δ t using equation <ref>. This gives, 


    (D^n+1- 2 D_n +D^n-1 )/Δ t^2   = - C_h M_ζ C_e M_ξ D^n.


Finally, assuming periodicity in time and applying Fourier harmonic ansatz, D^n=D̃ e^ i ω n Δ t, results in the following eigenvalue problem, 


    sin^2 (ωΔ t /2)   D̃ =  Δ t^2/4   C_h M_ζ C_e M_ξD̃.


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 <cit.>,


    Δ t ≤2/√(ρ (C_h M_ζ C_e M_ξ) ),


where ρ is the spectral radius of the matrix C_h M_ζ C_e M_ξ.

To show that SPD material matrices imply that the discrete curl-curl operator is diagonalizable and has non-negative eigenvalues (assuming,  the reciprocity relation C_e = C, C_h = C^T applies), first multiply both sides of equation (<ref>) by   
M_ξ^1/2. Then using the fact that the material matrices are SPD,  M_ξ=M_ξ^1/2 M_ξ ^1/2 , where the square root matrices themselves are SPD matrices <cit.>, equation (<ref>) can be written as 


    sin^2 (ωΔ t/2)  M_ξ^1/2D̃ = Δ t^2/4   (M_ξ^1/2)^ T  C^T (M_ζ^1/2)^T    M_ζ^1/2 C M_ξ^1/2 M_ξ ^1/2D̃.


Finally, introducing A = M_ζ^1/2 C M_ξ^1/2 and 
D̂ =  M_ξ^1/2D̃,  gives 
 

    sin^2 (ωΔ t/2) D̂  = (A^T A) D̂.


Since A^T A 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.



 §.§ 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, 3× 3, SPD electric permittivity or magnetic permeability matrices with blocks ordered according to the computational cell ordering  (i,j,k), where the indices are traversing along x, then y, and finally along the z directions, respectively, as illustrated in equation <ref> below, 



















    [ [      E_(i,j,k);    E_(i+1,j,k);              ⋮; E_(nx, ny, nz) ]]
        = [    ξ_(i,j,k)            0            …            0;            0  ξ_(i+1,j,k)            …            0;            ⋮            ⋮            ⋱            ⋮;            0            0            … ξ_(nx,ny,nz) ][ [      D_(i,j,k);    D_(i+1,j,k);              ⋮; D_(nx, ny, nz) ]]

where  E_(i,j,k), D_(ijk), and ξ_(i,j,k) denote the electric field components and the inverse permittivity matrix associated with the computational cell (i,j,k).  For example, 


    D_(ijk)=[ [ D_x; D_y; D_z ]]_ijk,       ξ_(ijk)= [ ξ_xx ξ_xy ξ_xz; ξ_yx ξ_yy ξ_yz; ξ_zx ξ_zy ξ_zz ]_i,j,k.

The indices (nx, ny, nz) label the last computational cell.





 §.§ SPD of the Material Matrix for the Averaged Method



Consider constitutive relation for electric field,  equation (<ref>),

    E = M_ξ  D,
  
where E and D 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 (i,j,k) by choosing the electric field components located on  the low-side faces of each computational cell as described in section <ref>, see figure <ref> and equation (<ref>).  

Now, consider the local vertex labeling for each computational cell as indicated in figure  <ref>. 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  <ref>, 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  E^(m), and D^(m) m=1,2,… , 8. 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 E in equation (<ref>) is the same as  E^(1). 

Equation (<ref>) (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 <ref> 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,  


    [ [ E_x; E_y; E_z ]]_i+1,j+1,k+1⟸[ ξ_xx ξ_xy ξ_xz; ξ_yx ξ_yy ξ_yz; ξ_zx ξ_zy ξ_zz ]_i,j,k[ [ D_x; D_y; D_z ]]_i+1,j+1,k+1.


Combining these relations for each computational cell into global vectors  E^(7) and  D^(7), gives 



















    E^(7)= [ [    Ẽ_(i,j,k)^(7);  Ẽ_(i+1,j,k)^(7);                ⋮; Ẽ_(nx,ny,nz)^(7) ]]
        = [    ξ_(i,j,k)            0            …            0;            0  ξ_(i+1,j,k)            …            0;            ⋮            ⋮            ⋱            ⋮;            0            0            … ξ_(nx,ny,nz) ][ [    D̃_(i,j,k)^(7);  D̃_(i+1,j,k)^(7);                 ⋮; D̃_(nx,ny,nz)^(7) ]],


where the corresponding material matrix M_ξ 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 m,  D^(m) and E^(m),  are ordered differently from the standard ordering. For example, for computational cell (i,j,k) the standard ordering electric field components attributed to E^(1) and E^(7) respectively are, 


    Ẽ_(i,j,k)^(1) = [ [ E_x(i,j,k); E_y(i,j,k); E_z(i,j,k) ]],      Ẽ_(i,j,k)^(7) = [ [ E_x(i+1,j+1,k+1); E_y(i+1,j+1,k+1); E_z(i+1,j+1,k+1) ]].



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

 
    E^(m)=M_ξ D^(m) ,  m=1,2,… , 8,

where the material matrix M_ξ 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 D^(m) and E^(m)  are distinct from the standard ordering, they  have to be permuted into a standard order denoted by  D and E^(m), respectively,  before they can be summed up to obtain the global E field.

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

    E =1/8∑_m=1^8 E^(m),



Substituting  D^(m)= P^(m)  D and E^(m)= P^(m) E^(m) into equation <ref> gives 


    P^(m) E^(m) = M_ξ P^(m)  D,
  
or   
 
    E^(m)=M_ξ,m  D.
  
where M_ξ,m =  (P^(m))^T  M_ξ  P^(m). Note, that   M_ξ,m matrices are  SPD matrices for each m, since they are orthogonally similar to an SPD matrix M_ξ. 
With all of the field components now being in standard order,  the electric field update in index form in equation <ref> can be written in a vector form as


    E =1/8∑_m=1^8 E^(m)=(1/8∑_m=1^8 M_ξ,m)  D =M_ξ  D,

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




§ 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, <cit.>.




 §.§ 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 a = 4.8 μ m 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 (Γ = 1, 50, 100, 144 in equation (<ref>)) embedded in a vacuum, 

    ϵ(Γ) = Γ[       10.225       -0.825 -0.55 √(3/2);       -0.825       10.225  0.55 √(3/2);  -0.55√(3/2)  0.55 √(3/2)         9.95 ]ϵ_0 , μ (Γ) = Γ[        3.75        0.75 -0.5 √(3/2);        0.75        3.75 -0.5 √(3/2);  -0.5√(3/2) -0.5 √(3/2)         3.5 ]μ_0.
 
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.



 §.§ 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


    U^n+1  =  AU^n,


where A is the update matrix, U is a vector containing the D and B fields at each spatial location in the computational domain, and the superscript n denotes the iteration number in time. 
We numerically construct the matrix A by applying a single iteration of the field update to each of the canonical basis elements of the vector U, e.g. U=(1,0,… , 0),  (0,1,… ,0)  (0,0,…, , 1). 
The corresponding eigenvalue equation for A is 

    Aν_i  =  λ_i ν_i ,

where ν_i is an eigenvector of A and λ_i is its corresponding eigenvalue.
For a N× N× N domain the vector contains 6 field components from D and B resulting in a matrix A with (6N^3)^2 elements, producing 6N^3 eigenvalues.
Computing the eigenvalues of the large matrix A puts a practical upper limit on the size of the computational domains on which this stability analysis can be performed. 

The update matrix, A, 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 <cit.>.
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, 12×12×12 grids are used, resulting in 10368×10368 update matrices, and the CFL used was 0.4.
The eigenvalues corresponding to each of the update matrices are calculated and plotted along with the unit circle in figure <ref>.
In each case we note that each of the 10368 eigenvalues of A lie on the unit circle within machine precision of the eigenvalue solver.






§ 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 <cit.>.
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. <cit.> 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 <ref>.







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 <ref>.
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 ϵ_xx and ϵ_xy  along a 1D cut through the computational grid   are plotted in figure <ref>.


Each cloak is simulated in a 3D computational domain of size 500nm×500nm×1000nm. 
The domain is discretized at various resolutions, quantified as points per wavelength (ppw),  with periodic boundary conditions employed in x and y, and uniaxial perfectly matched layer (UPML) boundary in z, the propagation direction. 
A polynomial of order 3 (m_pml =3) and 10  cells of UPML (n_pml = 10) were used in computing the grading of the UPML layers.
 The  σ_max and κ_max are fixed to 1 and  8 × (m_pml + 1)/(n_pmlΔ) respectively, where Δ is the cell size in the UPML region. 
A spherical cloaking region with an approximate radius of influence of 200nm 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 (λ = 200nm), 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 <ref>) 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 L_1 norm of the relative error using,

    Relative Error = 1/N∑_j | E_j - Ê_j|/|E_j|,

where N is the number of points in the sampling box, j indexes each point in the box, E_j is the reference vacuum electric field at point j,  Ê_j is the electric field of the test run containing the cloak at point j.
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 L_1 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 <ref>. 
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<cit.>.
In our implementation, the considerable increase in accuracy and higher-order accuracy of the averaged method only increases the computational run time by 20%.






§ 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%.


 



§ 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









§ ELECTROMAGNETIC CLOAK
 

To construct the cloaking objects used in section <ref>, we follow the general principles were outlined in <cit.>. 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

    r^' = {1-depth×exp(-( r/σ)^n ) } r,

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

For the non-smooth cloak, the non-smooth piecewise-linear transformation outlined in <cit.>,

    r = f(r') = R_1 r'/R'_1,      if  r' <R'_1, 
     R_1 - R_2/R'_1- R_2(r' - R'_1) + R_1,      if  R'_1 ≤ r' ≤ R_2,
    
       r'     otherwise.

This is used over the canonical piecewise-linear cloak <cit.> due to the infinite material parameters needed for the singularity of an ideal cloak by mapping the finite volume R_1 to another finite volume R^'_1 rather than a single point <cit.>.


We then obtain the material parameters for the cloak from the above transformations using the formulation outlined in <cit.> as follows:

    ϵ'    =    |Λ| Λ^-1ϵΛ^-T
    μ'         =    |Λ| Λ^-1μΛ^-T,

where  Λ is the Jacobian matrix, which we compute numerically. 

The parameters for the  smooth cloak in the accuracy test runs were n =3, depth = 0.8, and σ = 80nm in equation (<ref>). The parameters for the non-smooth cloak used in the accuracy test runs were R_1 = 8nm, R_2 = 130nm and R^'_1 = 40nm in equation <ref>. 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 <cit.>, though under plane wave illumination, the downstream solution is the desired plane wave.





§ NON-AVERAGED CONSTITUTIVE RELATIONS UPDATE
 

 
Each material and field component has a subscript i,j,k 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 <ref>.
Furthermore, each material also has a {xx,xy,xz,...} subscript denoting the specific matrix element. 
The E and H field component update are defined as follows.


 §.§ E-field update

E_x:

    cx    =   ξ_xx(i,j,k) D_x(i,j,k)
     
    cy    =   ξ_xy(i,j,k) D_y(i,j,k)
     
    cz    =   ξ_xz(i,j,k) D_z(i,j,k)
     
    E_x(i,j,k)   =    cx + cy + cz


E_y:

    cx    =   ξ_yx(i,j,k) D_x(i,j,k)
     
    cy    =   ξ_yy(i,j,k) D_y(i,j,k)
     
    cz    =   ξ_yz(i,j,k) D_z(i,j,k)
     
    E_y(i,j,k)   =    cx + cy + cz


E_z:

    cx    =   ξ_zx(i,j,k) D_x(i,j,k)
     
    cy    =   ξ_zy(i,j,k) D_y(i,j,k)
     
    cz    =   ξ_zz(i,j,k) D_z(i,j,k)
     
    E_z(i,j,k)   =    cx + cy + cz




 §.§ H-field update

H_x:

    cx    =   ζ_xx(i,j,k) B_x(i,j,k)
     
    cy    =   ζ_xy(i,j,k) B_y(i,j,k)
     
    cz    =   ζ_xz(i,j,k) B_z(i,j,k)
     
    H_x(i,j,k)   =    cx + cy + cz


H_y:

    cx    =   ζ_yx(i,j,k) B_x(i,j,k)
     
    cy    =   ζ_yy(i,j,k) B_y(i,j,k)
     
    cz    =   ζ_yz(i,j,k) B_z(i,j,k)
     
    H_y(i,j,k)   =    cx + cy + cz


H_z:

    cx    =   ζ_zx(i,j,k) B_x(i,j,k)
     
    cy    =   ζ_zy(i,j,k) B_y(i,j,k)
     
    cz    =   ζ_zz(i,j,k) B_z(i,j,k)
     
    H_z(i,j,k)   =    cx + cy + cz




§ AVERAGED CONSTITUTIVE RELATIONS UPDATE
 

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



 §.§ E-field update

E_x:

    cx    =   1/2(ξ_xx(i,j,k)+ξ_xx(i-1,j,k)) Dx_x(i,j,k)
     
    cy    =   1/4(ξ_xy(i,j,k) (D_y(i,j,k) + D_y(i,j+1,k)) + 
          ξ_xy(i-1,j,k) (D_y(i-1,j,k) + D_y(i-1,j+1,k)))
     
    cz    =   1/4(ξ_xz(i,j,k) (D_z(i,j,k) + D_z(i,j,k+1)) +
          ξ_xz(i-1,j,k) (D_z(i-1,j,k) + D_z(i-1,j,k+1)))
     
    E_x(i,j,k)   =    cx + cy + cz


E_y:

    cx    =   1/4(ξ_yx(i,j,k) (D_x(i,j,k) + D_x(i+1,j,k)) + 
          ξ_yx(i,j-1,k) (D_x(i,j-1,k) + D_x(i+1,j-1,k)))
     
    cy    =   1/2(ξ_yy(i,j,k) + ξ_yy(i,j-1,k)) D_y(i,j,k)
     
    cz    =   1/4(ξ_yz(i,j,k) (D_z(i,j,k) + D_z(i,j,k+1)) + 
          ξ_yz(i,j-1,k) (D_z(i,j-1,k) + D_z(i,j-1,k+1))) 
     
    E_y(i,j,k)   =    cx + cy + cz



E_z:

    cx    =   1/4(ξ_zx(i,j,k) (D_x(i,j,k) + D_x(i+1,j,k)) + 
          ξ_zx(i,j,k-1) (D_x(i,j,k-1) + D_x(i+1,k,k-1)))
     
    cy    =   1/4(ξ_zy(i,j,k) (D_y(i,j,k) + D_y(i,j+1,k)) +
          ξ_zy(i,j,k-1) (D_y(i,j,k-1) + D_y(i,j+1,k-1)))
     
    cz    =   1/2(ξ_zz(i,j,k) + ξ_zz(i,j,k-1)) D_z(i,j,k)
     
    E_z(i,j,k)   =    cx + cy + cz




 §.§ H-field update

H_x:

    cx    =   1/4(ζ_xx(i,j,k)+ζ_xx(i,j-1,k)+ζ_xx(i,j,k-1)+ζ_xx(i,j-1,k-1)) B_x(i,j,k)
     
    cy    =   1/8((ζ_xy(i,j,k) + ζ_xy(i,j,k-1)) (B_y(i,j,k) + B_y(i+1,j,k)) +
                (ζ_xy(i,j-1,k) + ζ_xy(i,j-1,k-1)) (B_y(i,j-1,k) + B_y(i+1,j-1,k)))
     
    cz    =   1/8((ζ_xz(i,j,k) + ζ_xz(i,j-1,k)  ) (B_z(i,j,k) + B_z(i+1,j,k)) +
                   (ζ_xz(i,j,k-1) + ζ_xz(i,j-1,k-1)) (B_z(i,j,k-1) + B_z(i+1,k,k-1)))
     
    H_x(i,j,k)   =    cx + cy + cz


H_y:

    cx    =   1/8((ζ_yx(i,j,k) + ζ_yx(i,j,k-1)) (B_x(i,j,k) + B_x(i,j+1,k)) +
                   (ζ_yx(i-1,j,k) + ζ_yx(i-1,j,k-1)) (B_x(i-1,j,k) + B_x(i-1,j+1,k)))
     
    cy    =   1/4(ζ_yy(i,j,k) + ζ_yy(i-1,j,k)+ζ_yy(i,j,k-1) + ζ_yy(i-1,j,k-1)) B_y(i,j,k)
     
    cz    =   1/8((ζ_yz(i,j,k) + ζ_yz(i-1,j,k)) (B_z(i,j,k) + B_z(i,j+1,k)) +
                   (ζ_yz(i,j,k-1) + ζ_yz(i-1,j,k-1)) (B_z(i,j,k-1) + B_z(i,j+1,k-1)))
     
    H_y(i,j,k)   =    cx + cy + cz


H_z:

    cx    =   1/8((ζ_zx(i,j,k)+ζ_zx(i,j-1,k)) (B_x(i,j,k) + B_x(i,j,k+1)) +
                   (ζ_zx(i-1,j,k) + ζ_zx(i-1,j-1,k)) (B_x(i-1,j,k) + B_x(i-1,j,k+1)))
     
    cy    =   1/8((ζ_zy(i,j,k)+ζ_zy(i-1,j,k)) (B_y(i,j,k) + B_y(i,j,k+1)) +
                   (ζ_zy(i,j-1,k)+ζ_zy(i-1,j-1,k)) (B_y(i,j-1,k) + B_y(i,j-1,k+1)))
     
    cz    =   1/4(ζ_zz(i,j,k) + ζ_zz(i-1,j,k)+ζ_zz(i,j-1,k) + ζ_zz(i-1,j-1,k)) B_z(i,j,k)
     
    H_z(i,j,k)   =    cx + cy + cz


