## Abstract

We outline a full-vectorial three-dimensional multi-mode matching technique in a cylindrical coordinate system that addresses the mutual coupling among multiple modes co-propagating in a perturbed whispering gallery mode microcavity. In addition to its superior accuracy in respect to our previously implemented single-mode matching technique, this current technique is suitable for modelling waveguide-to-cavity coupling where the influence of multi-mode coupling is non-negligible. Using this methodology, a robust scheme for hybrid integration of a microcavity onto a silicon-on-insulator platform is proposed.

© 2014 Optical Society of America

## 1. Introduction

Whispering gallery mode (WGM) microcavities display the uniquely combined features of ultrahigh quality factors (Q’s) and small mode volumes. This has attracted research interest towards a wide range of applications, spanning single-molecule biosensing to on-chip optical filters [1–8]. Utility with regards to silicon photonics, however, is limited due to the relatively low quality factors of the microcavities available on silicon photonic platforms (e.g. silicon-on-insulator or SOI). Due to fabrication inaccuracies and the state of today’s technology, the highest quality factor for a microcavity fabricated on an SOI wafer is reported to be 5×10^{6} [9]. To implement a cavity with a quality factor above 10^{8}, such as a silica microtoroid reported in [1], a hybrid integration method is required [10]. In the search for an optimal hybrid integration scheme with minimal alignment requirements, numerical modelling then becomes essential. Several techniques (with varied pros and cons) have been developed in the past to model cavity-to-waveguide coupling [11–15] as well as to analyze general properties of microcavities [16–20]. We have recently demonstrated a mode matching method (MMM) that models the above problems [21]. In this article, the matching of single whispering gallery modes propagating along the azimuthal direction of a cavity perturbed by nanoparticles was performed. Such a technique yields highly accurate results when the dominant portion of light is confined within a single whispering gallery mode on the course of propagation. Here, we generalize the MMM analysis by incorporating multiple WGM’s co-propagating in the perturbed region of a cavity. The validity of the newly developed algorithm is illustrated by investigating a tapered waveguide coupled to a silica microtoroid, wherein the mutual coupling between WGM’s invalidates the single-mode matching method. This modelling technique allows for the examination of various schemes for hybrid integration of a whispering gallery microcavity onto an SOI platform and for the identification of a robust integration approach.

## 2. Theoretical formulations

An ideal whispering gallery microcavity (with refractive index independent of the azimuth *ϕ*) may support multiple transverse modes of the same azimuthal order *M* distinguished by resonance wavelengths. Consequently, continuous wave (CW) laser light delivered to the cavity excites a specific mode if its optical wavelength coincides with the resonance wavelength of that mode. The corresponding electric field distribution of an *l*^{th}-order transverse mode can be expressed as **E**(**r**) = *Aẽ̂ _{M,l}*(

*ρ*,

*z*)

*e*

^{jm̃lϕ}. The term

*m̃*=

_{l}*M*+

*jm̃*is a complex number whose real part is an integer

_{i,l}*M*representing the azimuthal order of the WGM and whose imaginary part

*m̃*characterizes the intrinsic quality factor of the mode according to

_{i,l}*Q*=

*M*/2

*m̃*. The unit vector

_{i,l}*ẽ̂*(

_{M,l}*ρ*,

*z*) is the normalized

*ϕ*-independent mode field distribution such that the squared amplitude |

*A*|

^{2}represents the total energy of the light stored in the cavity (in units of Joules). For convenience, we have chosen a cylindrical coordinate system (

*ρ*,

*ϕ*,

*z*) as defined in [21]. In a whispering gallery microcavity where loss exists and the field extends to infinity, there is a resultant departure from ”real” modes as defined in traditional lossless straight waveguides. For a high-Q low-loss cavity

*m̃*≪ 1 and, with a good approximation, the discrete set of normalized

_{i,l}*M*

^{th}-order azimuthal whispering gallery modes

*ẽ̂*(

_{M,l}*ρ*,

*z*)

*e*

^{jm̃lϕ}span an orthogonal basis in Hilbert space. In a cylindrical coordinate system [22],

*ε*(

*ρ*,

*z*,

*ϕ*) =

*ε*(

*ρ*,

*z*)) and Eq. (1) can then be simplified to

Here, *δ _{lk}* is the Kronecker delta and bra-ket notation is used to define the inner product in the functional space spanned by the WGM’s as 〈

