Fast and accurate electromagnetic solutions of finite periodic optical structures

Alberto Serna; Mario F. Manzano; Luis Landesa; Diego M. Solis; Jose M. Taboada

doi:10.1364/OE.25.018031

1. Introduction

The electromagnetic analysis of large periodic structures is an interesting topic in the scope of computational electromagnetics, from antenna arrays, RF metamaterial and metasurfaces (based on periodic patterns of resonant or non-resonant cells) to new applications in the optical regime for plasmonic periodic structures, among others. Several methods have been proposed to deal with this kind of problems.

The Floquet’s harmonic theory [1, 2], based on series expansions, is useful for analysing the electromagnetic behaviour of periodic infinite body arrays. Although Floquet theory might also be used to approximate the electromagnetic response of a finite periodic structure, important issues in such systems, such as edge effects or stationary waves, would not be properly taken into account.

Solutions based on the discretization of the surface integral equations can be applied to solve large problems through iterative methods and acceleration techniques. Among them, the Multilevel Fast Multipole Algorithm (MLFMA) [3–9] is one of the most useful approaches. To obtain the solution, the MLFMA accelerates the Matrix Vector Product (MVP) in the framework of an iterative process. Its efficiency achieves 𝒪(N log(N)) (where N is the problem number of unknowns). The MLFMA does not take special advantage of periodic problems since it is designed for arbitrary shaped geometries. Another extended approach for fast MVP calculation is the precorrected Fast Fourier Transform (p-FFT) method [10, 11]. At the expense of a significant loss of accuracy, the p-FFT transforms a general problem in a new uniform three-dimensional (3D) grid alternative scheme. The solution can be accelerated using the Fast Fourier Transform (FFT). The CG-FFT method [12] is a similar procedure for accelerating a 3D uniform grid discretized problem using the FFT.

Apart from the use in general non periodic problems, the FFT has been applied to solve, accurately, some kind of periodic structures. For example, in some specific antenna array systems, the solution of the Method of Moments (MoM) [13] deals with multilevel block-Toeplitz matrices that can be solved in a fast way using the FFT [14–16]. For dealing with generic block-Toeplitz matrices, another FFT based approach is developed in [17, 18] that obtains fast solutions for periodic antenna arrays. By taking advantage of the periodic nature of the composite, techniques in [14–18] accelerate the MVP without any loss of accuracy. Also, the technique proposed in this paper has these same advantages in addition to be valid for general periodic composites.

By taking pieces (slots) of the MoM impedance matrix, the solution of the MVP can be partitioned in several portions (also slots). Each one of these partial MVP’s can be efficiently solved in a fast way using the FFT. This technique, so called slot FFT, is applied for obtaining, both, accurate and fast electromagnetic solutions of very and extremely large generic periodic composites, e.g. metamaterials or plasmonic periodic structures in the optical regime.

The slot FFT will be explained in detail in the following section. Two different applications in the optical regime are presented in the Results section. Finally, the significance of the accuracy in fast solutions of large periodic problems is explained in the Discussion section.

2. The slot FFT method for periodic problems

Let us consider a periodic electromagnetic problem of M bodies (Fig. 1) each of them identically modelled with the MoM and using n unknowns. Let Z be the N × N impedance matrix of the MoM, being N = nM the total number of unknowns of the problem [13,19].

Fig. 1 Periodic composite of M bodies.

Download Full Size | PDF

For the sake of simplicity, only the one dimensional (1D) periodicity will be explained. In this case, the impedance matrix can be block-wise expressed as:

Z = {\begin{matrix} Z_{1, 1} & Z_{1, 2} & Z_{1, 3} & Z_{1, 4} & \dots & Z_{1, M} \\ Z_{2, 1} & Z_{2, 2} & Z_{2, 3} & Z_{2, 4} & \dots & Z_{2, M} \\ Z_{3, 1} & Z_{3, 2} & Z_{3, 3} & Z_{3, 4} & \dots & Z_{3, M} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ Z_{M, 1} & Z_{M, 2} & Z_{M, 3} & Z_{M, 4} & \dots & Z_{M, M} \end{matrix}}_{{N \times N}},

