1. Introduction

Over the last decades, investigation of resonant properties of nanophotonic structures has attracted considerable research attention [1, 2]. This is not only due to fundamental interest, but also due to the potential applications of resonant structures for spectral and spatial filtering [3, 4], enhancing light-matter interaction [2], and analog optical computing [5, 6], to name a few. Resonant effects in such structures are associated with the excitation of eigenmodes. Moreover, it was recently shown that the knowledge of the eigenmodes enables fast and accurate prediction of optical properties of the structure under varying excitation conditions [7, 8]. Therefore, the ability to calculate the eigenmodes of nanophotonic structures is of crucial importance.

One of the most widely studied classes of resonant optical structures are periodic structures, in particular, diffraction gratings or photonic crystal slabs (PCSs) [1, 4]. PCSs are planar structures with 1D or 2D periodicity. For the structures with spatial periodicity, several efficient numerical methods for solving Maxwell’s equations were developed [9–12]. Among the most commonly used and efficient methods is the rigorous coupled-wave analysis (RCWA) also referred to as the Fourier modal method (FMM). In its simplest formulation, the FMM is intended for simulating plane wave diffraction on binary lamellar gratings with 1D periodicity [9]. Over the years, numerous improvements to the method have been proposed, concerning, on the one hand, numerical stability, convergence and computational efficiency [13, 14], and, on the other hand, extension of the applicability of the method to multilayer structures [15, 16], structures with 2D periodicity [17, 18], and anisotropic structures [18].

The FMM also allows one to calculate the eigenmodes of the structure. Within the FMM-based approach, the eigenmodes can be calculated as the poles of the scattering matrix [19, 20] using different numerical methods [7, 20–23]. These methods can be used to calculate complex wave numbers of real-frequency modes, complex frequencies of quasiguided modes (real-wave-number modes), and mode field distributions. The complex-wave-number modes decay in space and allow one to describe resonant features in the angular spectrum [24], whereas the complex-frequency modes decaying in time are widely used to describe resonances in the frequency spectrum [19, 21].

An important milestone in the evolution of the FMM family was the development of the aperiodic FMM (AFMM) [25, 26]. Within the AFMM, the so-called perfectly matched layers (PML) are added to the computational cell of the structure in order to optically isolate the adjacent periods [25–27]. This allows one to apply the method to the solution of integrated optics problems. The recently proposed contrast field formulation of the AFMM enables simulating plane wave diffraction on isolated non-periodic optical structures [28, 29]. It is also worth mentioning the modification of the AFMM for calculating localized complex-frequency eigenmodes of integrated optical resonators [30].

Surprisingly, the aperiodic FMM also enables calculating the eigenmodes of periodic structures. In paper [25], complex wave numbers of the real-frequency modes were found as the eigenvalues of the transfer matrix, which was calculated using the AFMM. A similar approach incorporating the scattering matrix was proposed in papers [31, 32]. However, calculation of complex-frequency modes is substantially more sophisticated and cannot be performed using the methods of Refs. [25, 31, 32].

In the present work, we propose an iterative AFMM-based method that enables calculating complex frequencies and field distributions of quasiguided PCS modes with 1D periodicity. The proposed method, along with the methods of papers [25, 28, 29, 31, 33], is computationally efficient when modeling long-period optical structures. In particular, the method proposed in this work allows one to calculate the eigenmodes of the PCS with a very long period (at least up to 500λ). To the best of our knowledge, the known methods do not allow one to calculate the modes of structures with such a long period. The ability to simulate long-period PCS is particularly useful for the analysis of random and quasiperiodic structures within the super-cell approach [34].

Following the Introduction, Section 2 defines the scattering matrix and gives the general dispersion relation of the PCS eigenmodes within the AFMM approach. In Section 3, we propose an iterative method for calculating complex frequencies of the eigenmodes. Section 4 is dedicated to the numerical validation of the method and its comparison to the known methods. In Section 5, we discuss the details of the application of the proposed method to long-period structures. Section 6 concludes the paper.

2. Eigenmodes and aperiodic Fourier modal method

2.1. The scattering matrix and the eigenmodes

