Resonant state expansion applied to planar waveguides

The resonant state expansion, a recently developed method in electrodynamics, is generalized here to planar open optical systems with non-normal incidence of light. The method is illustrated and verified on exactly solvable examples, such as a dielectric slab and a Bragg reflector microcavity, for which explicit analytic formulas are developed. This comparison demonstrates the accuracy and convergence of the method. Interestingly, the spectral analysis of a dielectric slab in terms of resonant states reveals an influence of waveguide modes in the transmission. These modes, which on resonance do not couple to external light, surprisingly do couple to external light for off-resonant excitation.


I. INTRODUCTION
Optical waveguides (WGs) are a basic building block for optical technology owing to their lossless guiding of light, enabled for example by total internal reflection. WGs provide confinement of the light in one or two dimensions, while allowing light waves to propagate along the remaining dimensions in which the waveguides are approximately invariant. Planar WGs with onedimensional (1D) confinement, such as a dielectric slab, and fiber WGs with two-dimensional (2D) confinement, are widely used, for example in fiber optic cables for telecommunication, photonic crystal fibers, 1 integrated optical circuits 2 , and terabit chip-to-chip interconnects. 3 The optical spectra of WGs, however, do not consist of only these bound modes called WG modes, but also contain unbound modes which couple to the outside, commonly known as leaky modes. An elegant and intuitive way to understand and describe the properties of optical systems is to use the concept of discrete resonant states 4,5 (RSs) which include all types of modes in the system and present a mathematically complete set of spatial functions. RSs are defined as eigensolutions of the Maxwell equation having outgoing wave boundary conditions. Their energies are generally complex reflecting the fact that these states decay in time and leak out of the system. RSs with complex energies are therefore characterized by exponentially growing tails outside the system that require an adapted normalization. 4-8 WG modes, which may exist in the system, are included in the set of RSs and are required for the completeness of the set, even though they have real energies and evanescent tails.
To calculate the RSs in optical systems in which analytic solutions are not possible, the resonant state expansion (RSE), a rigorous perturbation method in electrodynamics, has been recently developed 6 and applied to finite 1D and 2D systems, such as planar 7 and cylindrical 8 resonators. The RSE was shown to be particularly suited for the calculation of sharp resonances, such as whispering gallery modes in microcylinders [9][10][11] and microspheres, 12 where available computational techniques [13][14][15][16][17] need prohibitively large computational resources. 18,19 Up to now, the RSE has been applied only to modes with zero wave vector p along the translationally invariant direction of the system considered, corresponding to normal incidence of light without propagation along the waveguide. In this work, we extend the application of RSE to arbitrary wave vectors p, thus allowing to describe the propagation along waveguide structures. This introduces in the spectrum of RSs, which for normal incidence is dominated by lossy Fabry-Perot (FP) modes, WG and anti-waveguide (AWG) modes, as well as a continuum of modes due to a cut of the Green's function in the complex frequency plane appearing for p = 0. The modes on the cut contribute significantly to the optical spectra and are required for the completeness of the RS basis. They present a challenge in the technical implementation of the RSE which is dealing with discrete states. We have recently shown 8 that one can make an effective discretization of such continua for the RSE applied to 2D systems which show a cut already for p = 0. In the present work, we eliminate the cut in planar systems with p = 0 by going from the frequency representation of the system to the normal wave-vector representation.
We treat here planar WGs, while the application of the RSE to fiber WGs, generalizing our recent work on cylindrical resonators 8 to non-normal incidence, will be the subject of a future work. We verify our theory on exactly solvable structures such as a homogeneous dielectric slab and a Bragg-mirror microcavity, using the RSs of a reference slab as a basis for the RSE. The role of the different types of RSs is studied in detail, revealing the importance of WG modes in the transmission.
The paper is organized as follows. In Sec. II we study the transmission of a homogeneous slab in the complex frequency and normal wave vector plane, in order to analyze the contributions of different types of RSs to the optical spectra of planar WGs. In Sec. III we present a general formulation of the RSE for planar systems with non-zero in-plane momentum. In Sec. IV we demonstrate applications of the RSE to different systems and compare results with available exact solutions. In particular, we introduce in Sec. IV A the basis of RSs for a homogeneous slab in inclined geometry and then use it for calculation of optical modes of a homogeneous slab with a different refractive index in Sec. IV B and of a Bragg-mirror microcavity in Sec. IV C.

