Three-dimensional adaptive coordinate transformations for the Fourier modal method

The concepts of adaptive coordinates and adaptive spatial resolution have proved to be a valuable tool to improve the convergence characteristics of the Fourier Modal Method (FMM), especially for metallo-dielectric systems. Yet, only two-dimensional adaptive coordinates were used so far. This paper presents the first systematic construction of three-dimensional adaptive coordinate and adaptive spatial resolution transformations in the context of the FMM. For that, the construction of a three-dimensional mesh for a periodic system consisting of two layers of mutually rotated, metallic crosses is discussed. The main impact of this method is that it can be used with any classic FMM code that is able to solve the large FMM eigenproblem. Since the transformation starts and ends in a Cartesian mesh, only the transformed material tensors need to be computed and entered into an existing FMM code. © 2014 Optical Society of America OCIS codes:(050.1970) Diffractive optics; (160.5298) Photonic crystals; (160.3918) Metamaterials; (050.1755) Computational electromagnetic methods. References and links 1. K. Busch, G. von Freymann, S. Linden, S. F. Mingaleev, L. Tkeshelashvili, and M. Wegener, “Periodic nanostructures for photonics,” Phys. Rep. 444, 101–202 (2007). 2. L. Li, “Formulation and comparison of two recursive matrix algorithms for modeling layered diffraction gratings,” J. Opt. Soc. Am. A13, 1024–1034 (1996). 3. P. Lalanne and G. M. Morris, “Highly improved convergence of the coupled-wave method for TM polarization,” J. Opt. Soc. Am. A13, 779–784 (1996). 4. G. Granet and B. Guizal, “Efficient implementation of the coupled-wave method for metallic lamellar gratings in TM polarization,” J. Opt. Soc. Am. A13, 1019–1023 (1996). 5. L. Li, “Use of Fourier series in the analysis of discontinuous periodic structures,” J. Opt. Soc. Am. A 13, 1870– 1876 (1996). 6. G. Granet, “Reformulation of the lamellar grating problem through the concept of adaptive spatial resolution,” J. Opt. Soc. Am. A16, 2510–2516 (1999). 7. G. Granet and J.-P. Plumey, “Parametric formulation of the Fourier modal method for crossed surface-relief gratings,” J. Opt. A4, S145–S149 (2002). 8. T. Vallius and M. Honkanen, “Reformulation of the Fourier modal method with adaptive spatial resolution: application to multilevel profiles,” Opt. Express 10, 24–34 (2002). 9. T. Weiss, G. Granet, N. A. Gippius, S. G. Tikhodeev, and H. Giessen, “Matched coordinates and adaptive spatial resolution in the Fourier modal method,” Opt. Express 17, 8051–8061 (2009). 10. S. Essig and K. Busch, “Generation of adaptive coordinates and their use in the Fourier Modal Method,” Opt. Express18, 23258–23274 (2010). 11. J. Küchenmeister, T. Zebrowski, and K. Busch, “A construction guide to analytically generated meshes for the Fourier Modal Method,” Opt. Express 20, 17319–17347 (2012). #199335 $15.00 USD Received 11 Oct 2013; revised 9 Dec 2013; accepted 13 Dec 2013; published 14 Jan 2014 (C) 2014 OSA 27 January 2014 | Vol. 22, No. 2 | DOI:10.1364/OE.22.001342 | OPTICS EXPRESS 1342 12. M. Decker, M. Ruther, C. E. Kriegler, J. Zhou, C. M. Soukoulis, S. Linden, and M. Wegener, “Strong optical activity from twisted-cross photonic metamaterials,” Opt. Lett. 34, 2501–2503 (2009). 13. V. Liu and S. Fan, “S 4: A free electromagnetic solver for layered periodic structures,” Comput. Phys. Commun. 183, 2233–2244 (2012). 14. H. Kim, J. Park, and B. Lee, Fourier Modal Method and its Applications in Computational Nanophotonics (CRC Press, 2012). 15. L. Li, “Fourier modal method for crossed anisotropic gratings with arbitrary permittivity and permeability tensors,” J. Opt. A5, 345–355 (2003). 16. A. Vial, A.-S. Grimault, D. Macı́as, D. Barchiesi, and M. Lamy de la Chapelle, “Improved analytical fit of gold dispersion: Application to the modeling of extinction spectra with a finite-difference time-domain method,” Phys. Rev. B71, 085416 (2005).


