Analysis of guided and leaky modes of planar optical waveguides using transparent boundary conditions

Multilayer planar waveguides are of great interest in optics since they are basic parts of many photonic devices, such as semiconductor lasers,^1,2 sensors,³ Bragg reflectors, optical switches and modulators, spectral filters, waveguide polarizers,⁴ directional couplers,⁵ photonic crystal waveguides, plasmonic sensors, and metasurfaces.⁶ However, even for the simplest three-layer waveguide, the eigenmode equation is already transcendental, and numerical methods are required to find the waveguide modes.^7,8 A heuristic, graphical approach to the problem is often introduced in scientific textbooks for the purpose of demonstration.^7,9 For a planar waveguide of more than three layers of lossless or lossy materials, one can use the transfer matrix method that transforms the original mode solving problem into a problem of finding zeros of a complex-valued characteristic equation (the dispersion relation).^10,11 However, the matrix becomes numerically unstable when the number of layers increases and the floating error propagates across the layers and accumulates.¹² Worst scenarios happen if the transfer matrix method is applied to graded-index waveguides treated as a stuck of a large number of ultrathin layers of homogeneous materials. A more stable and versatile method to find the modes of a multilayer waveguide is a finite element method (FEM).

When using FEM, unbounded waveguides with semi-infinite superstrate and substrate are usually truncated by perfectly matched layers (PMLs), which are artificial anisotropic absorbing media,^13,14 or equivalently modeled by a complex coordinate stretching.¹⁵ The choice of the PML parameters, including the number of the elementary sublayers, their thickness, and the real and imaginary parts of the refractive index, affects the absorption efficiency of the PML. We emphasize that a PML designed to absorb the guided modes may cause significant errors for the leaky mode calculations.¹⁶ If, on the other hand, the PML is made to guarantee an effective absorption of the leaky modes, it fails to properly absorb the guided modes (see  Appendix B). Hence, a single choice of a set of the PML parameters is, in general, not suitable for both types of modes.

In this work, we generalize a one-dimensional FEM using transparent boundary conditions (TBCs)^12,17,18 to solve for the modes of an arbitrary unbounded multilayer waveguide. The semi-infinite regions of the waveguide are modeled rigorously by terminating them with TBCs at the two cladding interfaces, which avoids any numerical discretization and PML-related problems for the two unbounded cladding layers.¹² The TBC involves the transverse wave number that depends on the propagation constant, i.e., the eigenvalue we are searching for in the original problem.^12,17,18 Using TBCs, the equation may be formulated as an eigenvalue problem of the propagation constant^17,18 or the transverse wave number.¹² For the latter case, the original problem is transformed into a polynomial nonlinear eigenvalue problem written in terms of the two transverse wave numbers in the two semi-infinite regions. We develop a Python package to solve this nonlinear eigenvalue problem. The finite element part is implemented via the open-source library NGSolve.¹⁹ The nonlinear eigenvalue problem is solved by the Polynomial Eigenvalue Problem (PEP) solver of the Scalable Library for Eigenvalue Problem Computations (SLEPc).^20,21 In addition to numerical stability against possible PML-induced problems, the TBC approach we propose here permits one to perform a comprehensive and systematic analysis of the solved modes in terms of their power flow characteristics.

The paper is organized as follows: In Sec. II, we write the general equations for the electric and magnetic fields of a waveguide mode. Then, in Sec. III, we generalize a weak formulation expression of the waveguide mode problem that can be solved with FEM. After that, the approach is applied to a general planar multilayered waveguide, as described in Sec. IV. In Sec. V, we show how to use TBCs to simplify the numerical calculation process. Importantly, in Sec. VI, we systematically group modes into categories differing by their power flow characteristics. In Sec. VII, we analyze some typical examples of unbounded planar waveguides and discuss their properties. Finally, in Sec. VIII, our results are summarized and discussed.