II. ROLE OF WAVEGUIDE MODES IN TRANSMISSION SPECTRA
We study the role of RSs in the transmission of a dielectric slab, and in particular the influence of the WG modes on the slab transmission. The WG modes are RSs which have zero linewidth and are present in the spectrum of a planar system at non-normal incidence of the incoming light wave. We consider a dielectric slab with thickness 2a in vacuum, having the real dielectric constant where s is the permittivity of the slab and z is the coordinate normal to the slab. We assume a permeability of µ = 1 everywhere throughout this work. The electric field E satisfies Maxwell's equation, and Maxwell's boundary conditions on the dielectric/vacuum interfaces. For an incoming plane monochromatic wave with the transverse-electric (TE) polarization alongŷ (ŷ is the unit vector along the yaxis) and real frequency ω, the electric field in the system takes the form in which p is the in-plane projection of the wave vector. For the component E(z) of the electric field, Eq. (2) transforms to a 1D wave equation The electric field for z > a is given by the transmitted plane wave E(z) = t(ω)e ikz E 0 where E 0 is the amplitude of the incoming wave. The field transmission through the slab t(ω) has the analytic form t(ω) = 2ikqe 2ika 2ikq cos(2qa) + (k 2 + q 2 ) sin(2qa) = T (k), (5) in which are the z-components of the wave vector in vacuum and dielectric, respectively. Eq. (5) shows that the transmission t(ω) is a function of the real frequency ω. One can also express the transmission t(ω) as a function T (k) of the normal wave vector k, in which k takes only real positive values, as dictated by the outgoing character of the transmitted wave. The wave vector q inside the slab can be complex for a dielectric with dissipation and have an arbitrary sign, reflecting the fact that waves within the slab propagate in both directions. Hence the transmission is insensitive to the sign of q as seen in Eq. (5).
To study the influence of different modes on the transmission, we consider analytic continuations (ACs)t(ω) andT (k) of both functions in the complex ω-and kplane, respectively, in order to investigate their pole structure and for each of them apply the Mittag-Leffler theorem. [20][21][22] The AC of the transmission has different types of poles, which are shown in Fig. 1 for pa = 5. As in the case of normal incidence, 7 there is countable infinite number of FP modes having nearly equidistant real parts and finite imaginary parts. In addition there are two types of modes on the real ω-axis: WG and AWG modes, which are appearing for p = 0. The WG modes have an evanescent, i. e. exponentially decaying electric field into the vacuum, while the AWG modes are exponentially growing into the vacuum and are known in quantum-mechanics as anti-bound states. 23 Finally there is one leaky mode (LM) at the center of the spectrum which has zero real and negative imaginary part of ω. The functiont(ω) has two branch points at ω = ±pc connected by a cut, due to the square root in Eq. (6). We choose the cut going through ω = −i∞ and thus producing two vertical cut lines as shown in Fig. 1 other square root in the definition of q(ω) does not produce any cuts due to the above mentioned fact that t(ω) is an even function of q and thus independent of its sign. Integratingt(ω )/(ω−ω ) over a closed infinite-radius circular contour circumventing the cut, similar to that used in Ref. 8, we obtain the spectral representation in the frequency domaiñ (8) Here the first term represents a sum over residues at all poles oft(ω). The second term is the integral of the step ∆t(ω) in the transmission along the two parts of the cut shown in Fig. 1. Specifically, ∆t(ω) is defined as the difference between the values oft(ω) on the left and right sides of the cut for the given cut point ω.
Using the spectral representation Eq. (8) for real fre- quencies ω, we analyze contributions of the poles and the cut to the transmission. The transmission is usually studied for a fixed angle of incidence θ, motivated by experimental constraints. An example of the calculated transmission through a slab with s = 9 is shown for θ = π/4 in Fig. 2 (a). For a fixed θ, the in-plane wave vector p changes with frequency, so that the contributions of the poles (which are different for different p) are not constant across the spectrum. We therefore analyze the spectrum for a fixed p, as shown in Fig. 2 (b), in which the contributions of different pole types and the cut are shown individually, summing up to the analytic transmission Eq. (5). Note that the transmission t(ω) is defined over the angle range 0 < θ < π/2, corresponding to ω > pc. FP modes dominate for ω pc giving rise to the oscillations in the transmission, while the contribution of all other modes and the cut are significant only close to the threshold ω = pc, corresponding to grazing incidence θ ∼ π/2.
The cut contribution to the spectral representation Eq. (8) and to the transmission in Fig. 2 (b) produces a continuum of resonances. Such a continuum can be approximately treated in the RSE by replacing it with a series of poles, as done in Ref. 8. In the present case however, the cut can actually be removed by going into the wave-vector domain. Indeed, being treated as a function of the normal wave vector k, the AC of the transmissioñ T (k) has no cuts in the complex k-plane and its spectral representation obtained by using the Mittag-Leffler theorem has the following form: in which k n = ω 2 n /c 2 − p 2 , with n numbering the poles as in Eq. (8). On the real k-axis,T (k) coincides with the transmission T (k) given by Eq. (5) and is shown in Fig. 3 along with the contributions of the different types of modes. We see in particular that the WG modes, which are not emitting into an outgoing plane wave, and thus by reciprocity are expected not to be excitable by an incoming plane wave, have a finite contribution to the transmission, which is possible only due to their offresonant excitation. This contribution increases with decreasing the wave vector k, as the frequency of the incoming wave is getting closer to the resonant frequencies of the WG modes lying beyond the vacuum light cone.