Introduction
Periodic nanostructures gathered a tremendous amount of interest in the past decade [1].Evolving experimental techniques allowed for more and more complex structures.Alongside this experimental development, numerical tools to solve Maxwell's equations became much more elaborate, too.
One of these rapidly developing numerical techniques is the Fourier Modal Method (FMM).It is capable to predict the transmission properties of periodic photonic systems, both dielectric and metallic.The systems normally considered are periodic with respect to the xy-plane and finite in z-direction.The system is sliced into layers with constant permittivity in z-direction and in each of these layers an eigenvalue problem is solved which stems from Maxwell's curl equations.This allows expanding the fields into eigenmodes.The layers are then connected using a scattering matrix algorithm which ensures the fulfillment of the continuity conditions [2].
After significant advancements for lamellar gratings and on the topic of the correct Fourier factorization rules [3][4][5], the FMM still faced the problem of properly calculating the response of metallic systems.In-plane stair-casing for not-grid-aligned structures and most of all the Gibbs phenomenon remained a challenge.
A popular approach to tackle these problems is the application of coordinate transformations.Two kinds of transformations emerged.First, adaptive coordinates (AC) are used to transform the permittivity distributions in such a way that the surface of the considered structures becomes grid-aligned.Second, adaptive spatial resolution (ASR) increases the coordinate line density along the interfaces.Combining both drastically enhance the performance of the method [6][7][8].In recent years, different concepts emerged for the construction of the corresponding meshes [9][10][11].
So far, adaptive coordinates have been used in the xy-plane.However, complex structures occurring in different layers pose a problem since different adaptive meshes would be necessary.How to connect these different meshes optimally remains a challenging task since each mesh represents a different basis.Also, the incident plane waves need to be transformed which induces additional errors.These problems can be tackled by designing a three-dimensional adaptive coordinate transformation.This method trades an increased amount of slices in the method for an accurate representation of the structure's surface in all three dimensions.In this paper, such a three-dimensional transformation is designed for a system that has gathered an extensive amount of interest in recent years, two periodic layers of mutually rotated, metallic crosses [12].These structures are known for their strong optical activity.The transformation leads to fully anisotropic permittivity and permeability tensors.The key point is that the metallic crosses are twisted by the transformation so that they are both grid-aligned in the transformed space (thus allowing an ideal representation by the FMM) and at the same time the transformation between the crosses is performed continuously in z-direction.
Moreover, a way to implement a three-dimensional adaptive spatial resolution is discussed.
The overall result is a method that allows the implementation of arbitrary three-dimensional adaptive coordinates and adaptive spatial resolution in classical, open source FMM codes [13,14].
In section 2, the basics of Maxwell's equations in generalized coordinates are discussed.Section 3 covers the design of the three-dimensional adaptive coordinate transformation.Finally, the concept for three-dimensional adaptive spatial resolution is presented in section 4.

Covariant formulation of the Fourier Modal Method with generalized coordinates
In this section, we discuss how generalized coordinates are incorporated in the FMM.The system we investigate in this paper is displayed in Fig. 1.It consists of two periodic layers of mutually rotated, metallic crosses.Fig. 1.Schematic view of the system of interest, namely two layers of mutually rotated, metallic crosses.The system is finite in x 3 -direction and periodic with respect to the x 1 x 2 plane.k in denotes the wave vector of an incident plane wave.Ox 1 x 2 x 3 describes a Cartesian coordinate system.
Since generalized coordinates in the context of the FMM have been discussed before, we will only briefly discuss them here.The presentation naturally follows previous publications on the topic in content and notation, see [9][10][11]15].We distinguish between a curvilinear coordinate system Ox 1 x 2 x 3 and a Cartesian coordinate system Ox 1 x 2 x 3 .The three-dimensional adaptive coordinate transformations that are investigated in this paper have the form Eventually, we want to solve Maxwell's curl equations which read in covariant form, see [15].Here, ξ denotes the Levi-Civita symbol and E σ and H σ are covariant components of the electric and magnetic field.Throughout the manuscript, Greek indices run from 1 to 3. Furthermore, we use the Einstein sum convention, meaning that repeated indices are implicitly summed over.The vacuum wave number is denoted k 0 = ω/c with the frequency ω and the speed of light c.The metric tensor reads and g (as used in Eqs. ( 4) and ( 5)) denotes the reciprocal of its determinant.As illustrated in detail in [11], the coordinate transformation leads to a transformed permittivity of the form where ε τ χ is the permittivity tensor in the Cartesian system.The permeability transforms identically.It is noteworthy that the matrix that is Fourier transformed when using FMM with AC and/or ASR is not the permittivity itself given in Eq. ( 7) but rather √ gε ρσ .Therefore, we refer to √ gε ρσ from now on as the effective permittivity.By entering Eqs. ( 1)-(3) in Eq. ( 7) one can observe that both the effective permittivity and the effective permeability become fully anisotropic.Therefore, the full anisotropic FMM eigenvalue problem has to be solved.