Equations (26)–(30) make it possible to find β². To obtain β, we may choose a branch with the positive imaginary parts of β.^28,29 In the absence of optical gain media in the structure, our choice for β results in a picture that the overall energy carried by a mode must propagate in the positive y-direction [see the exponential factor $e−2κ′′k0y$ in Eqs. (31)–(34)]. If the gain is present, the overall energy propagation direction becomes dependent on whether or not the gain surpasses the loss. In the latter case, the mode still propagates in the positive y-direction. In the former case, however, the mode propagates in the negative y-direction, and its overall energy grows upon propagation. The real part of β, however, can be negative, showing a negative effective index of the mode. This can happen, for example, if the mode involves plasmons.^28–30 The modes for which β has equal values but different signs are the same, but counter-propagating modes, as dictated by the mirror symmetry of the structure with respect to the xz-plane. By examining the power flow in each layer, it becomes clear that not all of the obtained solutions are physically meaningful.

In a layer of a typical dielectric with ϵ′ ≫ϵ″≥ 0 [see the blue dot in Figs. 2(a) and 2(e)], three propagation regimes of the modes depending on the values of κ′ and κ″ can be recognized. In Fig. 2(b), the thick blue line drawn in the (κ′, κ″)-space corresponds to $S̄y(TM)=0$ and reflects the regime of a frozen power flow. This line divides the (κ′, κ″)-space into two parts. On the right-hand side of the line, the parameters make the power flow in the positive y-direction (red-colored regions U and W, in which $S̄y(TM)>0$⁠). On the left-hand side of the blue line, we have $S̄y(TM)<0$ (green-colored regions V and X). At κ″ > 0, the power in the layer (as well as the total power in the mode) is lower at larger values of y, independently of the light propagation direction. This is shown in the figure schematically by a reduced number of arrows. In contrast, when κ″ < 0, higher powers are observed at larger values of y. Note that in regions V and W, light power in the layer increases upon propagation, while in regions U and X it decreases. The longitudinal power flow is typically co-directed with the phase velocity. However, it becomes opposite to the phase velocity of the mode in the layer, for example, if $S̄y(TM)$ is positive and κ′ is negative [see the wedge in the red region U at negative κ′ in Fig. 2(b)]. The waveguide in this case behaves similarly to a negative-index material. In the green region at positive κ″, the Poynting vector in this specific layer points in the negative y-direction together with the phase velocity of the mode, even though the overall energy of the mode (integrated overall layers) propagates in the positive y-direction (if there is no optical gain in the structure, discussed at the beginning of this section). Typically, waveguides such as this contain metal layers.^28–30 

In the case of a metallic layer with −ϵ′≫ϵ″≥ 0 [indicated by the orange dot in Figs. 2(a) and 2(e)], the longitudinal power flow is usually opposite to the mode’s phase velocity [see the green and red regions in Fig. 2(f)]. The drawing in Fig. 2(f) is similar to that in Fig. 2(b), except for a different slope of the line $S̄y(TM)=0$ (the orange line) and interchanged green- and red-colored regions.

In practice, the modes to be used are often required to have a quality factor higher than a certain critical value, i.e., $κ′/κ′′>Qc$⁠, in order for the mode to propagate over at least Q_c wavelengths before its amplitude vanishes.²⁸ In Figs. 2(b) and 2(f), the gray dashed lines correspond to $κ′/κ′′=10$⁠. The points for modes with higher Q-factors are located closer to the κ′-axis. For a layer of an ideal lossless dielectric or metal with negligible ϵ″, the locations of the blue and orange dots in Figs. 2(a) and 2(e) approach the ϵ′-axis, and the blue and orange lines in Figs. 2(b) and 2(f) become vertical.

Using Eq. (36), we can also characterize the time-averaged transverse power flow in the substrate and cover layers. Different regimes of the power flow are identified, respectively, in Figs. 2(c) and 2(d) for a dielectric layer and Figs. 2(g) and 2(h) for a metal layer. Again, the blue and orange lines separate the regions of positive and negative values of $S̄x(TM)$ (red and green regions, respectively). In regions A and D, the power in the substrate or cover layer flows toward the waveguide core, while regions B and C correspond to the power flowing away from the core. In addition, the power density decreases with the distance from the core when $τs′′<0$ or $τc′′>0$⁠, and it increases otherwise.