*a*|

*b*〉 =

*π*∬

*εa*

^{*}

*·bρdρdz*where

*a*

^{*}is the complex conjugate of

*a*. For convenience, we also drop the subscript

*M*by limiting our discussion to the same azimuthal order

*M*. Also note that in the case of a large cavity, where the field is tightly focused into a small spot of mean radius

*R*and proper normalization is adopted, the orthonormality condition Eq. (2) can be approximated by $\iint {\widehat{\tilde{e}}}_{l}^{*}(\rho ,z)\cdot {\widehat{\tilde{e}}}_{k}(\rho ,z)d\rho dz={\delta}_{lk}$ [21]. The loss

*e*

^{−m̃i,lϕ}is extracted while rendering

*ẽ̂*(

_{l}*ρ*,

*z*) lossless. In this sense, the orthogonality condition for

*ẽ̂*(

_{l}*ρ*,

*z*) will continue to hold. The integral additionally extends to the full space while, in practice, one numerically calculates the modal field with a finite computation window and hence resolves the integral within that window. Such an approximation reduces accuracy due to the neglect of the field radiating to infinity; however, for cases in which the material and mode mismatch losses dominate over the radiation loss, this provides sufficient accuracy for high-Q cavities where the radiation field is small.

A non-ideal microcavity has a refractive index profile that depends on the azimuthal angle. To formulate the multi-mode matching technique, we assume light with wavelength *λ* close to a non-ideal cavity resonance wavelength propagates in the cavity. The non-ideal microcavity can be subdivided into many sections along the azimuthal direction, each of which can be approximately considered as part of an ideal cavity that shares the same refractive index profile as that section. At a cross section at an arbitrary angle *ϕ*_{0}, one may find a orthonormal set of *N*(*ϕ*_{0}) WGM’s {*ẽ̂*_{1}(*ρ*, *z*, *ϕ*_{0}), *ẽ̂*_{2}(*ρ*, *z*, *ϕ*_{0}),..., *ẽ̂*_{N(ϕ0)}(*ρ*, *z*, *ϕ*_{0})} corresponding to the approximated ideal microcavity. Thus, any field pattern at *ϕ*_{0} can be expressed as a linear superposition of these modes:

*δϕ*

_{0}, the field at

*ϕ*

_{0}+

*δϕ*

_{0}evolves to

*m′*is a

_{l}*ϕ*

_{0}-dependent complex number for the

*l*

^{th}mode at the operating wavelength

*λ*that can be estimated according to [21]. Meanwhile, the field

**E**(

*ρ*,

*z*,

*ϕ*

_{0}+

*δϕ*

_{0}) can be expanded onto the new set of orthonormal modes at

*ϕ*

_{0}+

*δϕ*

_{0}:

*ẽ̂*(

_{k}*ρ*,

*z*,

*ϕ*

_{0}+

*δϕ*

_{0})|, and utilizing the orthonormality condition (i.e. Eq. (2)), we obtain

*A⃗*(

*ϕ*) ≡ {

*A*

_{1}(

*ϕ*),

*A*

_{2}(

*ϕ*),...,

*A*

_{N(ϕ)}(

*ϕ*)}

*as an*

^{T}*N*(

*ϕ*)-dimensional vector whose

*l*

^{th}element represents the expansion coefficient of the

*l*

^{th}mode.

*C̃*(

*ϕ*

_{0}) is an

*N*(

*ϕ*

_{0}+

*δϕ*

_{0}) ×

*N*(

*ϕ*

_{0}) matrix whose (

*k*,

*l*) elements are

*e*

^{−mklδθ}= 〈

*ẽ̂*(

_{k}*ρ*,

*z*,

*ϕ*

_{0}+

*δϕ*

_{0}) |

*ẽ̂*(

_{l}*ρ*,

*z*,

*ϕ*

_{0})〉 similar to those defined in [21]. To obtain the resonance mode of a non-ideal cavity, one may divide it into

*P*slices at

*ϕ*= {

*ϕ*

_{0},

*ϕ*

_{1},...,

*ϕ*

_{P−1}} and propagate the field one round trip by progressively applying Eq. (7):

*A⃗*(

*ϕ*

_{0}+ 2

*π*) =

*e*

^{j2πm}