where the submatrix Z_i,j is composed of n × n elements, representing the electromagnetic interaction of the radiating basis functions of the body j over the testing functions of body i. In addition, since bodies are identically modelled and realizing that the interactions between elements depend only on their relative location, the submatrices fulfil the following properties:

Z_{i, i} = Z_{1, 1}

Z_{i, j} = Z_{1, j - 1 + 1} if j > i

Z_{i, j} = Z_{i - j + 1, 1} if j < i

Replacing Eqs. (2)–(4) into Eq. (1), Z can be expressed as a block Toeplitz Matrix as follows:

Z = {\begin{matrix} Z_{1, 1} & Z_{1, 2} & Z_{1, 3} & Z_{1, 4} & \dots & Z_{1, M} \\ Z_{2, 1} & Z_{1, 1} & Z_{1, 2} & Z_{1, 3} & \dots & Z_{1, M - 1} \\ Z_{3, 1} & Z_{2, 1} & Z_{1, 1} & Z_{1, 2} & \dots & Z_{1, M - 1} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ Z_{M, 1} & Z_{M - 1, 1} & Z_{M - 2, 1} & Z_{M - 3, 1} & \dots & Z_{1, 1} \end{matrix}}_{{N \times N}},

To simplify the notation let Z^k be the matrix given by:

Z^{k} \equiv Z_{i, j} where k = j - i, - M + 1 < k < M - 1

Replacing Eq. (6) into Eqs. (2)–(4) we get

Z^{0} \equiv Z_{i, i} = Z_{1, 1}

Z^{j - i} \equiv Z_{i, j} = Z_{1, j - i + 1} if j > i

Z^{j - i} = Z^{- (i - j)} \equiv Z_{i, j} = Z_{i - j + 1, 1} if j < i

Hence, matrix Z given by Eq. (5) can be expressed as:

Z = {\begin{matrix} Z^{0} & Z^{1} & Z^{2} & Z^{3} & \dots & Z^{M - 1} \\ Z^{- 1} & Z^{0} & Z^{1} & Z^{2} & \dots & Z^{M - 2} \\ Z^{- 2} & Z^{- 1} & Z^{0} & Z^{1} & \dots & Z^{M - 3} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ Z^{- (M - 1)} & Z^{- (M - 2)} & Z^{- (M - 3)} & Z^{- (M - 4)} & \dots & Z^{0} \end{matrix}}_{{N \times N}},

If the Z^k elements were scalar (1 × 1 matrices), the MVP of martrix Z by a vector x would be solved as a convolution. In such a case, this MVP could be solved in 𝒪(N log N) operations by using the FFT to accelerate the convolution. Nevertheless, this is only an hypothetical case and we need to deal with general block Toeplitz matrices in order to accelerate the MVP.

Transforming Z from block-Toeplitz to a standard Toeplitz matrix, the FFT algorithm is directly applied slotting the matrix and the solution. This transformation is the key point of the slot FFT method.

Let R_i be the block row matrix satisfying:

R_{i} = {\begin{matrix} Z^{1 - i} & Z^{2 - i} & Z^{3 - i} & \dots & Z^{M - i} \end{matrix}}_{{n \times N}} 0 < i \leq M

Hence, matrix Z given in Eq. (10) can be written as:

Z = {\begin{matrix} R_{1} \\ R_{2} \\ R_{3} \\ ⋮ \\ R_{M} \end{matrix}}_{{N \times N}}

Next, we define the special matrix R of dimensions n × (2N − n) which is a container of the R_i matrices:

R \equiv {\begin{matrix} Z^{- (M - 1)} & Z^{- (M - 2)} & \dots & Z^{- 1} & Z^{0} & Z^{1} & Z^{2} & - & Z^{M - 1} \end{matrix}}

It is easy to see that the R_i matrices, 0 < i ≤ M, can be obtained from R taking columns from N − in + 1 to 2N − in. This property can be expressed using MATLAB notation as R_i = R(:, (N − in + 1) : (2N − in)).