Since the value of τ_s or τ_c may appear in any of regions A–D, the number of all possible combinations of the power flow types is 16 (see the table in Fig. 3 describing the case of κ″ > 0). However, the regimes corresponding to the gray regions in both Figs. 2 and 3 require the energy to be fed into the waveguide core from an unrealistic source located infinitely far away from the waveguide.³¹ We can, therefore, exclude these modes from practical consideration, reducing the number of combinations to 9. Note that the Poynting vector corresponding to region A_s/c also points toward the core, but this can happen, for instance, in a lossy guided mode when the core is absorbing and the cladding is not, in which case the energy flows from the evanescent part of the mode in the cladding toward the core to compensate for the energy loss. In the limit κ″ → 0, the guided mode is lossless and the transverse fields are purely evanescent, where the power flow direction is parallel to the propagation direction (see the limit cases of A_s–A_c and B_s–B_c enclosed by the dashed rectangles in Fig. 3). In practice, for dielectric slab waveguides, modes with their energy flowing from the more to the less optically dense cladding should also be excluded.²⁸ For example, if $ϵs′>ϵc′>0$ and $ϵs′′=ϵc′′=0$⁠, the mode cannot simultaneously have a power flow of type A_s–C_c (see the drawing on the purple background in Fig. 3). A quality factor threshold in the propagation constant, for example $κ′/κ′′>1$⁠, further restricts the number of modes to be retained.³⁰ Due to the considered restrictions, only in a certain region of the complex λ-plane [see Eq. (26)] we can find physical modes, which can help us in searching for the eigenvalues with the PEP solver.

The media in the waveguide may be either lossless or lossy. If all the media are lossless, we have κ″ = 0 for guided modes. Otherwise, we have κ″ > 0 even for a guided mode that does not radiate to the surroundings. In the case of κ″ = 0, we obtain either a pure real or a pure imaginary transverse wave number, as required by Eq. (38). The former case implies that the wave propagates in the cladding either toward or away from the core, which contributes to a continuum of the radiation modes.⁷ The latter case results in an evanescent or a growing transverse field in the cladding, meaning that the mode is either bounded or leaky.^7,32 For a bounded mode, the effective mode index must exceed the refractive index of the cladding, i.e., $(κ′)2>ϵi′$⁠, such that the transverse wave number is purely imaginary, as required by Eq. (37). A leaky mode must radiate and lose its energy into the cladding layers, and therefore, we have κ″ > 0, which contradicts the assumption that κ″ = 0. In other words, the bounded and guided mode solution requires $τi′=0$ and $τs′′<0$ or $τc′′>0$ so that the transverse fields in the substrate and cover layer are evanescent.

When a lossy material is present in the waveguide core, even a guided mode becomes lossy. The power of such a mode must decrease as the mode propagates, and consequently, we have κ″ > 0. In addition, the transverse intensity of the mode is evanescent, and therefore, it decreases with the distance from the core, i.e., $τs′′<0$ or $τc′′>0$⁠. When κ′ > 0, Eq. (38) written for the cladding layers leads to the inequality τ′τ″ < 0, which in turn implies that the total time-averaged Poynting vector is directed toward the core (see the power flow types U and A_s/c in Fig. 2 and A_s–A_c in Fig. 3). The other case, κ′ < 0, results in an opposite inequality, τ′τ″ > 0, and consequently, the Poynting vector is directed away from the core (see the power flow types V and B_s/c in Fig. 2 and B_s–B_c in Fig. 3). For a power flow corresponding to the case A_s–A_c in Fig. 3 when κ′ > 0, the intensity of the evanescent field in the lossless cladding decreases upon the mode propagation. This is explained by the fact that the energy in the cladding is lost as a result of flowing into the core, in which it is absorbed.²⁸ For the power flow corresponding to the case B_s–B_c in Fig. 3 when κ′ < 0, the intensity of the evanescent field in the cladding grows as the field propagates downward. This growing field seems counterintuitive since we have assumed that both the claddings are made of lossless dielectric materials and no gain medium is present. In fact, the energy in the cladding is increased along the longitudinal propagation direction due to the power flow from the core to the cladding. However, the total power carried by the mode still decreases in the upward direction (recall that κ″ > 0), implying that the energy flow in the core must be opposite to that in the cladding. It is also interesting that, in the limit of κ″ = 0 (see the insets outlined with the dashed rectangles in Fig. 3), the power flow is parallel or anti-parallel to the mode propagation direction and the transverse field becomes purely evanescent.