*A⃗*(

*ϕ*

_{0}) where multiple values of

*e*

^{j2πm}can be obtained by solving for the eigenvalues of the overall matrix in Eq. (8), each of which correspond to a transverse mode of the non-ideal or perturbed cavity. By analogy to the case of an ideal cavity, we may also label the

*m*of each transverse mode as

*m*. Note that resonance only occurs at a wavelength

_{l}*λ*such that the real part of

*m*is an integer. One may estimate the resonance wavelength

*λ*and the quality factor

_{p}*Q*of a specific non-ideal or perturbed whispering gallery mode according to ${\lambda}_{p}=\lambda \frac{\text{Re}\{{m}_{p}\}}{M}$ and ${Q}_{p}=\frac{M}{2\text{Im}\{{m}_{p}\}}$ following the same argument in [21]. The corresponding modal field distribution at

_{p}*ϕ*

_{0}can be obtained from the eigenvector

*A⃗*(

_{p}*ϕ*

_{0}) as

*ϕ*-dependent. To reiterate, one may obtain the field distribution at any azimuthal angle by progressively applying Eq. (6) starting from

*ϕ*

_{0}.

It is worth pointing out that, when a single whispering gallery mode remains, *A⃗* becomes a scalar and Eq. (8) reverts to its simplified form as seen in the single-mode propagation equation (i.e. Eq. (11) of [21]).

## 3. Application and discussion

The aforementioned algorithm was used to bring about a rigorous hybrid integration scheme for stationing an ultrahigh-Q microtoroid on an SOI platform (as illustrated in Fig. 1a). The validation of this approach could lead to a two orders of magnitude improvement in the quality factor for SOI integrated optical cavities, therefore opening the way for a novel series of products based on such arrangements. As shown in the top right inset of the plot, the SOI waveguide under investigation has a width of 500 nm and is formed via an upper Si layer as well as a lower SiO_{2} layer of 220 nm and 2 *μ*m respective thicknesses. The refractive indices of silica (1.44462) and silicon (3.48206) are taken from the literature [23,24]. In our model, we employ a microtoroid that has a major radius of 45 *μ*m and a minor radius of 5 *μ*m.

In this test case we study the coupling scheme between a straight waveguide and a toroid, as is shown in the lower inset of Fig. 1a. Light propagation along the azimuthal direction is computed at a wavelength of 1500 nm. In the first simulation, we placed the straight SOI waveguide in contact with the microtoroid’s equator and a 289^{th} fundamental quasi-TE toroid mode was launched at *ϕ* = −0.34 rad, wherein the waveguide was 3 wavelengths away from the cavity. The launched transverse electrical field distributions are displayed in the two left insets of Fig. 1b. Beyond that point, the coupling between the cavity and the waveguide falls to a negligible level (with relative Q degradation below 10^{−10}) and so it is neglected. As shown in the main plot of Fig. 1b, the amplitude of the launched mode diminishes while that of the second lowest-order mode reaches its maximum when the light propagates towards the center of the coupling region. It is also evident from the two right insets of Fig. 1b that the cavity mode launched as the input gradually evolves to a strongly hybridized cavity-taper mode when the light is approaching *ϕ* = 0 rad. The hybrid mode eventually relaxes to the unperturbed cavity mode at the exit of the cavity-taper coupling region and most of the light energy is coupled back to the originally launched mode in a similar manner as that of a directional coupler. By computing the total energy loss of the launched mode subtracted by the contribution from the intrinsic Q, we obtain a coupling Q-factor of 10^{5}. The coupling Q calculated here takes the loss from the launched mode to both the tapered waveguide and higher-order cavity modes into account. To precisely estimate the energy delivered to the waveguide, one may simply compute the overlap factor between the output field at *ϕ* = 0.34 rad and the mode of the tapered waveguide. To characterize the convergence rate of our algorithm, we plotted the reciprocal of the quality factor as a function of azimuthal angle steps and the number of modes included (given by the blue triangles in Figs. 1c and 1d). The expectation value of this quantity at an infinitesimal azimuthal step or infinite mode numbers (i.e. blue dashed line) was extracted through the Richardson extrapolation procedure [25] and the relative error (i.e. red cross markers) was estimated with the extrapolated value as a reference. The black line in Fig. 1c depicts a convergence rate of *O*(Δ*ϕ*^{1.0}), while that of Fig. 1d yields a convergence rate of about *O*(*N*^{0.8}).

