A space-domain integral equation method accelerated with adaptive cross approximation (ACA) is presented for the fast and accurate analysis of electromagnetic (EM) scattering from multilayered metallic photonic crystal (MPC). The method directly solves for the electric field in order to easily enable the periodic boundary condition (PBC) in the spatial domain. The ACA is a purely algebraic method allowing the compression of fully populated matrices; hence, its formulation and implementation are independent of integral equation kernel (Green’s function). Therefore, the ACA is very well suited for accelerating integral equation analysis of periodic structure with the integral kernel of the periodic Green’s function (PGF). The computation of the spatial-domain periodic Green’s function (PGF) is accelerated by the modified Ewald transformation, such that the multilayered periodic structure can be analyzed efficiently and accurately. An effective interpolation method is also proposed to fast compute the periodic Green’s function, which can greatly reduce the time of matrix filling. Numerical examples show that the proposed method can greatly save the frequency sweep time for multilayered periodic structure.

1. Introduction

Photonic crystals [1] are periodic structures of great interest for their applications both in the microwave region and in the optical range. The main feature is the presence of frequency bands wherein the waves are highly attenuated and cannot propagate. It results from the removal of degeneracies of the free-photon states at the Bragg planes provoked by the periodicity, which produces forbidden frequency gaps so-called photonic band gaps (PBGs). This property is exploited in the electromagnetic and optical applications.

Many numerical methods for computing such band structures have been adapted from solid-state physics. They include, but are not restricted to, plane-wave expansions [2], the Korringa-Kohn-Rostoker method [3], and shell methodologies [4]. Over the past two decades, finite difference time domain (FDTD) [5] has been widely used in computational optics and photonics. Recent developments include specialized versions of the finite element method (FEM) [6] and the flexible local approximation method [7]. FDTD and FEM can cope with arbitrary complicated shapes and materials. However, the radiation into unbounded regions requires either absorbing boundary conditions or perfectly matched layers, and special measures need to be taken into account when employing these conditions to the scattered field formulation. For open-region problems, integral equation/method of moments (IE/MoM) has a great advantage.

The ACA is a purely algebraic method, which is independent of the kernel. The ACA was originally introduced in 2000 by Bebendorf [8] and since then has been successfully applied for integral formulations involved electromagnetic scattering [9, 10]. The integral kernel is a free space Green’s function. In addition, it is noted that most photonic crystal applications deal with two-dimensional (2D) structures that are invariant along a longitudinal axis and periodic in the transverse plane [11]. The manufacture of a 2D photonic crystal structure is easier than that of three-dimensional (3D) one [12].

In this paper, a space-domain integral equation accelerated with ACA is presented for the fast and accurate analysis of electromagnetic scattering from 2D and 3D multilayered periodic structures. The ACA is first introduced for surface integral equation with the kernel of periodic Green’s function (PGF). The ACA can be applied to compress impedance matrix blocks and accelerate matrix-vector multiplication (MVP) for reducing computational complexity. A modified Ewald transformation is proposed to efficiently compute the PGFs. Compared with traditional Ewald method, it has less computational cost and can eliminate the imbalance of two Ewald sequences for multilayered periodic structure. An efficient interpolation method is also proposed to fast compute the PGFs and then accelerate impedance matrix filling. Our proposed method can effectively reduce the frequency sweep time for multilayered periodic structure involving a lot of unknowns.

2. Formulation

Figure 1 shows the multilayered double periodic structure with identical metallic objects of arbitrary shape periodically repeated in the -plane. and are the primitive lattice vectors. is the cross-sectional area of the unit cell. is the incident electric field. is the polarization angle. The wave number is and . The wave vector is and phase shift . The array is assumed to be periodic in the -plane and the cell of the array is obtained by shifting the cell through the relation . is defined as the translation vector of the lattice.

2.1. Integral Equation

In this paper, we consider the metallic photonic crystal. Let be the surface of a metallic object. By enforcing the boundary conditions that the total tangential electric field is zero on the perfect electric conductor (PEC) surface , the electric field integral equation (EFIE) is given as [13]where is a doubly periodic Green’s function. The Rao-Wilton-Glisson (RWG) basis functions on triangular elements are employed to discretize the electric current [14]. In addition, the special basis and test functions derived from RWG basis function are introduced to ensure current continuity across the periodic boundaries [15]. The application of the Floquet-Bloch theorem reduces the computational domain of infinite periodic structures to a single unit cell but leads to the numerical evaluation of very slowly converging series.