Leaky modes of waveguides can be divided into two types. In a so-called substrate leaky mode, energy “leaks” out of the waveguide core into the substrate. The value of τ_s for such a mode appears in region C_s in Fig. 2(c). In the substrate, the wave has a transversely outgoing component, with the intensity increasing with the distance from the core. On the other hand, the wave in the cover region is still bounded to the waveguide core, which is transversely incoming and decaying in its amplitude with the distance from the core [see region A_c in Fig. 2(d)]. The energy of this type of mode is partly absorbed in the core and partly radiating into the substrate. If the waveguide is lossless, the energy coming from the cover region will be transmitted completely to the substrate through the waveguide core. Such a substrate leaky mode is demonstrated schematically in the drawing of C_s–A_c on the light yellow rectangle in Fig. 3.

The modes leaking through both the substrate and the cover layers are called air leaky modes. The value of τ_i for such a mode must be in region C also for the cover layer (see the power flow C_s–C_c on the dark yellow rectangle in Fig. 3). As a consequence of energy leakage, the power of the mode inside the core decreases as the mode propagates, even if the waveguide is made of lossless dielectric materials. On the other hand, the intensity of the field radiated into the cladding increases with the distance from the core.

If one of the cladding layers is made of a lossless dielectric and the other of a lossless metal, then the directions of the longitudinal power flow in the cladding layers are opposite to each other [see Figs. 2(b) and 2(f) in the limit of an ideal dielectric and an ideal metal, respectively]. Repeating a similar analysis for the case of lossless dielectric cladding layers, we conclude that the power flow types corresponding to B_s–A_c, A_s–B_c, B_s–C_c, and C_s–B_c in Fig. 3 can exist.

In this example, we consider a planar waveguide with a core consisting of four layers of lossless media. Such a waveguide has also been studied previously.^10,33,34 The refractive indices of the six media of the waveguide are n_c = 1, n₁ = 1.66, n₂ = 1.53, n₃ = 1.60, n₄ = 1.66, and n_s = 1.5. The thickness of each core layer is the same: d₁ = d₂ = d₃ = d₄ = 500 nm. The waveguide modes are studied at a wavelength of 632.8 nm. In the FEM calculation, we use 31 second-order Lagrange finite elements per wavelength, resulting in 201 degrees of freedom (DoFs) in the complex-valued finite element space.¹⁹ In addition, the SLEPC’s PEP module is configured with the tolerance and maximum number of iterations being 1 × 10⁻⁸ and 50, respectively, and to search for eigenvalues that outnumber the DoFs by four times. On a Linux laptop equipped with an 8-core 11th Gen Intel(R) i5-1135G7@2.4 GHz processor, the simulation for all the modes takes about 10 s. We noticed that when a large amount of the eigenvalues are requested, the actual iteration is small. In the case of using the default TOAR solver in the PEP module, it takes three iterations. In Fig. 4, the modes obtained with our FEM simulations are represented by the circles. For comparison, the red crosses correspond to the results extracted from the literature for the modes calculated using the transfer matrix method.^10,33,34 The gray circles stand for the solutions that are not physical because the modes have either of the following three properties: (1) τ_s corresponds to the regime D_s, (2) τ_c corresponds to the regime D_c, and (3) the power flow is of type A_s–C_c (when n_s > n_c). The physical modes with $κ′/κ′′>1$ are marked by the colored circles, with each mode having its own color. In Figs. 4(a)–4(d), we show all the revealed modes in their κ-plane, λ-plane, τ_s-plane, and τ_c-plane, respectively.