The mathematical interpretation of the convergence rate is as follows: in Eq. (4), it is assumed that *m′ _{l}* remains constant between

*ϕ*

_{0}and

*ϕ*

_{0}+

*δϕ*. A more generalized form which takes the continuous variation of

*m′*into account should be

_{l}*m′*is a well behaved function between

_{l}*ϕ*

_{0}and

*ϕ*

_{0}+

*δϕ*, we can expand it as a power series around

*ϕ*

_{0}: ${m}_{l}^{\prime}(\varphi ,\lambda )={{\sum}_{n=0}^{\infty}\frac{1}{n!}\frac{{\partial}^{n}{m}_{l}^{\prime}(\varphi ,\lambda )}{\partial {\varphi}^{n}}|}_{\varphi ={\varphi}_{o}}{(\varphi -{\varphi}_{o})}^{n}$. By retaining the first term of the power series, one can obtain Eq. (4). Note that for a sufficiently small angle evolution

*δϕ*where the second power series term dominates the error,

*m′*(

_{l}*ϕ*) =

*m′*(

_{l}*ϕ*

_{0}) +

*O*(

*δϕ*). The 1/

*Q*term characterizes the overall loss contributions from both

*m′*, defined as the imaginary part of

_{i,l}*m′*, and the mode mismatch loss

_{l}*m*. The accuracy of the latter term is determined by the transverse grid spacing, the number of modes involved in the computation, and the relative insensitivity to

_{kl}*δϕ*for small azimuthal steps. Consequently, the adequately small azimuthal steps will procure a convergence rate for 1/

*Q*as a function of

*δϕ*determined by the convergence of

*m′*, which is of

_{i,l}*O*(

*δϕ*). At larger angles, this will be influenced by higher-order power series terms and the mode mismatch loss will deviate from the convergence rate of

*O*(

*δϕ*).

The convergence rate as a function of the number of modes that are included, on the other hand, is case dependent. One may expect that for the extreme condition of an ideal cavity where no mode coupling occurs, the convergence rate should be *O*(1). In the particular simulation of our structure, the convergence rate is *O*(*N*^{0.8}). The least square fit excludes the *N* = 1 point, as it has larger relative error. This is caused by the odd parity at *N* = 1 and sampling at the center of the structure, whereas the rest of the simulations have even parity and are sampled symmetrically about the center.

To reach the critical coupling condition, wherein the coupling Q-factor (*Q _{c}*) is equal to the intrinsic Q-factor of the cavity (assumed to be a practical value of 2×10

^{8}, as specified by the dot-dashed line of Figs. 2a and 2b), we gradually displaced the tapered waveguide horizontally away from the toroid. As displayed in Fig. 2a, a

*Q*on the order of 10

_{c}^{5}is computed when the waveguide touches the cavity surface and 10

^{10}is computed at a gap size of 2.5 wavelengths. For a waveguide sitting in the equatorial plane, a gap size of 0.75

*μ*m is also determined to be desirable in establishing critical coupling. The dependence of the Q-factor as well as the coupling parameter K on the gap is akin to the data in [26] for a microsphere to fiber-taper coupling system. Our K data points (i.e. red cross markers, for a 67

*μ*m-diameter microsphere and 1.35

*μ*m-diameter fiber at

*λ*=1550 nm) for the low-Q/accuracy-inhibiting regime in Fig. 2c are in good agreement with the experimental results of [26] (i.e. blue curve, for identical specifications), thus appropriately substantiating our methodology. It is evident from Fig. 2a that the coupling Q-factor is sensitive to the gap distance. More specifically, a gap distance fluctuation on the order of 300 nm may cause a drop of the coupling Q-factor from 10

^{8}to 10

^{7}. Therefore, highly precise alignment is required in order to integrate a toroid onto an SOI platform with this scheme.