Three-dimensional adaptive coordinates
The overall aim of this section is to obtain three-dimensional adaptive coordinate transformation of the form in Eqs. ( 1) and (2) for our system depicted in Fig. 1.Since three-dimensional meshing is a complex procedure, the discussion is structured the following way: First, an example of a two-dimensional mesh is discussed.Second, this planar mesh is utilized to create a three-dimensional mapping.Third, the impact on the transformed permittivity is illustrated and discussed.

Two-dimensional mesh for a rotated cross
The objective is to find a two-dimensional, planar mesh for a rotated cross.Since the procedure how to find such meshes is discussed in great detail in [11], we only briefly review the construction process depicted in Fig. 2. In Fig. 2(a), we show how the unit cell is divided for the mapping.The four points P, Q, R and S define specific coordinate lines (blue and red).The mapping in between these specific coordinate lines is given by a linear interpolation.The resulting mesh is depicted in Fig. 2(c).
One may notice that we could have chosen different specific coordinate lines which would also lead to a grid-aligned cross in the effective permittivity.The reason we choose the ones shown in Fig. 2 becomes clear in the next paragraph.

Constructing the three-dimensional transformation
The lower cross in Fig. 1 is already grid-aligned.This simplifies the discussion at this point but does not impose a restriction.The upper cross can be transformed such that its effective permittivity is also grid-aligned.However, we have to change the coordinate system in between continuously such that artificial reflections are avoided.This means that we obtain a different planar mesh for every value of x 3 .This value of x 3 directly translates into a rotation angle which in turn translates to a mesh like in section 3.1.Explicitly, this means that we perform coordinate transformations in the space between the crosses, too.This explains why we constructed the mesh in Fig. 2 the way we did-for a given x 3 value we only compute the rotation angle and easily obtain the planar mesh.In particular, we hereby make sure that the grid-aligned crosses are directly above each other in the transformed space.
The lower cross is grid-aligned and we choose the origin of our coordinate system to be in the center of the lower surface of this cross.The height of each cross is denoted h and the distance between the crosses is denoted b.Since we want the upper cross to be rotated by the angle ϕ 0 , we obtain for the rotational dependence of the planar mesh on the x 3 coordinate.In Fig. 3, we visualize the real parts of the effective permittivity tensors for two different values of x 3 .The system is a square lattice with lattice constant 600 nm.The cross is 250 nm in diameter and the width of the arms of the crosses is 50 nm.The height of the crosses is h = 25 nm and the spacing between them is b = 50 nm.The angle by which the upper cross is rotated is ϕ = 15 • .We assume the crosses to consist of gold, described by a Drude model with the parameters ε ∞ = 9.0685, a plasma frequency ω D = 1.3544 • 10 16 Hz and a damping coefficient γ = 1.1536 • 10 14 Hz, see [16].The wavelength used for Fig. 3 is 1000 nm.The color scale in Fig. 3(b) has been saturated at 0 for √ gε 11 and √ gε 33 in order to see more features.The real part of the dielectric function of the gold crosses is about − 42.We depict the √ gε 11 , √ gε 12 , √ gε 13 and √ gε 33 components of the effective permittivity.This suffices since the permittivity tensor is symmetric, e.g., √ gε 12 = √ gε 21 and the mapping in Fig. 2 is C 4 symmetric.This results in √ gε 22 and √ gε 23 being √ gε 11 and √ gε 13 rotated counter clockwise by 90 • , respectively.In Fig. 3(a) we display the effective permittivity at x 3 = h + b/2, i.e., between the crosses.As shown, the effective permittivity of this layer of air becomes fully anisotropic due to Eqs. ( 1), ( 2) and (7).In Fig. 3(b) we display the effective permittivity at x 3 = h + b, i.e., in the layer with the rotated cross.Due to the form of the coordinate transformation in Fig. 2, the gold cross is gridaligned in the transformed space.The origin of the discontinuities in the effective permittivity is the fact that the meshes above are not differentiable.This however does not affect the overall performance of the method as long as the discretization parameters are chosen wisely, see [11].