Four guided modes, for which κ′ > n_s, are found by both methods [see the four blue circles that overlap with the large red crosses in the zoomed-in part of Fig. 4(a)]. The corresponding transverse fields of the four guided modes are purely evanescent outside the core, as their τ_s values lie on the negative $τs′′$-axis in Fig. 4(c) and the τ_c values on the positive $τc′′$-axis in Fig. 4(d). As follows from Eq. (26), the four guided modes must have positive values of λ″, which is the case in Fig. 4(b). In the range of n_c < κ′ < n_s, five modes are found by both the methods, as indicated by the overlapping colored circles and red crosses in the region between the two vertical dashed lines in Fig. 4(a). It is evident that they are substrate leaky modes that reside in the C_s region of the τ_s-plane and the A_c region of the τ_c-plane. We discover that the transfer matrix method failed to find the other four physical modes with n_c < κ′ < n_s, which are marked by the four dark green circles without red crosses in Figs. 4(a)–4(d). They behave as air leaky modes because τ_s lies in the region C_s and τ_c is in the region C_c. These modes were overlooked in Refs. 10, 33, and 34, since only one of two possible transverse wave numbers in the substrate/cover layer has been chosen [for example, when using Eq. (25) in Ref. 10]. The missed modes have an opposite sign in the transverse wavenumber τ_c with respect to the other modes in the same range. Finally, in the range of κ′ < n_c, three air leaky modes are found by both the methods, but three substrate modes [see the light green markers that are not overlapping with any red crosses in Fig. 4(a)] are overlooked by the transfer matrix method for the same reason of selecting only one sign in the τ_c values.

We next consider another four-layer structure that is the same as the previous one, except for the two core layers that are made weakly absorbing by changing their refractive indices to n₁ = 1.66(1 + 10⁻⁴i) and n₂ = 1.53(1 + 10⁻⁴i).^10,34 The simulation is performed using the same discretization as in the corresponding lossless case, and the computational time is approximately the same. The inserted low loss causes no significant disturbance to the modes in the waveguide [see Figs. 5(a) and 5(b) that greatly resemble those of Figs. 4(a) and 4(b)]. Again, the unphysical modes and the modes of a low quality factor with $κ′/κ′′<1$ are marked by gray circles. A slight difference is revealed in the inset of Figs. 5, 5(c), and 5(d), where six extra modes (colored circles not overlapping with any red crosses) are found in the range of κ′ > n_s. Three of them behave as substrate leaky modes, with their τ_s values located in the region C_s and τ_c values in the region A_c. The other three modes have the τ_s values in the region C_s and τ_c values in C_c, and therefore, they are air leaky modes. The former modes appear because the energy transmitted from the cover to the core is partially absorbed and partially leaked into the substrate. In the latter, air leaky modes, the energy is absorbed in the core and radiated into both the cladding layers. This implies that the condition κ′ > n_s does not guarantee that the mode will be guided by the structure, especially for waveguides with multiple layers of close refractive indices.

Our FEM simulation results are plotted in Fig. 6. They agree with the analytical solutions presented in Ref. 30. At frequency ω = 0.74ω_p, the κ value of the TM₀ mode (the blue star or blue curve in Fig. 6) is located in the V region for the cladding layers [see Figs. 2(b) and 6(d)] but in the U region for the gold core [see Figs. 2(f) and 6(d)]. The anti-parallel power flows with $S̄y(TM)<0$ in the claddings and $S̄y(TM)>0$ in the core are also reflected in the normalized $S̄y(TM)$ profile [see the green curve in Fig. 6(b)]. The τ_c value of the TM₀ mode in the region B_c of Fig. 6(e) indicates a B_s–B_c power flow type, where the energy is radiated from the core to the claddings [see the purple dotted curve in Fig. 6(b)]. As discussed earlier, one may expect that the leaked energy in the claddings should give rise to a growing power flow, which contradicts the assumption of no gain in the structure. Indeed, the power flow in the claddings increases in the direction opposite to the mode propagation direction, which is the same as that in the gold core. Hence, the power flow in the claddings actually decreases as the mode propagates [see the power flow type B_s–B_c in Fig. 3 and assume the upward directions of the mode propagation and power flow in the gold core], which removes the contradiction. On the other hand, the TM₁ mode has its κ value in the X region of the cladding and W region of the core [see the orange star in Fig. 6(d)], and its τ_c value is in the A_c region [see the orange star in Fig. 6(e)]. Therefore, the TM₁ mode has a power flow type A_s–A_c, similar to that of the guided modes of the four-layer lossy waveguide, except for the negative direction of the power flow inside the gold core.