In this paper, we consider periodic optical structures [diffraction gratings or photonic crystal slabs, see Fig. 1(a)]. According to the Bloch–Floquet theorem, for the structure periodic along the x axis, the solutions to Maxwell’s equations satisfy the following quasiperiodicity condition:

 

Fig. 1 (a) Photonic crystal slab geometry. (b) “Rotated” PCS geometry. The adjacent periods are optically isolated by means of PML. Both structures are uniform in the y direction.

Download Full Size | PDF

Using Eq. (1), the field in the substrate and superstrate regions can be represented as Rayleigh expansions. Within this expansion, the field is represented as a sum of incident and scattered plane waves (diffraction orders) with the x components of the wave vector k_x + mG, where G = 2π/d is the magnitude of the reciprocal lattice vector and m ∈ ℤ is the diffraction order number.

The scattering matrix of a periodic structure relates the complex amplitudes of the incident plane waves to the complex amplitudes of the scattered plane waves:

The modes of the structure correspond to the poles of the scattering matrix [20]. Thus, the dispersion relation of the eigenmodes is given by the following equation:

2.2. Eigenmode description in the AFMM-based approach

Let us consider the aperiodic Fourier modal method (AFMM). According to this method, we take one period of the structure and rotate it by 90 degrees [see Fig. 1(b)]. Then, we periodize the structure with respect to the $\tilde{x}$ direction. Finally, we optically isolate the adjacent periods with the so-called perfectly matched layers (PML) [26]. In this paper, we use anisotropic PML [26]; however, one can also use gradient absorber [26] or PML based on complex coordinate transform [27].

In order to describe a mode of the periodic structure shown in Fig. 1(a) within our AFMM-based approach, we need to match the field components at $\tilde{z}=0$ and $\tilde{z}=d$ [see Fig. 1(b)] according to the quasi-periodicity condition (1) with the wave number k_x. We can do this by matching the complex amplitudes of the upward (downward) propagating modes:

Assuming ω = ω_p is a mode frequency, we can use Eq. (4) to rewrite Eq. (2) in the following form:

Equation (5) defines the dispersion of the modes of the PCS shown in Fig. 1(a). This equation is the AFMM version of Eq. (3). Let us remind that S in Eq. (3) is the scattering matrix of the structure shown in Fig. 1(a), while $\tilde{\mathbf{S}}$ in Eq. (5) is the scattering matrix of the “rotated” structure shown in Fig. 1(b).

Let us now discuss the choice of the substrate and superstrate of the “rotated” structure shown in Fig. 1(b). In the problem considered in this paper, the geometry of these regions only determines the basis, in which Eq. (4) matches the fields at the upper $(\tilde{z}=d)$ and lower $(\tilde{z}=0)$ boundaries of the “rotated” structure. Therefore, the geometry of substrate and superstrate can be chosen arbitrarily, provided that it is the same for both regions. In the conventional AFMM, the substrate and superstrate regions contain PML layers [26]. In the proposed method, however, it is useful to consider uniform substrate and superstrate as it is shown in Fig. 1(b). In this case, instead of solving an eigenvalue problem in these regions, we use simple Rayleigh plane wave expansion [9] and a conventional RCWA/FMM implementation [9, 13, 15] for calculating $\tilde{\mathbf{S}}$.

3. The iterative method

3.1. Formulation

In this section, we propose an iterative method that allows one to solve Eq. (5) with respect to complex ω_p, assuming that the mode wave number k_x is real. Note that the case of modes with complex wave number and real frequency is substantially simpler; it was considered in papers [25, 31, 32].

Assume that we are given a certain approximation of the mode frequency ω_n. We use the AFMM to calculate the scattering matrix $\tilde{\mathbf{S}}({\omega }_n)$ of the “rotated” structure. Then, by calculating the finite-difference derivative approximation of the scattering matrix, we expand it into Taylor series in the vicinity of ω_n:

Assuming ω = ω_p is the eigenfrequency of a mode of the structure, we multiply the last equation by $\tilde{\mathbf{\psi }}$ and, according to Eq. (5), we obtain

When the next iteration of the method does not change the value of ω_n (i. e. ω_n+1 = ω_n), Eq. (7) becomes Eq. (5). Hence, the found frequency corresponds to a mode of the structure. In practice, the iterative process is repeated while |ω_n+1 − ω_n| exceeds some defined tolerance.