2.2. Modified Ewald Transformation

We propose the following strategy to compute the spatial periodic Green’s function in an arbitrary point .

If , PGF is computed by its spectral representation with exponential convergence, and thus no acceleration technique is required [16]. Considerwhereis the reciprocal lattice vector. Considerwhere , .

If , PGF is computed by the representation of Ewald transformation [17], which are two exponential converging series. Considerwherewhere is the complementary error function of complex argument and .

By analyzing the asymptotic behavior of the series terms, the approximation to the optimal value for general skewed lattices is given bywhere is the maximum exponent permitted. In almost all practical cases it is sufficient to include only nine terms in (6) and (7) (i.e., the summation limits are from −1 to +1), in which the error is usually less than 0.1%.

For multilayered periodic structure, when , the PGFs are computed by their spectral representation rather than direct Ewald transformation. In doing so there are two advantages.(1)The spectral representation just requires evaluation of the exponential function, which is much simpler than the computation of complex-argument complementary error function (). Compared with , the computational complexity of the exponential function can be ignored.(2)In addition, it can eliminate the imbalance of two Ewald sequences when is too large. For the multilayered structure, the complementary error functions in the first terms of both Ewald series will take large imaginary arguments when is too large. These two terms will have very large values that are approximately equal in amplitude, but of opposite signs, so that summing them up will lead to a severe loss of accuracy due to the finite machine precision.

2.3. Singularity Extraction

When the distance between the observation point and the source point or one of its periodic extensions is electrically small, the contribution to Green’s function is essentially quasistatic. Usually the center distance between source and field triangles is less than ; then both are considered to be “near interaction.” Singular terms in (5) and (7) often arise, for example, from one of nine spatial terms with or +1, depending on which term represents the source nearest to the observer . By extracting its singularity [18], the PGF can be expressed as . The singularity subtraction technique [19] used in this paper can compute the singularity of () type.

Supposing the source point nearest to the observer point appears at term rather than term as shown in Figure 2, the primary source point can be transferred to lattice. Then, the singular term always appears at term . The singularity of can be extracted aswhere is the nonsingular partAs , the result of applying L’ Hôpital’s rule on (10) is as follows:

2.4. Interpolation of PGF

PGFs are very smooth and amenable to interpolation. Fast interpolation of PGFs can greatly reduce the time of matrix filling. The reason is that interpolating the PGF values from the table is much faster than directly computing the Ewald sums (since this has to be done for every pair of source and field triangles in the numerical integration using cubature formulas). Trilinear interpolation [20] suited for the calculation of complex-valued function with very high efficiency is applied here. As shown in Figure 3, the cuboid grouping the unit cell is equally divided into many subgrids with side length . Yellow points are the precomputed PGFs at the corner of the subgrid. Then the three-dimensional PGF values can be obtained by trilinear interpolation: where , .

The relative error obtained by interpolating the PGFs is defined aswhere and are the interpolated and the exact data, respectively. The exact data are obtained by the above accelerating method with an accuracy . The relative errors are less than and have, therefore, been neglected.

Take the multilayered structure with orthogonal lattices as an example. The relative error is shown in Figure 4 for the interpolation of PGF involving a series of interpolation points through the trilinear interpolation method, with various values of . In Figure 4, positions () are the sampling points for interpolation. A very high level of accuracy is easily obtained with a relatively small number of interpolation points. The side length of the subgrid is usually fixed to , and the higher interpolation accuracy can be obtained.

3. ACA Technique

Let the rectangular matrix represent the coupling between two well-separated groups in the MoM computation. The ACA algorithm aims to accurately approximate by a low-rank matrix . In particular, the ACA algorithm constructs the approximate matrix through a product form. Namely,where is the rank of or the effective rank of the matrix ; and are two dense rectangular matrices; is the column th of ; and is the row th of , respectively. The goal of ACA is to achievefor a given tolerance . Note that in this paper refers to the matrix Frobenus norm [9].

The ACA technique is numerically based on the principle of the approximation above but does not need to know all terms of . Indeed, only rows and columns are calculated during the algorithm allowing a significant reduction in the number of integral computations to perform. Therefore, not only is the memory requirement reduced but the assembly time is also saved. Yet, if we agree to deal with the MVP instead of the exact value of , the number of operations needed is also remarkably reduced by performing . Thus, if an iterative solver is used, one can accelerate matrix vector multiply and save once again the computation time [10].