To circumvent the alignment restriction, one can keep the waveguide in contact with the cavity and move it vertically by adjusting the thickness of the insulator layer as to tune the value of the coupling Q-factor. In such a configuration, one may pinpoint the location of the waveguide via an angle *θ* in respect to the equatorial plane, as illustrated in the inset of Fig. 2b. The Q-factor as a function of *θ* is shown in the main plot as solid circle markers, wherein a Q-factor of 3 × 10^{5} is calculated at *θ* = 0. It is also worth mentioning that the *θ* = 0 position (corresponding to the top surface of the SOI waveguide being aligned with the equator) is different from that in the previous case (i.e. the midplane of the silicon layer being aligned with the equator). At a larger *θ*, the local field intensity of the cavity mode is smaller. As a result, coupling between the cavity and the waveguide is weaker and a larger coupling Q-factor is observed. The right axis of Fig. 2b represents the magnitude of the electric field squared |*E*|^{2} on the surface of the cavity at different *θ*. It is important to recognize that the Q-factor function inversely resembles the |*E*|^{2} curve. For a straight SOI waveguide physically touching the toroid surface and that is placed below the equator, a waveguide location at *θ* = 65° is favourable in achieving critical coupling. In contrast to the previous case, a waveguide height fluctuation of 1 *μm* yields a drop of coupling Q from 10^{8} to 10^{7}. Control of such misalignment tolerance can be easily achieved with conventional technology.

## 4. Conclusion

A three-dimensional full-vector multi-mode matching method was formulated in cylindrical coordinates and tested. This technique solves generalized whispering gallery mode cavity problems, yet in this paper was specifically applied to a toroid-SOI coupling geometry where no precise control of the gap distance between the cavity and the waveguide was required. Furthermore, K parameters for many microsphere-fiber taper separations were assessed for accuracy restraining criteria. Once again, the multi-mode matching method proposed in this article is applicable to numerous classes of WGM related scenarios where the energy transfer between different modes is non-negligible. Two constraints remain in the current implementation: the assumption of neglecting the field radiating to infinity and the backscattering. These restrictions, however, can be lifted by, e.g., considering perfectly matched layer modes similar to those in [27] as well as adopting a bidirectional mode matching method. These additions, which have been used in conventional waveguide modelling, are currently under investigation in conjunction with the experimental verification of the integration scheme. Advancements of this kind will appear in subsequent publications.

## References and links

**1. **D. Armani, T. Kippenberg, S. Spillane, and K. Vahala, “Ultra-high-Q toroid microcavity on a chip,” Nature **421**, 925–928 (2003). [CrossRef] [PubMed]

**2. **F. Vollmer and S. Arnold, “Whispering-gallery-mode biosensing: label-free detection down to single molecules,” Nature Methods **5**, 591–596 (2008). [CrossRef] [PubMed]

**3. **V. R. Dantham, S. Holler, C. Barbre, D. Keng, V. Kolchenko, and S. Arnold, “Label-free detection of single protein using a nanoplasmonic-photonic hybrid microcavity,” Nano Letters **13**, 3347–3351 (2013). [CrossRef] [PubMed]

**4. **T. Lu, H. Lee, T. Chen, S. Herchak, J.-H. Kim, S. E. Fraser, R. C. Flagan, and K. Vahala, “High sensitivity nanoparticle detection using optical microcavities,” Proc. Nat. Acad. Sci. USA **108**, 5976–5979 (2011). [CrossRef] [PubMed]

**5. **Y. Sun, J. Liu, G. Frye-Mason, S.-J. Ja, A. K. Thompson, and X. Fan, “Optofluidic ring resonator sensors for rapid dnt vapor detection,” Analyst **134**, 1386–1391 (2009). [CrossRef] [PubMed]

**6. **G. Bahl, X. Fan, and T. Carmon, “Acoustic whispering-gallery modes in optomechanical shells,” New J. Phys. **14**, 115026 (2012). [CrossRef]

**7. **V. Ilchenko and A. Matsko, “Optical resonators with whispering-gallery modes-part II: applications,” IEEE J. Sel. Top. Quantum Electron. **12**, 15–32 (2006). [CrossRef]

**8. **M. Santiago-Cordoba, S. Boriskina, F. Vollmer, and M. Demirel, “Nanoparticle-based protein detection by optical shift of a resonant microcavity,” Appl. Phys. Lett. **99**, 073701 (2011). [CrossRef]

**9. **M. Borselli, T. Johnson, and O. Painter, “Beyond the Rayleigh scattering limit in high-Q silicon microdisks: theory and experiment,” Opt. Express **13**, 1515–1530 (2005). [CrossRef] [PubMed]