In an MDM waveguide, we consider two semi-infinite gold claddings separated by an air gap of 80 nm. We use 22 second-order Lagrange finite elements for the core layer, and the spectrum of 65 wavelengths is obtained in about 3 s. Due to the negative ϵ′ of the gold claddings, the longitudinal power in the claddings points in the negative y-direction. This can be seen in Fig. 7(d), where the κ values of the two modes at normalized frequency ω/ω_p = 0.57 reside in the V region of the (κ′, κ″)-plane, and the normalized values of $S̄y(TM)$ in the claddings are negative [see the green curves in Figs. 7(b) and 7(c)]. In the transverse direction, the power flows from the core to the claddings, which is seen from Figs. 7(b) and 7(c); the normalized values of $S̄x(TM)$ shown by the dashed purple curve are positive in the cover medium and negative in the substrate. The power flow type is then B_s–B_c, but it is slightly different from that of the mode TM₀ of the DMD waveguide. Here, the mode propagation direction is parallel to the power flow in the dielectric core and anti-parallel to that in the gold claddings. In addition, the energy loss in the core is due to the radiation into the gold claddings, where light is absorbed.

In so-called metal-clad waveguides,³⁵ the core is separated from the substrate by a metal film, allowing the refractive index of the core to be independent of that of the substrate. Here, we investigate a planar metal-clad waveguide, which consists of a semi-infinite silver substrate, a 200-nm thick silica slab as a core with a refractive index of 1.4537, and an unbounded air cover. We use 25 second-order Lagrange finite elements in the core layer, and the spectrum of 60 wavelengths is calculated in 20 s. The refractive index of silver is taken from Johnson and Christy’s tabulated data.³⁶ At normalized frequency ω/ω_p = 0.57, where ω_p = 1.376 × 10¹⁶ (rad/s), the simulation results reveal four modes of relatively high quality factors, i.e., $κ′/κ′′>1$⁠, which are represented by colored stars in Figs. 8(a), 8(f), 8(g), and 8(h). It is evident that the modes TM₀ and TM₁ have nearly the same κ and τ_s values [where the blue star is hidden under the orange star in Figs. 8(f) and 8(g)], but opposite τ_c values [the blue star is in the region A_c, whereas the orange star is in the region C_c of Fig. 8(h)]. The TM₀ mode has a power flow type B_s–A_c, and the energy in the air cover is transmitted through the lossless silica slab to the silver substrate, where it is absorbed. On the other hand, the energy of the TM₁ mode leaks from the silica core into both the silver substrate and the air cover, manifesting a power flow type B_s–C_c. In the substrate, the leaked energy is absorbed. However, the energy leaked into the air cover is directed parallel to the mode propagation direction, resulting in an energy accumulation and, therefore, a growing transverse field. This is in contrast to the bounded transverse field in the substrate, where the leaked energy is directed anti-parallel to the mode propagation direction and the transverse field decreases away from the cover interface. The power flow types of the TM₂ and TM₃ modes are similar to those of TM₀ and TM₁, respectively. It should be noted that the τ_s value of the TM₂ mode is still located in the region B_s, where $S̄x(TM)<0$⁠. Therefore, the TM₂ mode does not have a power flow type A_s–A_c, even though it obeys the inequality $τs′<0$⁠.