The compression rate η of matrix block , which denotes the memory saved by using ACA, can be expressed as Obviously, the approximation by ACA approach above is useful only if , where presents similar values. For the electromagnetic application, the effective rank . Therefore, instead of storing entire entries, the algorithm only requires storing entries.

4. Numerical Results

The EFIE and ACA approaches to PEC are implemented using Fortran-95 language. Numerical experiments run on a PC with Intel i5-2400 3.1 GHz CPU and 16 GB RAM. The resulting impedance matrices are iteratively solved by using the transpose-free quasi-minimal residual (TFQMR) solver with diagonal preconditioning [21]; the relative error tolerance is set to be 10−3.

4.1. Infinite Metallic Rod

As a 2D example, we study the scattering properties of a photonic band gap (PBG) structure consisting of infinite metallic rod [17]. Each unit cell consists of two infinitely long metallic rods. The inset of Figure 5 is two meshed cylinders in a unit cell. The axis of each cylinder is along -direction: 6 mm and 6 mm. Each cylinder has the radius 0.6 mm and length 6 mm. The space between two cylinders is 6 mm. Each cylinder was discretized into 494 triangle facets with 741 unknowns to ensure accurate results obtained in the entire frequency band. Figure 5 shows the magnitude of the transmission and reflection coefficients for the Floquet TM00 mode with the electric field oriented in the -direction. Our results agree well with the reference results obtained using 2D IE approach [15]. The results converged for 9 terms in both Ewald sums and there are 1482 unknowns for the geometry. The time for matrix filling and solving is 0.58 s and 0.21 s per frequency point, respectively. In ACA, the box size is mm3; that is, each box includes one cylinder. If the trilinear interpolation method is not used, the time of matrix filling would increase dramatically to 128.6 s. Therefore, trilinear interpolation of PGFs can reduce the time of matrix filling by more than 99.5%.

Next, the number of layers increases from 2 to 5, 10 and 21, respectively. Figure 6 shows the transmission coefficient of multilayered structure. The transmission coefficient curves exhibit the typical resonances of photonic band gap materials. It can be seen that the null of transmission becomes deeper as the layer number increases. The total solution times for per frequency point are 5.1 s, 12.4 s, and 32.5 s, respectively.

The 21-layer structure is meshed into 10374 triangle facets with 15561 unknowns. The compression rate of matrix block is shown in Table 1. The times for matrix filling and matrix solving are 23.8 s and 8.7 s per frequency point, respectively. The corresponding times of direct solution by MoM without ACA acceleration are 187.6 s and 75.2 s per frequency point, respectively. Obviously, ACA can greatly reduce the total solution time, especially for frequency sweeping.

4.2. Woodpile

A 3D example showing the inset in Figure 7 is a woodpile MPC with  nm,  nm and in-layer lattice constant μm. The incident wave is a TE polarized (i.e., parallel to the top metallic layer) plane wave through the stacking direction. The number of layers are 4, 8, and 16, respectively. Two intersecting metallic pipes are considered to be fully contacting with each other in numerical modeling. The ACA box size is . The reflection coefficients are shown in Figure 7. Almost the same result is obtained for TM polarization. Our numerical results agree well with that of Ansys HFSS. Previous results reported by other studies [12, 22] for a similar kind of the structural geometry also validate our calculated results. For the 16-layer MPC with 13149 unknowns (8766 triangle facets), the times for matrix filling and linear system solving are 40.3 s and 15.2 s per frequency point, respectively. The corresponding times of direct solution by MoM without ACA acceleration are 132.3 s and 92.5 s, respectively.

5. Conclusion

An integral equation accelerated with adaptive cross approximation (ACA) is proposed for the analysis of multilayered metallic photonic crystal. The ACA is a purely algebraic and kernel independent algorithm. It is verified to be suited for accelerating integral equation analysis of periodic structure. The ACA can be applied to compress impedance matrix blocks and accelerate MVP for reducing computational complexity. A modified Ewald transformation and an efficient interpolation method are also proposed to fast compute the periodic Green’s function; hence, that can greatly reduce the time of matrix filling. Numerical results confirm the validity and accuracy of the proposed method. In our next stage of research, ACA will be applied to accelerate hybrid field integral equation for arbitrarily complex multilayered periodic structure composed of metal and dielectric materials with more unknowns.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


The work was supported in part by the National Science Foundation of China (NSFC) under Grants no. 61331002, no. 61102011, and no. 61201082 and in part by Excellent Innovation Team of CUC under Grant no. yxtd201303.