Let us note that the proposed iterative method (8) is similar to the iterative method presented in papers [21, 35]. In these works, the method was based on linearizing the inverse of the scattering matrix, S⁻¹, of the structure shown in Fig. 1(a), while the proposed method (8) was obtained by linearizing the scattering matrix $\tilde{\mathbf{S}}$ of the “rotated” structure shown in Fig. 1(b).

3.2. Convergence

The proposed method is a Newton-type iterative method. Unfortunately, the methods of Newton type sometimes diverge. For example, if the matrix $\tilde{\mathbf{S}}(\omega )$ has a pole at some frequency ω near the initial approximation, the Taylor series in Eq. (6) has a small radius of convergence. This may result in divergence of the iterative method. Let us note that the poles of $\tilde{\mathbf{S}}(\omega )$ do not correspond to the modes of the considered structure of Fig. 1(a). The matrix $\tilde{\mathbf{S}}(\omega )$ is likely to have poles when the period of the grating d is very long. Besides, the position of the poles (and, hence, the convergence of the method) depend on the used substrate/superstrate geometry [see Fig. 1(b)].

A simple way to obtain convergence in this case consists in using the so-called damped Newton method. To do this, we multiply the minimal eigenvalue in Eq. (8) by the damping factor α, with 0 < α < 1. This, as a rule, ensures the convergence of the method, albeit with a lower convergence rate than in the non-damped method.

Another factor affecting the convergence is the choice of the unit cell. Indeed, for the structure with period d we can arbitrarily choose the borders of the unit cell in the following way: [x, x + d], where x ∈ [0, d]. If the structure supports several modes with frequencies located near the initial approximation ω₀, different choices of the unit cell may result in the convergence of the method to different modes. The simulations show that the method usually converges to the nearest mode that has non-zero field distribution at the border of the unit cell.

3.3. Mode field

Within the AFMM approach, the field distribution of the mode can be calculated using the method presented in paper [32], which remains numerically stable when the light frequency is a complex number.

For the description of the far-field behavior of the mode, one can calculate the scattered diffraction orders amplitudes [R and T in Eq. (2)] in the following way. Firstly, we calculate the field distribution along the boundaries of the PML. Then, we expand the field into the Fourier series to calculate the reflected diffraction orders amplitudes:

3.4. Bérenger modes

Along with the “physical” modes, the proposed method can sometimes converge to non-physical modes localized inside the PML. These modes are sometimes referred to as Bérenger modes. One approach to identify a Bérenger mode is to calculate and analyze its field distribution, since such a mode will not decay towards the depth of the PML.

Let us consider an alternative approach to identifying Bérenger modes that requires neither the calculation of the field distribution nor its human-aided analysis. As mentioned in Subsection 2.2, solving Maxwell’s equations using the AFMM requires one to formally define the wave number $k_{\tilde{x}}$ with respect to the $\tilde{x}$ direction. This wave number gives the phase shift between the left and right boundaries of the period $\tilde{d}$ shown in Fig. 1(b). If the field obtained via AFMM decays inside the PML, the field distribution will not change with changing $k_{\tilde{x}}$. However, for the Bérenger mode, which is localized inside the PML, the change in $k_{\tilde{x}}$ will result in a substantial change of the mode field distribution and the mode frequency. Thus, after calculating the mode frequency ω_p using the proposed iterative approach, one can change the value of $k_{\tilde{x}}$ and run the iterative method again with the frequency ω_p taken as the initial approximation. If the method converges to a frequency which is very close to ω_p, then the mode in question is a “physical” mode of the structure; otherwise, it is a Bérenger mode.

4. Numerical validation

In this section, we study the performance of the proposed method. We compare the new method with a known iterative method presented in paper [20], which is based on the calculation of the poles of the scattering matrix S. The comparison is performed in terms of the computation time required to calculate the complex frequency of the mode with a required accuracy. We define the accuracy by the relative tolerance value δ = |ω_n − ω_p|/|ω_p|, where ω_p is the exact frequency of the mode and ω_n is the numerically calculated frequency obtained using one of the considered methods.