Three-dimensional adaptive spatial resolution
Up to this point we created an effective permittivity where the gold crosses are grid-aligned in both layers.Since we designed the three-dimensional adaptive coordinates such that the crosses are right above each other, it would suffice to apply the same two-dimensional, planar adaptive spatial resolution (ASR) transformation function in each layer.This means that the coordinate lines are compressed at the vicinity of the metallic surface.Mathematically, this is just another coordinate transformation that transforms the effective permittivity.The design of such a transformation is described in detail in [8,11].However, this requires an FMM code that is capable to process coordinate transformations directly, including a switch in the basis.Also, it is difficult to start a plane wave in physical space, since this plane wave is also transformed when we change the basis functions due to the coordinate transformation we perform in every Fig. 3. Real parts of the transformed effective permittivity tensors and the corresponding meshes.x 1 is displayed on the horizontal axis and x 2 on the vertical.Panel (a) displays the real part of the effective permittivity of the test system at x 3 = h + b/2, i.e., in the middle between the crosses.Panel (b) depicts the real part of the effective permittivity at x 3 = h + b, i.e., in the layer of the rotated cross.All results were obtained with a real space discretization of 1024 × 1024 points.The color scale in panel (b) has been saturated at 0 for √ gε 11 and √ gε 33 .Here, the real part of the permittivity of the crosses reaches values of about − 42.layer.Therefore, a three-dimensional ASR is desirable to avoid these problems.
As indicated above, the aim of this paper is to demonstrate a way to incorporate threedimensional coordinate transformation in any classic FMM code.To do so, the basic idea is switching on the adaptive spatial resolution smoothly with increasing x 3 -coordinate.Thereby, the basis functions of the problem represent the real, physical space.Therefore, we can easily start an ordinary plane wave in the incoming, Cartesian, physical half-space.Then, we can introduce several intermediate layers to start the adaptive spatial resolution.The general procedure is sketched in Fig. 4(a).In this sketch, we start by a Cartesian layer at the bottom, cf.Fig. 4(b).In the next three layers, we gradually introduce the ASR as discussed in great detail in [11]."1/3 ASR" figuratively means that the ASR has reached a third of its desired strength, see Figs. 4(c)-4(e).
The layer in which the ASR is fully introduced at its desired strength (cf.Fig. 4(e)) is the layer where the first cross is located (shaded in yellow).Then, like before, we rotate the mesh, only with an ASR applied before that, see Fig. 4(f).Once the mesh is rotated up to the cross rotation angle ϕ 0 , we can compute the layer of the rotated cross, again shaded in yellow in Fig. 4. Three-dimensional meshing with a Cartesian layer as first and last layer.Starting from Cartesian, physical, untransformed space, the density of coordinate lines is gradually increased using a x 3 -dependent compression function.Then, after the layer with the gridaligned cross, the mesh is rotated up to the angle ϕ 0 in the layer with the second cross.After this, the mesh is rotated back and the coordinate line density is reduced until the mesh is Cartesian again.The corresponding meshes are depicted in panels (b) to (g).From those meshes the transformed permittivity and permeability tensors can be computed.Fig. 4(a).The mesh that is used to compute the effective permittivity in this layer is depicted in Fig. 4(g).Like above, we then gradually reverse the mesh changes-first, the mesh is rotated back, then the ASR is decreased until we reach the outgoing layer with a Cartesian mesh.
Mathematically, this looks the following: when we want the increase of the density to happen on the interval x 3 ∈ [0, c], then the mapping has to obey x Here, ASR denotes the compression function (see x(x) in Section 8 in [11]).A linear introduction of the ASR seems most reasonable.Therefore, a suitable function fulfilling the require-ments is The x 2 mapping is constructed similarly.Conceptually, the whole coordinate transformation still has the form of Eqs. ( 1)-(3).For any given value of x 3 we first compress the coordinate lines and then apply the adaptive coordinate transformation.The result is meshes like in Fig. 4. The great advantage of such a procedure is that it can be easily incorporated into any classical FMM code which can solve the large eigenproblem.Since the incoming and outgoing layer are Cartesian, this is perfectly compatible.So any classical FMM code that can solve the large eigenproblem can just be given the transformed permittivity and permeability and can, thereby, incorporate three-dimensional coordinate transformations.Moreover, the issues of two-dimensional transformations that were discussed above are avoided.

Conclusion
This work dealt with the enhancement of the Fourier Modal Method towards three-dimensional adaptive coordinate transformations.We demonstrated how a three-dimensional mesh can be constructed and how this transformation translates into fully anisotropic effective permittivity and permeability tensors.The presented approach can be used to extend any classical FMM code that is able to solve the large eigenproblem such that it uses coordinate transformations, namely by simply transforming the tensors in the presented fashion.While it increases the number of layers to be solved, the overall Fourier representation of the entire structure is highly improved since the structures are transformed to be grid-aligned and the transition in between the structures is performed continuously.This concept expands the range of possibilities for the FMM greatly, especially for complex systems which vary in propagation direction.
Figure 2(b) depicts how these coordinate lines are mapped.

Fig. 2 .
Fig.2.Illustration of the construction principle for a mesh of a cross with given rotation angle.The points P, Q, R and S in panel (a) divide the unit cell into several zones which are mapped according to panel (b).The mapping in between the specific coordinate lines (blue and red) is obtained by linear interpolation.The resulting mesh is depicted in panel (c).