We presented a numerical approach that allows one to investigate and comprehensively characterize both guided and leaky modes of multilayered planar optical waveguides. The approach is based on FEM using TBCs. We outlined a general weak formulation for the FEM-based analysis and derived a linear eigenvalue problem, Eq. (18), for the modes of planar waveguides composed of multiple layers of different materials. We modeled the two unbounded cladding regions analytically, which removed the necessity to treat them numerically using, e.g., PMLs. This was done by applying the TBCs to the cover-core and substrate-core interfaces. We have then transformed the weak formulation into a nonlinear eigenvalue problem, Eq. (30), that can be solved using the PEP module of the numerical library SLEPc. We have also categorized the modes in terms of their possible transverse and longitudinal power flow directions and densities, $S̄x$ and $S̄y$ (see Fig. 3). We have shown how to find physical and unphysical solutions and how to select the modes that can be relevant for practical applications, e.g., the ones with |κ′/κ″| > Q_c. Several typical power-flow types were considered in more detail using the examples of four-layer lossless and lossy waveguides, MDM and DMD plasmonic waveguides, and an asymmetric metal-clad waveguide, which are all common optical components in integrated optics. We have made our approach accessible as an open-source Python package at GitHub.³⁷ It provides researchers with a framework to systematically analyze and characterize eigenmodes of planar optical waveguides.

In contrast to FEM without TBCs, where the unbounded claddings are usually truncated with PMLs, the TBC approach is particularly suitable for the eigenmode analysis in planar optical waveguides. First, the TBC approach avoids the numerical discretization of the unbounded claddings, and therefore, it is numerically more stable. Moreover, usually, the absorption performance of a special PML design is determined by the transverse wavenumber, which is different for waveguide modes of different propagation constants. Therefore, it is challenging or even impossible to accurately calculate both guided and leaky modes using a single PML design. Another advantage of applying TBCs is that the nonlinear eigenvalue solutions can be intrinsically mapped to the parameter space of complex-valued τ_s and τ_c, using Eqs. (27) and (28), making it possible to uniquely determine the power flow properties of all the calculated modes. This is not the case, however, when using a PML-based approach, in which the transverse wavenumber τ is obtained from Eqs. (22) and (23) and the sign of the square root is undetermined.

It remains to figure out how to use TBCs for studying the scattered fields in the coupling of light from a source to a waveguide because the scattered field usually involves several modes with different transverse wavenumbers. In addition, it is yet to explore the application of FEM with TBCs to a three-dimensional periodic lattice of planar layers. It is critical to have an analytical form for the field profile beyond the interface where the TBC is to be applied. In addition, the challenge to be expected is that more terms will emerge in the nonlinear eigenvalue problem since the Bloch modes will exhibit all three vector components of the electric and magnetic fields.

We acknowledge the Research Council of Finland Flagship Program, Photonics Research and Innovation (PREIN; decision number 346529, grant number 320167). We thank the developers of the open-source finite element library NGSolve and the nonlinear eigenvalue computation library SLEPc. We also appreciate the computational resources provided by the Aalto Science-IT project.

The authors have no conflicts to disclose.

Xiaorun Zang: Software (equal); Writing – original draft (equal); Writing – review & editing (equal). Andriy Shevchenko: Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data and codes that support the findings of this study can be found in Github at https://github.com/shelvon/ngslab/.

We now calculate the mode effective indices κ of the same four-layer lossless waveguide as in Fig. 4 by means of FEM simulations using both PML-based and TBC-based methods and then compare them with the reference results obtained previously by using the T-matrix method.³³ In Fig. 9, we show the values of κ and their relative differences obtained with the PML scaling profile s′ + is″ = 1 + 10|Δx/L|² + 10i|Δx/L|². In addition, the PML layer is placed 1.5 wavelengths away from the cladding interface and has a thickness of 1 vacuum wavelength. It is discretized by 17 sublayers of equal thickness. This PML-based design results in approximately the same guided modes as the TBC approach. However, the PML-based simulation for leaky modes with the effective indices close to the refractive indices of the claddings degrades in accuracy. In an attempt to increase the accuracy for the leaky modes, we try another scaling profile s′ + is″ = 1 − 2|Δx/L|² + 10i|Δx/L|² while keeping all the other parameters unchanged. Although the accuracy slightly increases, the PML design results in relatively large errors in the guided mode calculations (Fig. 10).