III. RESONANT STATE EXPANSION FOR NON-NORMAL INCIDENCE
The RSE formulated in our previous works 6-8 is based on three key elements which are: (i) Dyson's equation for the Green's function (GF), (ii) spectral representation of the GF, and (iii) completeness of RSs used for expansions of perturbed states and the GF itself. We have previously applied the RSE to infinite 1D and 2D systems at normal incidence. The non-normal incidence, characterized by p = 0, is treated here. The previously used spectral representation of the GF in the frequency domain contains a cut for p = 0, which however can be removed by mapping the problem onto the complex normal wave-vector space k, as demonstrated in Sec. II. We therefore reformulate the RSE in the complex k-plane, for which the spectral representation of the GF of an infinite planar system with an in-plane momentum p = 0 for TE polarization 24 has the form where E n (z) is the electric field of a RS, defined as an eigensolution of Eq. (4) with an arbitrary profile of ε(z) within a finite interval |z| < a, satisfying the outgoing wave boundary conditions and orthonormality conditions 25 The GF satisfies the equation and thus has the asymptotics G k (z, z ) ∝ k −2 for k → ∞.
Applying it to Eq. (10), and using the fact that the GF has no pole at k = 0 corresponding to ω = ±pc, as exemplified in Fig. 1, we obtain the following sum rule: For p = 0 the right-hand side of the above sum rule is replaced by i, due to the k = 0 pole of the GF. 6 Using Eq. (14), one can write Eq. (10) as where F (k) is an arbitrary function which will be appropriately chosen later, in order to linearize a resulting matrix eigenvalue problem of the RSE. We now consider an arbitrary perturbation ∆ε(z) of the dielectric constant inside the layer |z| < a. The new, perturbed GF G k (z, z ) is related to the unperturbed one via the Dyson equation Substituting Eq. (15) into Eq. (16) and using a similar spectral representation for the perturbed GF in terms of the perturbed modes E ν (z), we equate, following Ref. 6, residua at the perturbed poles k = κ ν in Eq. (16). This results in the following relationship between unperturbed and perturbed modes Note that the perturbed modes E ν (z) satisfy Eq. (4) with ε(z) replaced by ε(z) + ∆ε(z) and the BCs Eq. (11) with k n replaced by κ ν . In the interior region |z| < a which contains the perturbation, the perturbed RSs can be expanded into the unperturbed ones, exploiting the completeness of the latter: Substituting this expansion into Eq. (17) and equating coefficients at the same basis functions E n (z) results in the matrix equation where is the matrix of the perturbation in the basis of unperturbed RSs. Equation (19) is a matrix eigenvalue problem which can be solved numerically in order to find the wave vectors κ ν and the corresponding eigenfrequencies of the perturbed RSs, as well as their expansion coefficients b nν in terms of the unperturbed ones. However, this problem is generally nonlinear in κ ν , as can be seen by choosing F (k) = 0. Nonlinear eigenvalue problems are known to lead to numerical instabilities and can produce spurious solutions. In order to avoid these issues, we choose explicitly depending on the in-plane wave vector p, which linearizes the eigenvalue problem. Indeed, with the substitution c nν = b nν k n /κ ν , the eigenvalue problem is given by which is linear and can be solved by inverting the matrix on the right-hand side of Eq. (22) and diagonalizing the resulting non-symmetric matrix on the left-hand side, in order to obtain its eigenvalues 1/κ ν . Alternatively, one can solve Eq. (22) by employing a variety of software libraries available for generalized linear matrix eigenvalue problems. Note that the matrix equation of the RSE for normal incidence previously derived in Ref. 6 is restored by choosing p = 0 in Eq. (22).