With the decrease in δ, firstly, the required number of harmonics N in the Fourier field expansions increases. Secondly, a larger number of iterations n might be required for the methods to converge. Let us note that for a given tolerance δ, the N and n values are, in general, different for the two considered methods. Besides, the performance of the proposed method depends on the parameters of the PML layer, which are different for different values of δ. The details of the performed numerical experiments are presented in the Appendix.

In the rest of this section, we consider two different examples. Subsection 4.1 is dedicated to the study a subwavelength PCS, which is considered as a structure with a super-period. In Subsection 4.2, we study long-period PCSs, which are perturbations of the subwavelength PCS.

4.1. Subwavelength PCS

We start our analysis with a subwavelength PCS (periodically perforated slab waveguide) with the following parameters: period d₀ = 1000 nm, height h = 500 nm, slit width w = 200 nm (fill factor 80%), dielectric permittivity ε = 2, surrounding medium permittivity ε = 1. We will denote this structure by ${\mathcal{S}}_1$. The calculations using the method of paper [20] show that the considered structure supports a TM-polarized mode with the frequency ω_p = (1.615628−0.002594i)·10¹⁵ s⁻¹ and the wave number k_x = 0 (the details of the numerical experiment are given in the Appendix). In order to investigate the performance of the method for long-period structures, we will study the PCS ${\mathcal{S}}_m$ that has the same geometry as ${\mathcal{S}}_1$, but is considered as the structure with the (super-)period d = md₀.

Numerical simulations show that the use of iterative method from paper [20], which is based on solving Eq. (3), requires approximately the following number of Fourier harmonics: N = 2·2m + 1 for relative tolerance δ = 10⁻³, N = 2·3m + 1 for δ = 10⁻⁴, and N = 2·8m + 1 for δ = 10⁻⁵. Since FMM requires O(N³) operations, the total computation time grows cubically with increasing the period [36]. As for the memory, FMM operates with matrices of O(N²) size; hence, the memory consumption grows quadratically with increasing the period [9].

In contrast, the required number of harmonics in the proposed AFMM-based method does not depend on the period of the investigated structure due to the alternative discretization direction (rotation of the structure) [36]. Obviously, N still depends on the tolerance δ; numerical simulations give the following values: N = 2·4 + 1 for δ = 10⁻³, N = 2·13 + 1 for δ = 10⁻⁴, and N = 2·30 + 1 for δ = 10⁻⁵. However, the number of steps in the S-matrix algorithm [15] used for the calculation of $\tilde{\mathbf{S}}$ depends on the m value defining the number of layers of the “rotated” structure. This allows us to expect approximately linear growth in the computation time with increasing the period [36]. At the same time, the memory requirements for the proposed method do not depend on the period.

Figure 2 shows the computation time vs. m for the two compared methods. For the method of paper [20], the computation time grows smoothly with increasing the period d = md₀. It is evident that at small m, this method is faster comparing to the proposed one. At the same time, the lines corresponding to the proposed method exhibit pronounced “jumps”, each occurring due to an increase in the number of iterations required to calculate the mode frequency with a given tolerance δ. These “jumps” violate the expected linear dependency of the computation time on the period. Nevertheless, at m > 12 the proposed method turns out to be more efficient for all considered δ values. For the values of m as large as 40, the proposed method is several times faster than the known method.

 

Fig. 2 Computation time vs. the super-period size m for the proposed method (dashed lines) and the method of paper [20] (solid lines). Different lines correspond to different relative tolerances δ.

Download Full Size | PDF

4.2. Long-period PCS

Now let us consider the calculation of the modes of the structure ${\overline{\mathcal{S}}}_m$ shown in Fig. 1. In this structure, the height of each m-th rod is halved. We will refer to the structure ${\overline{\mathcal{S}}}_m$ as the perturbation of the structure ${\mathcal{S}}_m$. Let us note that d = md₀ is the minimal period of the perturbed structure ${\overline{\mathcal{S}}}_m$, while d₀ is no longer its period.

Table 1 shows the complex frequencies of the modes calculated using the proposed method (ω_n) and the method of paper [20] $({\omega }_n^{\prime })$ along with the required computation time. The results are presented for the structures with various periods d = md₀ and for different relative tolerances δ. Let us note that the obtained values of ω_n and ${\omega }_n^{\prime }$ satisfy the relation $|{\omega }_n-{\omega }_n^{\prime }|\le \delta \left|{\omega }_n\right|$ which confirms the validity of the proposed method.