**10. **M. Hossein-Zadeh and K. J. Vahala, “Free ultra-high-Q microtoroid: a tool for designing photonic devices,” Opt. Express **15**, 166–175 (2007). [CrossRef] [PubMed]

**11. **D. Rowland and J. Love, “Evanescent wave coupling of whispering gallery modes of a dielectric cylinder,” Optoelectronics, IEEE Proceedings J **140**, 177–188 (1993). [CrossRef]

**12. **M. L. Gorodetsky and V. S. Ilchenko, “Optical microsphere resonators: optimal coupling to high-Q whispering-gallery modes,” J. Opt. Soc. Am. B **16**, 147–154 (1999). [CrossRef]

**13. **M. A. C. Shirazi, W. Yu, S. Vincent, and T. Lu, “Cylindrical beam propagation modelling of perturbed whispering-gallery mode microcavities,” Opt. Express **21**, 30243–30254 (2013). [CrossRef]

**14. **C.-G. Xu, X. Xiong, C.-L. Zou, X.-F. Ren, and G.-C. Guo, “Efficient coupling between dielectric waveguide modes and exterior plasmon whispering gallery modes,” Opt. Express **21**, 31253–31262 (2013). [CrossRef]

**15. **A. Kaplan, M. Tomes, T. Carmon, M. Kozlov, O. Cohen, G. Bartal, and H. G. L. Schwefel, “Finite element simulation of a perturbed axial-symmetric whispering-gallery mode and its use for intensity enhancement with a nanoparticle coupled to a microtoroid,” Opt. Express **21**, 14169–14180 (2013). [CrossRef] [PubMed]

**16. **M. Oxborrow, “Traceable 2-d finite-element simulation of the whispering-gallery modes of axisymmetric electromagnetic resonators,” IEEE Trans. Microw. Theory Tech. **55**, 1209–1218 (2007). [CrossRef]

**17. **J. Wiersig, “Boundary element method for resonances in dielectric microcavities,” J. Opt. A: Pure Appl. Opt. **5**, 53 (2003). [CrossRef]

**18. **Y.-F. Xiao, Y.-C. Liu, B.-B. Li, Y.-L. Chen, Y. Li, and Q. Gong, “Strongly enhanced light-matter interaction in a hybrid photonic-plasmonic resonator,” Phys. Rev. A **85**, 031805 (2012). [CrossRef]

**19. **I. Teraoka, S. Arnold, and F. Vollmer, “Perturbation approach to resonance shifts of whispering-gallery modes in a dielectric microsphere as a probe of a surrounding medium,” J. Opt. Soc. Am. B **20**, 1937–1946 (2003). [CrossRef]

**20. **M. R. Foreman and F. Vollmer, “Theory of resonance shifts of whispering gallery modes by arbitrary plasmonic nanoparticles,” New J. Phys. **15**, 083006 (2013). [CrossRef]

**21. **X. Du, S. Vincent, and T. Lu, “Full-vectorial whispering-gallery-mode cavity analysis,” Opt. Express **21**, 22012–22022 (2013). [CrossRef] [PubMed]

**22. **H. A. Haus, *Electromagnetic Noise and Quantum Optical Measurements* (Springer, 2000). [CrossRef]

**23. **I. H. Malitson, “Interspecimen comparison of the refractive index of fused silica,” J. Opt. Soc. Am. **55**, 1205–1208 (1965). [CrossRef]

**24. **M. Bass, C. DeCusatis, J. Enoch, V. Lakshminarayanan, G. Li, C. Macdonald, V. Mahajan, and E. Van Stryland, *Handbook of Optics, Volume I: Geometrical and Physical Optics, Polarized Light, Components and Instruments*, 3rd ed. (McGraw-Hill, 2010).

**25. **W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, *Numerical Recipes in C*, 2nd ed. (Cambridge University, 1992).

**26. **S. Spillane, T. J. Kippenberg, O. J. Painter, and K. J. Vahala, “Ideality in a fiber-taper-coupled microresonator system for application to cavity quantum electrodynamics,” Phys. Rev. Lett. **91**, 043902 (2003). [CrossRef] [PubMed]

**27. **H. Derudder, F. Olyslager, D. De Zutter, and S. Van Den Berghe, “Efficient mode-matching analysis of discontinuities in finite planar substrates using perfectly matched layers,” IEEE Trans. Antennas. Propag. **49**, 185–195 (2001). [CrossRef]