IV. RESULTS
To apply the method developed in Sec. III we first construct a convenient basis of unperturbed states. We use the RSs of a homogeneous dielectric slab discussed in Sec. II. We calculate the wave functions E n (z) of the RSs and investigate the dependence of their eigenvalues k n on the in-plane wave vector p. Then, using the RSE, in particular Eqs. (20) and (22), we calculate the perturbed eigenvalues κ ν for the simplest perturbation, which is constant across the slab, and compare the RSE results with the available exact solution. Finally, we use the RSE to treat a structured perturbation simulating a Bragg-mirror microcavity (MC). We specifically discuss the lowest-energy cavity mode (CM) and compare results with the transfer matrix calculation of the MC transmission and with an available analytic approximation for the CM linewidth.

A. Unperturbed resonant states
The solutions of Eq. (4) which satisfy the outgoingwave boundary conditions Eq. (11) in TE polarization take the form where the eigenvalues k n satisfy the secular equation with q n = s k 2 n + ( s − 1)p 2 . We use here an integer index n which takes even (odd) values for symmetric (antisymmetric) RSs, respectively. The normalization constants A n and B n are found from the continuity of E n across the boundaries and the normalization condition Eq. (12). They take the form where ω 2 n /c 2 = k 2 n + p 2 . The frequencies ω n of the RSs of a dielectric slab for pa = 5 and s = 9 were shown in Fig. 1. The normal wave vectors k n of the RSs for a slab with s = 3 versus p are given in Fig. 4. All states in the range |Re k n a| < 5 and |Im k n a| < 5 for |pa| < 5 are shown in Fig. 4 (a) and separated into mode types in Fig. 4 (b) and (c). For WG and AWG modes Re k n = 0, therefore Fig. 4 (c) shows only their imaginary part, which is positive for WG modes, corresponding to evanescent waves, and negative for AWG modes, corresponding to exponentially growing waves outside the slab. The WG and AWG modes continuously transform into each other and produce branches similar to those seen also for FP modes. These branches cross each other at certain points [shown in Figs. 4 (b) and (c) by magenta dots] where two FP modes are transformed into two AWG modes. The AWG mode branch which starts at p = 0 has no connection to any WG or FP branches; a mode on this branch was identified in Fig. 1 as the leaky mode.
The RSs of the homogeneous slab shown, similar to those shown in Fig. 4, are used as a basis for the RSE in the two examples given below. In general, for any local perturbation ∆ε(z) which does not change the translational symmetry of the slab, i. e. does not depend on x of y, the in-plane momentum p remains a good quantum number. In other words, ∆ε(z) does not mix states with different p, so that in any such problem solved by the RSE, one can use the basis of RSs with a given fixed value of p.

B. Full-width perturbation
To illustrate the accuracy and convergence of the RSE, we consider a homogeneous full-width perturbation of the slab, which is given by and for which the exact solution can be obtained by solving the transcendental Eq. (24) with s replaced by s + ∆ . We denote these exact perturbed wave numbers as κ (exact) ν and compare them with the perturbed values κ ν obtained by using the RSE for different basis sizes N . We choose as basis of given size all poles with |k n | < k max (N ), using a suitably chosen wave-number cutoff k max (N ).
In Fig. 5 we compare the RSE wave numbers with the exact wave numbers for our system in the case of pa = 5. We can see in Fig. 5(a) that the RSE is reproducing the exact solution to a high accuracy, which is quantified by the relative error |κ ν /κ (exact) ν − 1| shown in Fig. 5(b) for the FP modes with Re κ ν > 0 and in Fig. 5 (c) for the WG and AWG modes. We see that the relative error scales as N −3 , which was also observed in planar systems at normal incidence 6,7 and in homogeneous micro-cylinders 8 and microspheres. 6 We find in the simulation used to generate Fig. 5 for a basis of N = 2000 that the RSE reproduces about 300 modes with a relative error below 10 −8 . This error can be further improved by 1-2 orders of magnitude using the extrapolation method described in Ref. 7.

