On the use of the Fourier modal method for calculation of localized eigenmodes of integrated optical resonators

We propose the generalization of the Fourier modal method aimed at calculating localized eigenmodes of integrated optical resonators. The method is based on constructing the analytic continuation of the structure's scattering matrix and calculating its poles. The method allows one to calculate the complex frequency of the localized mode and the corresponding field distribution. We use the proposed method to calculate the eigenmodes of rectangular dielectric block located on metal surface. We show that excitation of these modes by surface plasmon-polariton (SPP) results in resonant features in the SPP transmission spectrum. The proposed method can be used to design and investigate optical properties of integrated and plasmonic optical devices.


Introduction
In recent years, study of optical properties of resonant diffraction structures has been given considerable attention [1][2][3][4][5][6]. A subwavelength diffraction grating may serve as a simplest example of a periodic resonant diffraction structure (Fig. 1a). In such structures sharp resonant maxima and minima in reflection and transmission spectra are observed. Such resonances are caused by the excitation of the structure eigenmodes [1,5,6]. These modes are described by the complex-valued frequency [1,2,6]. In periodic structure, modes can be either quasiguided modes, which propagate along the periodicity axis [1], or localized modes, which are supported by the ridges or grooves of the structure [4]. The most universal approach to calculating diffraction of light by periodic structures is through the rigorous coupled-wave analysis (RCWA), which is also referred to as the Fourier modal method (FMM) [7]. Eigenmodes of the periodic structures can be calculated using a modified RCWA approach, which is based on calculating poles of the analytic continuation of the S-matrix of the structure [1,6].  (a) Geometry of a periodic structure: an array of ridges on a metal layer (the structure is infinite in the y-direction); (b) geometry of an aperiodic structure: 2D rectangular cavity on a thick metal layer (the structure is infinite in the y-direction); (c) surface plasmon-polariton transmission spectrum of a rectangular cavity (length w = 900 nm, height h = 600 nm, permittivity ε = 5.5) on a silver layer with thickness 200 nm For practical purposes, of great interest are resonances that occur in aperiodic diffraction structures. Diffraction by aperiodic structures can be calculated using the Aperiodic Fourier Modal Method (AFMM) [8,9]. This method analyzes a periodic structure whose adjacent periods are optically isolated. The isolation is achieved either using anisotropic perfectly matched layers (PML) or gradient absorbing layers [8], or through complex coordinate transform [9].
When dealing with aperiodic diffraction structure, resonances similar to those observed in diffraction gratings occur, except for them being exclusively caused by the excitation of localized modes. Calculation of localized modes in aperiodic structures is important when designing and investigating integrated and plasmonic optical components, such as cavities of vertical-cavity surface-emitting lasers (VCSEL) and spasers, light couplers and out-couplers. Methods for calculating localized modes have been proposed in a number of papers [3,[10][11][12]. The most popular approach uses the FDTD-based analysis of time-dependence of the electromagnetic field amplitude [10]. The major drawback of this approach is low accuracy of calculation of modes with low quality factor. Besides, within the FDTD-based approach it is extremely difficult to calculate the field distribution of a particular mode. In paper [11] the authors calculate the modes of laser resonator using the FMM. However, this method allow one to calculate only real-frequency modes. In paper [12] the FMM was reformulated for cylindrical coordinate system, which allowed the authors to calculated the eigenmodes of body-of-revolution cavities.
In this work, based on the aperiodic Fourier modal method, we propose a rigorous approach to calculating modes of aperiodic diffraction structures. The proposed method by rigorously constructing the analytic continuation of the scattering matrix allows one to calculate the complex frequency and field distribution of the modes. Although the periodic and aperiodic Fourier modal methods are similar, the problem of calculating modes in the aperiodic structures is essentially more challenging, involving a number of aspects which are discussed in this paper.

Modes and resonances of periodic diffraction structures
Let us consider diffraction of the electromagnetic wave by a periodic structure (diffraction grating) shown in Fig. 1a. Usually, the incident waves are assumed to be plane waves. However, in the general case, diffraction of incident periodic electromagnetic waves with period equal to that of the considered structure can also be analyzed. Moreover, we will consider waves that are incident from both substrate and superstrate regions.
For periodic structure, according to the Bloch-Floquet theorem, the transmitted and reflected fields can be represented as a Rayleigh expansion, i.e. an expansion into plane waves (propagating and evanescent diffraction orders). The same expansion can be used for the incident field. To solve diffraction problem we need to find amplitudes of scattered plane waves (diffraction orders) from known amplitudes of incident plane waves. The solution to this diffraction problem can be represented as an S-matrix [2,6,13]. The grating's scattering matrix S relates the complex amplitudes of incident waves (Ψ inc ) with the amplitudes of the scattered waves (Ψ scatt ) as where Here, R and T are the vectors composed of the complex amplitudes of the reflected and transmitted diffraction orders, while I 1 and I 2 are the vectors composed of the complex amplitudes of the plane waves incident on the structure from the substrate and superstrate regions. An S-matrix element with indices (i, j) defines the scattering amplitude of the incident wave with number j into the scattered wave with number i. The S-matrix has a physical meaning only at real light frequencies. Let us analyze the analytic continuation of the S-matrix onto a complex-ω plane. Assume that the determinant of matrix S(ω) has a complex pole at ω = ω p . In this case, det S −1 = 0 and there exist nontrivial solutions to the following homogeneous equation: Thus, at ω = ω p there is a nontrivial solution to Maxwell's equations that does not contain incident waves [see Eq. (1)], which means that the structure has an eigenmode at frequency ω p [6]. Now let us assume that there is an incident plane wave; let us consider the corresponding element of the S-matrix, which is a complex transmission/reflection coefficient of the structure. If the mode corresponding to the pole ω p of the S-matrix is excited by the considered incident wave, ω p will also be the pole of the transmission/reflection coefficient T(ω). Then, the following approximate relation will hold for T(ω) [14]: This equation has the meaning of a truncated Laurent series in the vicinity of point ω p . The transmission coefficient representation (3), which holds in the vicinity of the structure's eigenmode frequency, describes resonant features in the transmission spectrum. Thus, the problem of analyzing and understanding the resonant properties of the structure is reduced to calculating the complex frequencies of the structure's eigenmodes. The real part of the complex frequency defines the resonance frequency, while the imaginary part defines the resonance quality factor Q = Re ω p /(−2 Im ω p ).

Modes and resonances of aperiodic diffraction structures
In this section, we analyze an aperiodic diffraction problem. We consider diffraction of a slab waveguide mode by a non-uniformity or resonator, located near the waveguide. In this case, the incident and scattered fields propagate not in a free space but in a medium, containing the waveguide. In particular, Fig. 2b shows the considered plasmonic waveguide with rectangular block on its interface.
Analysis of aperiodic diffraction structures can be reduced to the analysis of periodic structure as follows. The considered aperiodic structure is periodically continued (see Fig. 1b), with the adjacent periods either being separated by perfectly matched layers (PML) (dashed area in Fig. 1b) [8], or via introducing the complex coordinate transform [9]. This results in optical isolation of adjacent periods and, hence, the solutions to the problem of diffraction by a periodically continued structure and an aperiodic structure become identical. Let us note, that instead of a homogeneous substrate (Fig. 1a), one finds periodic medium above and under the structure (Fig. 1b). Therefore, instead of using Rayleigh expansion, the scattered and incident fields should be expanded as a sum of the modes of periodic medium.
By way of illustration, consider a plasmonic mode supported by a thick silver film diffracted by a dielectric step ( Fig. 1b; the structure parameters are given in the figure caption). The transmission spectrum is shown in Fig. 1c. Note that the transmission coefficient is interpreted as the ratio of the complex amplitude of the transmitted plasmonic mode to the amplitude of the incident mode. The transmission spectrum has pronounced resonance features, bringing to the forefront the problem of calculating the eigenmodes (complex eigenfrequencies) of an aperiodic structure. The calculation of the eigenmodes will allow us to explain the resonant features, as well as to find numerical parameters of the resonances, such as its frequency and quality factor.
For calculating poles of the S-matrix (or of transmission coefficient), one needs to be able to calculate the S-matrix (transmission coefficient) at complex values of the frequency. One of the most universal methods for calculating S-matrix is the Fourier modal method [7,15]. The following section deals with basic formulae of the Fourier modal method for the case of real-valued frequency. In Section 5 we will proceed to construct the analytic continuation.

Rigorous coupled-wave analysis
The RCWA is based on the Bloch-Fouquet theorem, which suggests that for periodic structure, solutions to Maxwell's equations can be represented as a quasi-periodic function. Within the RCWA approach, the structure is assumed to consist of layers whose permittivity is z-independent. In this case, in each layer, the field can be decomposed into a Fourier series in terms of variable x, which denotes the direction of structure's periodicity [7,15] (see Fig. 1a, b).
Let the vector Φ i (z) consist of Fourier coefficients of the electromagnetic field's tangential components in the i-th layer (i = 1, . . . , L). The vector Φ i (z) has dimension 4N, where N is the number of Fourier harmonics employed. From Maxwell's equations, Φ i (z) is seen to satisfy the following matrix differential equation: where the matrix A i ∈ C 4N×4N is defined by the i-th layer's geometry and the incident wave parameters (frequency and quasi-wavenumber k x ) [7,15]. The solution to Eq. (4) describes the field in the i-th layer: where W i is a matrix with columns made up of eigenvectors of the matrix A i ; Λ i is a diagonal matrix of the corresponding eigenvalues; and C i is a vector of arbitrary constants. Relation (5) should be considered as an expansion of the i-th layer's field in terms of the layer's eigenmodes propagating along the z-axis. Note that the propagation constant for the j-th mode in the z-direction is derived from the eigennumber λ j (defining the time-dependence as e −iωt , we have k z, j = −iλ j ). The transverse field distribution of the mode is defined by the j-th column of the matrix W i . The elements of the j-th column are the Fourier coefficients of the field expansion in terms of variable x.
To solve Maxwell's equations for a multilayered structure, we need to equate the tangential field components (or, equivalently, their Fourier coefficients) at the boundaries of the layers. Thus, we obtain the following set of equations: where z = z i is the boundary between layers i and i + 1. System (6) contains Φ 0 (z 0 ) and Φ L+1 (z L ), which denote the Fourier coefficients of the field at the upper and under lower boundaries the structure. Let us consider the diffraction of a waveguide mode (or of a plane wave in the case of periodic structure) by the structure. Assume that the structure is illuminated from top by a mode whose field distribution is described by a vector of Fourier coefficients V inc . We assume V inc to be a column of matrix W 0 . Diffraction of the mode by the structure produces a set of scattered (reflected and transmitted) modes. Let matrix V r consist of matrix W 0 columns that define the reflected modes and matrix V t is composed of matrix W L+1 columns that define the transmitted modes. Thus, the field in the superstrate region can be described as superposition of reflected modes and the incident mode: where R is the vector composed of complex amplitudes of the reflected modes. The field in the substrate region is defined as superposition of transmitted modes: where T is the vector composed of complex amplitudes of the transmitted modes. In the current paper we do not present the expressions for the matrices A, V r , V inc , V t , Φ i . The reader can find the general form of these matrices in [7,15].
Equations (6)-(8) form a system of linear equations in unknown amplitudes of reflected and transmitted modes (R and T ). The diffraction problem is solved by solving this system of linear equations. It is worth noting that in practice the numerical instability of the system's solution can be avoided using a variety of techniques [7,13,15,16]. Besides, when calculating the matrices A i a number of special techniques need to be used [15,17].
If there is a homogeneous medium above and under the structure (see Fig. 1a), the modes of the 0-th and (L + 1)-st layers are plane waves (propagating and evanescent diffraction orders). Thus, matrices A 0 , A L+1 take a simple form, and matrices V t , V r , V inc are defined analytically. In this case, expansions (7) and (8) are referred to as the Rayleigh expansions.
In the general case, matrices V t , V r , V inc are built of the corresponding columns of matrix W . Here, the term "corresponding" is understood as follows. Matrix V t is supposed to contain matrix W L+1 columns that describe transmitted modes, matrix V r contains matrix W 0 columns that describe reflected modes, and vector V inc is matrix W 0 column that describes the incident mode. In most cases (e.g. in the case of reciprocal materials), each mode with propagation constant k z has a "matching" (reciprocal) mode with propagation constant −k z . It is worth noting that one of the two modes is incident whereas the other is scattered. Thus, for the diffraction problem to be solved, all the superstrate and substrate modes need to be classified as either incident or scattered ones. The modes can be classified with respect to either the propagation direction (sign of the real part of k z ) or the evanescence direction (sign of the imaginary part of k z ). Classification of modes into the incident and scattered ones is a fairly delicate issue, which has to be addressed with regard for the physical meaning of a diffraction problem.
If the substrate and superstrate (see Fig. 1a) are homogeneous lossless medium, there will be either propagating or evanescent modes (plane waves) above and below the structure. Then, the scattered modes are chosen so that the propagating modes propagate away from the structure, whereas the evanescent ones decay with increasing distance from the structure. For instance, with the time-dependence of e −iωt , the superstrate scattered waves comprise propagating modes (diffraction orders) with Im k z = 0, Re k z > 0, and evanescent modes (diffraction orders) with Re k z = 0, Im k z > 0. Figure 2 shows the values of where d = 5300 nm, for incident and scattered waves with crosses and circles, respectively. Given an aperiodic, rather than homogeneous, medium under the structure (an infinite grating of Fig. 1b), the propagation constants k z are usually complex numbers with nonzero both real and imaginary parts. Note that the direction of mode evanescence defined by the sign of Im k z and that defined by the sign of Re k z may differ. Hence, the problem of determining the mode propagation direction needs a thorough investigation. This problem becomes especially important when calculating localized modes in plasmonic elements, in which case a dielectric, a metal, and anisotropic PML-material are present both in substrate and superstrate regions.
In this case, modes should be classified as incident and scattered as follows. Modes with real k z can, as before, be classified according to the propagation direction. The rest modes (with complex k z ) should be classified according to its evanescence direction [8]. Hence, the reflected modes will be characterized by Im k z = 0, Re k z > 0 and Im k z > 0 (see Fig. 2b). In this case, it is guaranteed that the solution of the diffraction problem gives field that does not increase with increasing the distance from the structure.
Thus, the diffraction problem can only be solved after accurately choosing the reflected, transmitted, and incident modes from the entire set of the substrate/superstrate modes. For example, this allow one to find the mode transmission coefficient as a function of real frequency (T(ω)). In order to understand the mechanisms behind the resonances, in the following section we construct the analytic continuation of function T(ω) onto the complex ω-plane.

Constructing the analytic continuation
In the case of complex frequency, similarly to the RCWA discussed above, to calculate the transmission coefficient we need to calculate eigendecomposition of matrix A for the substrate/superstrate regions. Then, the modes (i.e. the eigenvectors and the corresponding eigenvalues) should be classified into two groups, namely, the incident and scattered ones. To make sure that T(ω) is analytic, the elements of matrices V t , V r , and V inc should be analytic functions of frequency. The key difference from the real-frequency case is that with complex frequencies we are unable to rely upon physical considerations when classifying modes as the incident or scattered ones. Figure 2a shows that the propagation constants for the incident (blue curves with crosses) and scattered (red curves with circles) modes above the structure change once the imaginary part of the frequency appears (Im ω < 0). The plot suggests that once the frequency is complemented with the imaginary part, the scattered (reflected) modes, which used to be propagating (k z ∈ R + ), become exponentially increasing: k z acquires the negative imaginary part. Thus, if the frequency is complex, we are unable to rely upon the mode propagation or evanescence direction for the V t , V r , and V inc matrices to be constructed.
In a simplest case of a homogeneous dielectric substrate/superstrate (Fig. 1a), the field in the substrate/superstrate regions can be represented analytically [7]. In this case, constructing the analytic continuation of the matrices V t , V r , and V inc (and, hence, the analytic continuation of T(ω)) presents a trivial problem [1].
In the general case of a nonhomogeneous periodic medium found above and under the structure (Fig. 1b), the analytic relations for the eigenvectors of the matrix A are unknown. This makes the problem of constructing the analytic continuation much more challenging. The only condition based on which the modes should be classified as incident/scattered (which is equivalent to constructing the above matrices) is formulated as follows. If the frequency is complex, the mode is incident (scattered) if and only if with imaginary part decreasing to zero, the mode of interest becomes incident (scattered) in the sense of the real-frequency definition of the incident (scattered) mode. Thus, the modes can be classified as the incident or scattered ones by, first, classifying them for a real frequency and, then, "tracking" each of the modes while gradually decreasing the (negative) imaginary part of the frequency. In practice, it will suffice to "track" the change of the propagation constant using the following technique. First, the substrate/superstrate modes need to be calculated for real frequency ω 0 = Re ω and classified as incident/scattered modes. Then, the modes need to be calculated for complex frequency ω 1 = Re ω + i Im ω. Each mode at frequency ω 1 is put into one-to-one correspondence with the "nearest" mode at frequency ω 0 . The mode at frequency ω 1 is assumed to be incident (scattered) if the corresponding mode at frequency ω 0 is incident (scattered).
The proposed method allows each mode at frequency ω 1 to be put in correspondence with the mode at frequency ω 0 , with the total change of the propagation constants being minimized. Note that building a one-to-one correspondence of modes is the most time-consuming operation. The following approach is proposed. We build a "proximity" matrix with its elements defined as where k z,i are the propagation constants for the modes at frequency ω 0 andk z, j is the set of the propagation constants at frequency ω 1 . Then, for the cost matrix P, an assignment problem is solved using the Hungarian algorithm [18]. Generally speaking, if for the mode of interest the imaginary part of frequency is large enough, even the use of the proposed algorithm can give an incorrect result because the propagation constants of the modes at frequency ω 1 may be essentially different from those of the modes at frequency ω 0 . In this case, it is advisable to utilize the proposed technique K times, while consecutively increasing the imaginary part of the frequency: ω k = Re ω + i k K Im ω, k = 0 . . . K. The "evolution" curves in Fig. 2 that illustrate the change of the propagation constant have been calculated using the above-described technique with a large value of K being used. To calculate the modes in the following part of the paper we use K = 1 + − Im ω · 10 −13 s , where square brackets denote the floor operation.
Note that in some particular cases the modes can be classified as incident/scattered ones in an approximate way [3], say, based on the sign of Re k z + Im k z . However, such an approximate approach is not always suitable. In particular, the classification appears to be incorrect at large values of the imaginary part of frequency (see Fig. 2b). Besides, an approximate classification of the modes results in incorrect calculation of the scattered field.
In this section, we have discussed constructing the analytic continuation of the transmission coefficient T(ω). The analytic continuation of the S-matrix is built in a similar way.

Rayleigh anomalies
In the previous section, we described a general approach to classification of the modes as the incident and scattered ones. However, one important case has been omitted. If one mode has the propagation constant k z = 0, two modes with k z = 0 are supposed to be present in the set of modes. One of them should be considered as incident, while the other one as scattered. Thus, in this case the incident mode is indiscernible from the scattered one.
For the diffraction grating, (Fig. 1a), this special case corresponds to Rayleigh anomalies. These anomalies are observed at the Rayleigh frequencies of one of diffraction orders. Below the Rayleigh frequency, the diffraction order is evanescent (ik z ∈ R), above the Rayleigh frequency it is propagating (k z ∈ R). Speaking in terms of the analytic continuation, the Rayleigh anomalies are branch points of the transmission function T(ω). Presence of the branch points leads to a number of remarkable features in the transmission spectrum of diffraction gratings [14].
Similar anomalies are also observed in perfectly conducting optical elements [19], in particular, when coupling perfectly conducting slab waveguides or putting inhomogeneities in perfectly conducting slab waveguides. This is because the dispersion relation for modes of a slab waveguide with perfectly conducting interfaces is identical to the dispersion relation for plane waves in the Rayleigh expansion. For integrated optical elements operating in the visible range, dispersion laws for the incident and scattered modes are much more complex. Because of this, analogs of the Rayleigh anomalies can only be observed accidentally, if two complex propagation constants of the incident and scattered waves get coincident (at a real or complex frequency). It seems unlikely that such a situation may occur at practice.

Calculating modes of an aperiodic structure
According to Eq. (2), the structure modes correspond to the poles of the S-matrix and, hence, it is most natural to find them by calculating the poles of the S-matrix of the structure. In the simplest case, the calculation of the S-matrix poles is reduced to solving numerically the following equation: Note that the straightforward solution of Eq. (10) leads to the numerical instability. Because of this, a number of numerically stable methods have been developed for the calculation of poles of the S-matrix [3,6]. Note that some techniques discussed in [6] require to calculate the derivative of the S-matrix, which should be done with caution: when taking the derivative using finite-difference formulae one should check that the corresponding modes in the substrate and superstrate regions are ordered in the same way at different frequencies [3]. The calculation of modes based on calculating the S-matrix poles is the most universal approach. However, with this approach, alongside finding all localized modes of the structure, one may also derive unphysical Bérenger modes localized in the PML-layers, which have no physical meaning. Thus, when calculating modes as poles of the S-matrix, one has to calculate the field distributions for all modes and then discard nonphysical solutions. If we are only interested in modes that directly affect the transmission spectrum (Fig. 1c) it is advisable to calculate modes as poles of the transmission coefficient. In this case, the eigenfrequencies of the modes are derived numerically as complex roots of the following equation: By way of illustration, using the above-described technique, we calculated modes of a rectangular cavity on a metal substrate (Fig. 1b). The metal substrate was made of silver, with its permittivity described by the Lorentz-Drude model [20]. The cavity width was w = 900 nm, height h = 600 nm, permittivity of the cavity material ε = 5.5, and metal layer thickness 200 nm. The modes were calculated for the frequency range Re ω ∈ 1.0 × 10 15 ; 1.6 × 10 15 s −1 . Note that we analyzed only high-Q modes (Im ω ∈ −1.5 × 10 14 ; 0 s −1 ). Frequencies of the modes were derived as roots of Eq. (11). For complex arguments, the transmission coefficient T(ω) was calculated using the above-discussed algorithm based on the Fourier modal methods [7,13,17]. We used N = 75 · 2 + 1 Fourier harmonics during in the following calculations. Figure 3 shows the electromagnetic field distribution of the modes, which was derived using a numerically stable method [11]. Corresponding complex frequencies are given in Table 1. Note that the use of an approximate technique proposed in [3], in which the modes are classified into the incident and scattered ones according to the sign of Re k z + Im k z , gives an incorrect field distribution, which sharply increases outside the cavity region. The transmission spectrum of Fig. 1c is affected by the modes TM 2 and TM 3 . These modes have the highest Q-factor and their excitation produces sharp minima in the transmission spectrum at frequencies ω = 1.2 × 10 15 s −1 and ω = 1.5 × 10 15 s −1 (see Fig. 1c).
From the field distribution, the modes TM 1 , TM 2 , and TM 3 are seen to be rectangular-cavity modes of the first, second, and third orders, respectively. The mode TM 4 is a plasmonic Fabry-Pérot-type mode. The scattered field distribution shows that mode TM 1 is scattered into a wave directed normally to the surface, whereas modes TM 2 and TM 3 are scattered in several directions. Mode TM 4 is mostly scattered into a surface plasmon-polariton. Thus, the analysis of the mode field distribution makes it possible to determine directions (channels) in which the electromagnetic energy is scattered. Using the reciprocity theorem, it can be shown that these are the same possible directions in which the considered mode can be excited [21]. Moreover, the FMM enables the coupling coefficient between the mode and incident plane wave to be numerically calculated. This analysis is important when designing nanophotonic elements to perform high-efficiency excitation of waveguide modes and surface plasmon-polaritons [21]. In particular, the method for calculating modes proposed in this work was applied by the present authors in combination with the reciprocity theorem to investigate the controlled excitation of surface plasmon-polaritons in rectangular cavities made of a magneto-optical material [22]. As we mentioned above, the calculation of the modes based on Eq. (11) enables obtaining only the modes that can be excited by the considered incident field. With a surface plasmon-polariton used in this work as the incident wave, the set of the calculated modes is restricted just to TM-modes. In addition to polarization-related limitations, the proposed approach may impose limitations on the symmetry of the calculated modes [2,5]. To calculate all structure modes within the considered frequency range, we calculated the modes as poles of the S-matrix [6]. Based on this approach, all modes of Fig. 3 and TE-modes of Fig. 4 were derived. The field distribution of the TE-modes shows them to be uncoupled with the surface plasmon-polariton; these modes scatter into free space above the structure. Table 1 shows the complex frequencies, complex wavelengths, and quality factors of all calculated eigenmodes. As a matter of comparison, in parentheses the values obtained using the MEEP package [23] and library Harminv are presented. MEEP and Harminv implement the FDTD approach and the method of paper [10]. The presented table confirms the validity of the proposed method. The relative difference in frequencies does not exceed 7%. Let us note, that the developer of the MEEP 1.3 package suggests low accuracy when calculating the low-quality modes (Q < 50). Besides, the method of paper [10] does not allow one to calculate the field of one particular mode, whereas within the proposed approach the field distribution is easily calculated (see Figs. 3, 4).

Conclusion
Summing up, we have proposed a generalized RCWA-based approach intended to calculate localized modes of integrated optical cavities. The generalization consists in constructing the rigorous analytic continuation of the S-matrix, thus enabling modes in aperiodic structure to be calculated as poles of the S-matrix. It is of particular importance to build accurate analytic continuation when dealing with low-Q modes. The results of numerical simulation of localized modes in the dielectric structures on a metal substrate demonstrated high efficiency of the proposed method.
In this work, the problem of mode calculation has been reduced to calculating poles of the S-matrix. In an alternative approach suggested in Refs. [11,24] the modes were calculated on the assumption of generating the Fabry-Pérot type resonances. Note that methods discussed in Refs. [11,24] can also be generalized using the proposed method for constructing the analytic continuation of the S-matrix.