1. 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–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–8]. In recent years, different concepts emerged for the construction of the corresponding meshes [9–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.

2. 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̄³-direction and periodic with respect to the x̄¹x̄² plane. k_in denotes the wave vector of an incident plane wave. Ox̄¹x̄²x̄³ describes a Cartesian coordinate system.

Download Full Size | PDF

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–11, 15]. We distinguish between a curvilinear coordinate system Ox¹x²x³ and a Cartesian coordinate system Ox̄¹x̄²x̄³. The three-dimensional adaptive coordinate transformations that are investigated in this paper have the form

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

3.1. 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). Figure 2(b) depicts how these coordinate lines are mapped. 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.

 

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

Download Full Size | PDF

3.2. 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³. This value of x³ 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³ 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 φ₀, we obtain

In Fig. 3, we visualize the real parts of the effective permittivity tensors for two different values of x³. 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¹⁶ Hz and a damping coefficient γ = 1.1536 · 10¹⁴ 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 $\sqrt{g}{\epsilon }^{11}$ and $\sqrt{g}{\epsilon }^{33}$ in order to see more features. The real part of the dielectric function of the gold crosses is about −42. We depict the $\sqrt{g}{\epsilon }^{11}$, $\sqrt{g}{\epsilon }^{12}$, $\sqrt{g}{\epsilon }^{13}$ and $\sqrt{g}{\epsilon }^{33}$ components of the effective permittivity. This suffices since the permittivity tensor is symmetric, e.g., $\sqrt{g}{\epsilon }^{12}=\sqrt{g}{\epsilon }^{21}$ and the mapping in Fig. 2 is C₄ symmetric. This results in $\sqrt{g}{\epsilon }^{22}$ and $\sqrt{g}{\epsilon }^{23}$ being $\sqrt{g}{\epsilon }^{11}$ and $\sqrt{g}{\epsilon }^{13}$ rotated counter clockwise by 90°, respectively.

 

Fig. 3 Real parts of the transformed effective permittivity tensors and the corresponding meshes. x¹ is displayed on the horizontal axis and x² on the vertical. Panel (a) displays the real part of the effective permittivity of the test system at x³ = h + b/2, i.e., in the middle between the crosses. Panel (b) depicts the real part of the effective permittivity at x³ = 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 $\sqrt{g}{\epsilon }^{11}$ and $\sqrt{g}{\epsilon }^{33}$. Here, the real part of the permittivity of the crosses reaches values of about −42.

Download Full Size | PDF

In Fig. 3(a) we display the effective permittivity at x³ = 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³ = 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 grid-aligned 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].

4. 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 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 three-dimensional coordinate transformation in any classic FMM code. To do so, the basic idea is switching on the adaptive spatial resolution smoothly with increasing x³-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).

 

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³-dependent compression function. Then, after the layer with the grid-aligned cross, the mesh is rotated up to the angle φ₀ 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.

Download Full Size | PDF

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 φ₀, we can compute the layer of the rotated cross, again shaded in yellow in 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³ ∈ [0, c], then the mapping has to obey

Conceptually, the whole coordinate transformation still has the form of Eqs. (1)–(3). For any given value of x³ 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.

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