C. Microcavity perturbation
To evaluate the RSE for inclined geometry in the presence of sharp resonances in the optical spectrum, we use a Bragg-mirror MC, which consists of a FP cavity of thickness L C and dielectric constant C = 9 surrounded by distributed Bragg reflectors (DBRs). The DBRs consist of P = 5 pairs of dielectric layers with alternating high H = 9 and low L = 2.25 susceptibility, as illustrated by the inset in Fig. 6. The alternating layers have a quarter-wavelength optical thickness and the cavity has a half-wavelength optical thickness. The nominal wavelength which determines the layer thickness is that of the lowest-frequency CM at normal incidence. As unperturbed system we used a dielectric slab with s = 9 as in Sec. IV B.
The unperturbed modes of the slab and the perturbed modes of the MC are shown in Fig. 6 (a) for pa = 5. One can see how the nearly equidistant FP modes of the unperturbed system are redistributed in the MC, transforming into a sharp CM in the middle of a wide stop-band and modes outside of the stop-band. The link between the peaks in the transmission in Fig. 6 (b) and the poles in Fig. 6 (a) is also exemplified by the real part of the poles giving the position of the peaks in transmission and the imaginary part giving their linewidth. This is discussed in Ref. 7 in terms of the GF which is related to transmission via T (k) = 2ikG k (a, −a).
The transmission T (k) for a layered planar structure can be calculated using the transfer matrix method leading to the following explicit result: in which ξ + M is found from the recursive formula and the normal component of the wave vector in the j-th layer Here j and a j are, respectively, the dielectric constant of the j-th layer and its width, so that layers j = 0 and j = M correspond to the vacuum before and after the MC, respectively, so that q 0 = q M = k, and q j 0 for real j . M gives the total number of interfaces in the structure, in the present case M = 2(2P + 1).
In Fig. 7 we show the evolution of the perturbed poles with p. We see that one of the modes is separated in the middle of a gap and has an imaginary part well below the others. This mode is know as the CM. The perturbed Green's function G k (z, z ) which has a spectral representation equivalent to Eq. (10) and the corresponding transmission T (k) are dominated by the single term from the CM in this frequency region, therefore a sharp isolated peak is seen in the center of the stop-band in Fig. 6 (b). Interestingly, the modes in Fig. 7 show an almost circular behavior, indicating that the frequency of each mode ω ν = c κ 2 ν + p 2 is approximately constant versus angle θ.
Indeed, we can see in Fig. 8 (a) that the CM frequency ω C has a weak dependence on θ, while the corresponding wave vector κ C changes more strongly. In parallel, the linewidth given in Fig. 8 (b) shows a similar behavior both in the ω-and k-representations, though at θ → π/2 the imaginary part of ω C is one order of magnitude smaller than that of κ C . Figures 8(a) and (b) demonstrate a good agreement between κ C obtained using the RSE and κ We can also compare the results in Fig. 8(b) with an analytic approximation for the CM linewidth which we have derived by generalizing the approximation for normal incidence of light available in the literature. 7,26,27 Here n j is the refractive index of layer j, η j = n j cos(θ j ), and θ j is the angle to the normal in layer j, given by sin(θ j )n j = sin(θ). The layers j used are: the external region (ext) which is vacuum in our case, the high-index (H) layer, the low-index (L) layer of the Bragg mirror, and the cavity layer (C). The cavity wavelength is given by λ C = 2L C cos(θ C ). Equation (32) is exact in the limit P → ∞, for a structure with Braggmirror layer widths strictly equal to a quarter-wavelength (c) Relative error of κC determined by RSE for different basis sizes N as given. All data are shown as a function of the angle of incidence θ, and all symbols are connected by lines as a guide to the eye. and the cavity layer width to a half-wavelength optical thickness. This condition depends on the incident angle, and in our fixed structure is fulfilled for normal incidence only. Nevertheless, Eq. (32) reproduces the exact result reasonably well over the whole angle range, as shown in Fig. 8(b).

V. CONCLUSION
We have generalized the resonant state expansion to planar optical systems with inclined geometry. The method is based on the spectral representation of the Green's function of Maxwell's equation and expansion of the optical modes of a perturbed system into a com-plete set of resonant states of a simple dielectric slab. In inclined geometry, the spectrum of a planar system contains a continuum of resonances originating from a cut of the Green's function, which we have eliminated by mapping the frequency into the normal wave vector. The optical modes and spectra of a perturbed planar system are then calculated by solving a linear matrix eigenvalue problem containing matrix elements of the perturbation in the basis of discrete resonant states only. We have verified the method on full-width homogeneous and Braggmirror microcavity perturbations and compared results with obtained analytic solutions, demonstrating fast convergence of the method towards the exact result. We have recently demonstrated the application of RSE to two-dimensional open optical systems at normal incidence. We expect that we can extend this treatment to inclined geometry using a similar approach, which would provide an efficient algorithm to calculate the optical modes in fibers and waveguides, including photonic crystal fibers having a complex structure.