Let 𝒵_k be the n injected N × N matrices, with 0 < k ≤ n, defined as:

{𝒵_{k}}_{i, j} \equiv R_{k, i - j - k + N + 1}

In general, i − j − k + N + 1 > 0 in the above equation. For the very few elements in which the evaluation of the injected matrix 𝒵_k element is out of the range of R we can take, for example, zero or Not-a-Number values.

Matrix 𝒵_k is built by using only the k-row of R. Notice that the k + in row (with i = 0, . . . , M − 1) of the 𝒵_k matrix is equal to the k + in row of the real original Z matrix, that is:

{𝒵_{k}}_{k + in, j} = Z_{k + in, j} 0 \leq i < M, 1 \leq j \leq N

In addition, each row of 𝒵_k can be obtained by shifting the precedent row as:

{𝒵_{k}}_{i, j} = {𝒵_{k}}_{i - 1, j - 1} 1 < i, j \leq N

which implies that each one of the 𝒵_k matrices (k = 1 . . . n) is a Toeplitz matrix.

Let S_k be the slotting matrix given by

{S_{k}}_{i, j} = {\begin{array}{l} 1 & if i = j AND i = k + \dot{n} \\ 0 & elsewhere \end{array}}

where ṅ represents a multiple of n. S_k is a matrix with M periodically distributed ones in the diagonal. It satisfies the following properties:

The product of two slotting matrices is: $S_{k} S_{l} = {\begin{array}{l} S_{k} & k = l \\ 0_{N \times N} & k \neq l \end{array}$
The sum of all the slotting matrices gives the identity matrix: $I = \sum_{k = 1}^{n} S_{k}$ hence, any matrix, A, can be written as: $A = I A = \sum_{k = 1}^{n} S_{k} A$
The S_k operator produces the same result when it is applied to Z and when it is applied to the applicable injected matrix. Using Eq. (15): $S_{k} Z = S_{k} 𝒵_{k}$

Summing over k in Eq. (19) we obtain:

\sum_{k = 1}^{n} S_{k} Z = \sum_{k = 1}^{n} S_{k} 𝒵_{k}

Finally, replacing Eq. (18) in Eq. (20), the matrix Z is obtained from slotted versions of the injected matrices 𝒵_k as follows

Z = \sum_{k = 1}^{n} S_{k} 𝒵_{k}

Then the MVP of the matrix Z by a vector x can be calculated using the Eq. (21) as:

Z x = \sum_{k = 1}^{n} [S_{k} (𝒵_{k} x)]

The computational cost of Eq. (22) is dominated by the product 𝒵_kx (the product of S_k by a vector is a trivial replacement by zeros). Given that 𝒵_k is a Toeplitz matrix, 𝒵_k x is a convolution and it can be solved in a expedited way using the FFT in 𝒪(N log N) operations. Since there are n MVP operations in Eq. (22), the total cost of evaluating the product Zx is 𝒪(nN log N).

Notice that no calculation of injected matrices, 𝒵_k, is required in the slot FFT method. Also, the calculation of the MoM matrix, Z, is not required either. The method computes only one matrix, R, expressed in Eq. (13), i.e. it needs to store n × (2N − n) elements in memory.

In short, if n, the number of unknowns of each cell element, were negligible compared with M, the method would be virtually of order 𝒪(N log N) and 𝒪(N) in terms of CPU time and memory requirement, respectively.

3. Results

The following results were performed in a High Performance Computing (HPC) system composed by 4 Intel Xeon E7-4820 at 2.00GHz (eight cores each one, 32 cores in all), with 512GB of RAM memory.

In Fig. 2(a) the memory cost of the whole MoM and the slot FFT using periodic structures in 1D, 2D and 3D uniform grids is showed. These results are obtained for the scattering of a λ = 540nm wavelength electromagnetic wave impinging on 1D, 2D and 3D periodic compositions of gold nanospheres of radius 100nm and 280nm separation between centres; each sphere is modelled with n = 240 unknowns and using a number of spheres, M, from just a few units to 40000. A significant memory saving can be observed in the slot FFT, given that only R in Eq. (13) needs storage. In Fig. 2(b) the mean computing time of the MVP is shown. The differences between the 𝒪(nN log N) cost of slot FFT and the 𝒪(N²) direct cost of carrying out the MVP are observed. For more than 100 · 10³ unknowns, costs for the direct MoM MVP solution are extrapolated for both memory usage and MVP time. Notice that in order to avoid bias in the comparison, these results were performed without parallelization.

Fig. 2 Comparison of the memory and computing time needed with MoM and slot FFT methods.

Download Full Size | PDF

The first example is the electromagnetic scattering analysis at a wavelength of λ = 550nm from a periodic dense metastructure composed of 110 × 110 × 110 nanocubes in a uniform 3D grid. The edge dimension of the nanocubes is 55nm and the separation between centres is 110nm. The composite parameters are ∊_r = 3 and μ_r = 1 and the nanocubes are suspended in vacuum medium. The problem is composed of nearly 50 million unknowns. In Fig. 3 the Radar Cross Section (RCS) of the metastructure is showed. The RCS is calculated for normal plane wave incidence (θ = 90°) using the slot FFT to speed up the MVP. This result is compared to the RCS of two different homogeneous dielectric cubes of edge dimension 12.1μm and ∊_r = 1.1579 and ∊_r = 1.1638 (effective permittivity obtained with the Maxwell-Garnett and Bruggeman homogenization method [20], respectively). The computation of the RCS for the two entire homogeneous cubes were performed using the MLFMA algorithm [5]. The results of the homogeneous effective cube using MLFMA needed 1, 7 million unknowns. This result proves that both homogenization formulas (Maxwell-Garnett and Bruggeman) provide accurate results for general purposes. Nevertheless, in this kind of large problems, we can make use of the full-wave solution based on the slot FFT for a more accurate result. The three results show a good agreement. Notice that 30 degrees apart from the principal directions a slight difference between the electromagnetic scattering of the two homogeneous cubes and the real electromagnetic response of the large periodic composite object is revealed. The properties of the real and homogeneous material differs at the microscopic scale, which is important to describe the exact behaviour of the material.

Fig. 3 Scattering diagram of a 110³ cube-shaped array and its equivalent homogenized by the Maxwell-Garnett and Bruggeman equivalent permittivities.

Download Full Size | PDF

As second example, a plane-shaped array composed by gold nanospheres is analysed. In this example we study the surface-enhancement Raman spectroscopy (SERS) [21–23] to show the effects that appear in large SERS substrates. Some important effects in SERS may not be discovered with other widespread analysis techniques, e.g. approximating the structure by an infinite array or analysing it by a small portion of it, this giving a wrong prediction of the behavior of the substrate in a real scenario. This plays an important role in these kind of applications, as the SERS intensity distribution may present stationary wave-patterns effects.

The studied array was composed of 200 × 200 gold nanospheres located in an uniform 2D grid and the simulation was performed at a wavelength of 600nm, in a medium filled with water (∊_r = 1.7689) and illuminated by two incident plane waves perpendicular to the array, each of them linear polarized along the principal axes. The diameter of the nanospheres was 55nm and the separation between centres was 56nm. The discretization of the problem resulted in nearly 10 million unknowns.

For the sake of simplicity, SERS intensity was calculated using no Raman shift. This approach is valid for small Raman shift and yields close-fitting results in a faster way through the following expression [23]:

SERS enhancement \approx {| E (ω_{in}) |}^{4} {| E_{i} (ω_{in}) |}^{- 4},

where E is the total electric field in the presence of the SERS substrate and E_i the field of the incident planewave. Figure 4(a) shows the SERS enhancement obtained after smoothing with a 2D Gaussian filter—to simulate the effect of a lens—. The representation is saturated so as to give a better illustration of the effects that appear along the xy plane. It is clearly discernible how the edge effect and the stationary wave pattern appear. This makes highly recommended the analysis of real case substrates.

Fig. 4 Representation of the SERS intensity of a small part of the composite object.

Download Full Size | PDF

Figure 4(b) shows the SERS of a small part of the structure without filtering. The hotspots between nanospeheres can be appreciated, with peaks of 1.0824 · 10⁵ arbitrary units. This figure also shows a saturated representation of the SERS intensity to give a better representation of the hotspots.

4. Discussion

The presented approach, based on slotting partial FFT solutions (Eq. 22), is a very efficient way to compute the solution of very or extremely large electromagnetic finite periodic problems. The amount of required memory of 𝒪(nN) and the 𝒪(nN log N) computing time allow us to analyse complicated periodic problems in the optical regime without accuracy loss. Although, for the sake of clarity and conciseness, only the 1D approach is formulated in this paper, their 2D and 3D counterparts can be formulated in an analogous fashion. In fact, in the results Section the use of 2D and 3D approaches were necessary. The method is straightforward to apply for any MoM implementation. Also, parallelization of the slot FFT does not require too much effort. As a matter of fact, the OpenMP parallelization was used to obtain the results in Figs. 3 and 4. Furthermore, future combination with MPI can be used to obtain electromagnetic solutions to extremely large problems in supercomputer systems.

These results support the advantages of slot FFT, suggesting this method as a feasible approach to tackle the upcoming challenges in the fields of nanomaterials and other areas where periodic structures are present.

The importance of having full-wave techniques for analysing very large periodic structures is showed in the results section. Using the slot FFT, the edge effects are well modelled and stationary waves are well predicted. These effects might not be discovered approximating the structure by an infinite array or a by a small portion of it. In addition, some scattering deviations, between equivalent homogeneous models and the real metamaterial, are revealed with the slot FFT technique.

Funding

European Regional Development Fund (ERDF) and the Spanish Ministerio de Economia y Competitividad (MAT2014-58201-C2-1-R); ERDF and the Extremadura Regional Government (GR15052); Spanish Ministerio de Educacion, Cultura y Deporte (FPU 00022/15).

References and links

1. M. E. Veysoglu, R. T. Shin, and J. A. Kong, “A Finite-Difference Time-Domain analysis of wave scattering from periodic surfaces: Oblique incidence case,” J. Electromagn. Waves Appl. 7, 1595–1607 (1993). [CrossRef]

2. P. Yeh, A. Yariv, and C.-S. Hong, “Electromagnetic propagation in periodic stratified media. I. General theory,” J. Opt. Soc. Am. 67, 423–438 (1977). [CrossRef]

3. O. Ergul and L. Gurel, The Multilevel Fast Multipole Algorithm (MLFMA) for Solving Large-Scale Computational Electromagnetics Problems (Wiley-IEEE Press, 2014). [CrossRef]

4. O. Ergul and L. Gurel, “Accurate solutions of extremely large integral-equation problems in computational electromagnetics,” Proc. IEEE 101, 342–349 (2013). [CrossRef]

5. J. M. Taboada, M. G. Araujo, F. Obelleiro, J. L. Rodriguez, and L. Landesa, “MLFMA-FFT parallel algorithm for the solution of extremely large problems in electromagnetics,” Proc. IEEE 101, 350–363 (2013). [CrossRef]

6. O. Ergul and L. Gurel, “Fast and accurate analysis of large-scale composite structures with the parallel Multilevel Fast Multipole Algorithm,” J. Opt. Soc. Am. A 30, 509–517 (2013). [CrossRef]

7. J. M. Taboada, M. G. Araujo, J. M. Bertolo, L. Landesa, F. Obelleiro, and J. L. Rodriguez, “MLFMA-FFT parallel algorithm for the solution of large-scale problems in electromagnetics,” Prog. Electromagn. Res. 105, 15–30 (2010). [CrossRef]

8. O. Ergul, “Analysis of composite nanoparticles with surface integral equations and the Multilevel Fast Multipole Algorithm,” J. Opt. 14, 062701 (2012). [CrossRef]

9. M. G. Araujo, J. M. Taboada, J. Rivero, D. M. Solis, and F. Obelleiro, “Solution of large-scale plasmonic problems with the multilevel fast multipole algorithm,” Opt. Lett. 37, 416–418 (2012). [CrossRef] [PubMed]

10. X.-C. Nie, L.-W. Li, and N. Yuan, “Precorrected-FFT algorithm for solving combined field integral equations in electromagnetic scattering,” J. Electromagnet. Wave 16, 1171–1187 (2002). [CrossRef]

11. V. Okhmatovski, M. Yuan, I. Jeffrey, and R. Phelps, “A three-dimensional precorrected FFT algorithm for fast Method of Moments solutions of the mixed-potential integral equation in layered media,” IEEE Trans. Microw. Theory Tech. 57, 3505–3517 (2009). [CrossRef]

12. M. F. Catedra, R. P. Torres, J. Basterrechea, and E. Gago, The CG-FFT Method: Application of Signal Processing Techniques to Electromagnetics (Artech House, 1995).

13. R. F. Harrington, Field Computation by Moment Methods (Wiley-IEEE Press, 1993). [CrossRef]

14. B. E. Barrowes, F. L. Teixeira, and J. A. Kong, “Fast algorithm for matrix-vector multiply of asymmetric multilevel block-Toeplitz matrices in 3-D scattering,” Microw. Opt. Tech. Lett. 31, 28–32 (2001). [CrossRef]

15. Y. Zhang, J. Y. Li, C. H. Liang, and Y. J. Xie, “A fast analysis of large finite-slot phased arrays fed by rectangular waveguides,” J. Electromagnet. Wave 18, 715–727 (2004). [CrossRef]

16. H. Zhai, Q. Chen, Q. Yuan, K. Sawaya, and C. Liang, “Analysis of large-scale periodic array antennas by CG-FFT combined with equivalent sub-array preconditioner,” IEICE Trans. Electron. E89B, 922–928 (2006). [CrossRef]

17. E. H. Bleszynski, M. K. Bleszynski, and T. Jaroszewicz, “Block-Toeplitz fast integral equation solver for large finite periodic and partially periodic antenna arrays,” in “2003 IEEE Topical Conference on Wireless Communication Technology,” (2003), pp. 428–429.

18. E. H. Bleszynski, M. K. Bleszynski, and T. Jaroszewicz, “Block-Toeplitz fast integral equation solver for large finite periodic and partially periodic array systems,” IEICE Trans. Electron. E87C, 1586–1594 (2004).

19. J. M. Taboada, J. Rivero, F. Obelleiro, M. G. Araujo, and L. Landesa, “Method-of-moments formulation for the analysis of plasmonic nano-optical antennas,” J. Opt. Soc. Am. A 28, 1341–1348 (2011). [CrossRef]

20. A. J. Pickles and M. B. Steer, “Effective permittivity of 3-D periodic composites with regular and irregular inclusions,” IEEE Access 1, 523–536 (2013). [CrossRef]

21. S. Schlucker, “Surface-enhanced raman spectroscopy: Concepts and chemical applications,” Angew. Chem. Int. Ed. Engl. 53, 4756–4795 (2014). [CrossRef] [PubMed]

22. R. A. Alvarez-Puebla and L. M. Liz-Marzan, “SERS-based diagnosis and biodetection,” Small 6, 604–610 (2010). [CrossRef] [PubMed]

23. D. M. Solis, J. M. Taboada, L. Landesa, J. L. Rodriguez, and F. Obelleiro, “Squeezing Maxwell’s equations into the nanoscale,” Prog. Electromagn. Res. 154, 35–50 (2015). [CrossRef]

Fast and accurate electromagnetic solutions of finite periodic optical structures

Abstract

1. Introduction

2. The slot FFT method for periodic problems

3. Results

4. Discussion

Funding

References and links

Cited By

Figures (4)

Equations (25)

Optics Express