Table 1. Performance of the proposed method in comparison to the method of paper [20]

View Table

According to Table 1, for relatively short-period PCSs (m = 5, m = 10) the considered methods require comparable amounts of time. For the case of m = 50, the proposed method is by an order-of-magnitude faster comparing to the known method. To calculate the modes of the structures with the period as large as 100d₀, the use of the proposed method is the only option, since the methods based on solving Eq. (3) require enormous amounts of memory and computation time. Therefore, the proposed method allows one to efficiently calculate the modes of long-period PCSs.

5. Mode dispersion and mode field analysis

As it was demonstrated in the previous section, the proposed method is useful for calculating modes of long-period PCSs. In this section, we will consider the issues arising when calculating modes of a PCS with a very long period. As a matter of illustration, we will study long-period PCSs that are perturbations of a subwavelength PCS. In Subsection 5.1, we discuss how the dispersion law ω_p(k_x) of the modes changes when the period of the PCS increases. In Subsection 5.2, we investigate the field distribution of the modes.

5.1. Mode dispersion

Using the notations of Section 4, let us first consider the subwavelength structure ${\mathcal{S}}_1$ with the period d₀ < λ. The parameters of this structure are given above in Section 4. The dispersion curves of the modes of ${\mathcal{S}}_1$ are shown with solid lines in Fig. 3(a). Following Subsection 4.1, we consider the same structure, but interpret it as a structure ${\mathcal{S}}_m$ with the period d = md₀. The dispersion curves of ${\mathcal{S}}_m$ can be obtained by repeating the dispersion curves of ${\mathcal{S}}_1$ with a shift along the k_x axis by $\frac{2\pi }{m{d}_0}\cdot p$, where p is an integer. As a result, the dispersion curves of ${\mathcal{S}}_1$ will fold to the first Brillouin zone of ${\mathcal{S}}_m$. As a matter of example, Fig. 3(a) shows the dispersion curves of ${\mathcal{S}}_5$ with dotted lines.

 

Fig. 3 Mode dispersion of the structures ${\mathcal{S}}_1$ (a) (solid lines), ${\mathcal{S}}_5$ (a) (solid and dashed lines), and ${\overline{\mathcal{S}}}_5$ (b).

Download Full Size | PDF

Now let us consider the structure ${\overline{\mathcal{S}}}_m$, which is a perturbation of the structure ${\mathcal{S}}_m$. By perturbation we mean a small local change in the geometry or the refractive index of the structure (see Subsection 4.2). We will assume that the perturbed structure ${\overline{\mathcal{S}}}_m$ has the period d = md₀, while d₀ is no longer its period. The dispersion curves of the perturbed structure ${\overline{\mathcal{S}}}_m$ are somewhat similar to ${\mathcal{S}}_m$, and so are the modes. As a matter of example, Fig. 3(b) shows the dispersion of the modes of ${\overline{\mathcal{S}}}_5$ [the geometry of the perturbed structure is shown in Fig. 1(a)].

According to Fig. 3(a), at k_x = 0 the structure ${\mathcal{S}}_1$ has two modes: the even mode, A₂, and the odd one, A₁ [37]. With the increase in the m value, the first Brillouin zone $\left[-\frac{\pi }{m{d}_0},\frac{\pi }{m{d}_0}\right]$ “collapses”. Hence, being folded to the first Brillouin zone, the dispersion curve of the structure intersects the line k_x = 0 many times. Therefore, the structure ${\overline{\mathcal{S}}}_m$ has several modes at k_x = 0. These modes are denoted by A_1,2, B_1,2, C_1,2, D_1,2, and E_1,2 in Fig. 3(b).

As it was noted in the previous section, in order to calculate the mode, one fixes the wave number k_x, and for a given mode frequency approximation ω₀ one uses the iterative method to calculate the nearest complex eigenfrequency. However, according to Fig. 3(b), when the period of the structure is long, the frequencies of the modes become close to each other. Hence, the iterative method can converge to any of the modes in the vicinity of the given approximation.

Let us now investigate whether these modes can be efficiently excited, i. e. whether they result in resonances in the transmission (or reflection) spectrum. Dashed line in Fig. 4 shows the transmission spectrum of the structure ${\mathcal{S}}_1$ for the case of a normally incident plane wave. The spectrum has a pronounced resonance, which corresponds to the excitation of the even mode A₂. Needless to say, the structure denoted as ${\mathcal{S}}_m$ has the same transmission spectrum. However, the transmission spectrum of the perturbed structure ${\mathcal{S}}_{20}$, shown with a solid red line, is a bit different. It has many low-amplitude resonances, which correspond to the intersections of the line k_x = 0 with the folded dispersion curves [similar to the points B_1,2, C_1,2, D_1,2, and E_1,2 shown in Fig. 3(b) corresponding to the structure ${\overline{\mathcal{S}}}_5$]. However, only the mode A₂ (calculated in Section 4), which is the perturbation of the mode of the non-perturbed structure ${\mathcal{S}}_1$, corresponds to a pronounced resonant dip. Therefore, of particular interest is the calculation of complex frequencies and field distributions of the modes that are the perturbations of the modes of the non-perturbed structure. One way to check whether a mode of this kind was found is to analyze its field distribution.

 

Fig. 4 Transmission spectra of the structures ${\mathcal{S}}_1$ (dashed line) and ${\overline{\mathcal{S}}}_{20}$ (solid line).

Download Full Size | PDF

5.2. Mode field

If the structure ${\mathcal{S}}_1$ has low optical contrast, the modes of ${\mathcal{S}}_1$ can be considered as perturbations of the modes of a slab waveguide [37]. Inside the waveguide, the mode field can be represented as a pair of propagating plane waves. Outside the waveguide, the mode field is an evanescent wave that has the same in-plane wave vector component k_x as the waves inside the waveguide. The value of k_x is referred to as the wave number of the mode.

Now let us consider a mode of the PCS ${\mathcal{S}}_1$, which is a periodically perforated slab waveguide. Being similar to the field of a slab waveguide mode, the field of the PCS mode contains high-amplitude evanescent waves in the scattered field. This means that the vectors R and T in Eq. (2) contain a large number at some position. This position is sometimes referred to as the number of the diffraction order that excites the mode. As a matter of example, the modes lying on the dispersion curve containing A₁ at k_x > 0 [see Fig. 3(a)] have strong −1-st diffraction order, while the modes on the dispersion curve containing A₁ at k_x < 0 have high-amplitude diffraction order with the number +1.

Let us return to the structure ${\mathcal{S}}_m$, which is the structure ${\mathcal{S}}_1$ considered as a PCS with a super-period md₀. As it was noted in Subsection 5.1, the dispersion curves of ${\mathcal{S}}_m$ are the shifted dispersion curves of ${\mathcal{S}}_1$. The shift of the dispersion curve corresponds to the shift in the vector of the diffraction order amplitudes. Hence, the modes B–E [see Fig. 3(b)] have substantially different scattered field amplitudes as compared to the modes A_1,2, which have maximal ±m-th diffraction orders. In the case of the mode B, the maximal-amplitude diffraction orders have the numbers ±(m + 1). Hence, by analyzing the amplitudes of the diffraction orders, one can distinguish the modes of a long-period PCS.

The same approach can be used to analyze the modes of the perturbed structure ${\overline{\mathcal{S}}}_m$. For the modes, which were calculated using the iterative approach presented in Section 3, one can use Eq. (9) to calculate the amplitudes of the scattered diffraction orders R and T. By finding the maximal-amplitude diffraction order number, one can identify which of the modes marked in Fig. 3(b) was found. If the found mode was not the mode of interest, one can correspondingly increase or decrease the value of the initial approximation ω₀ and run the iterative method again.

6. Conclusion

Within the framework of the aperiodic Fourier modal method, we presented a new iterative approach, which allows one to calculate complex frequencies of quasiguided modes of periodic optical structures. The method is particularly useful for investigating long-period diffraction gratings and photonic crystal slabs. Possible applications include investigating the impact of defects on resonant properties of photonic crystal slabs, and studying optical properties of random and quasiperiodic photonic structures. We have shown that with the increase in the period of the structure, the number of eigenmodes increases as well. To distinguish between these modes, one needs to analyze its field distribution.