## Abstract

We simulate light propagation in perturbed whispering-gallery mode microcavities using a two-dimensional finite-difference beam propagation method in a cylindrical coordinate system. Optical properties of whispering-gallery microcavities perturbed by polystyrene nanobeads are investigated through this formulation. The light perturbation as well as quality factor degradation arising from cavity ellipticity are also studied.

© 2013 Optical Society of America

## 1. Introduction

Whispering-gallery mode (WGM) microcavities are at the frontier of research on subjects ranging from biosensing, nonlinear optics, and laser physics, to fundamental physics such as cavity quantum electrodynamics [1–15]. Contrary to its rapid experimental advances, numerical exploration of WGM’s has been largely lagging behind with a limited number of available options [16–21]. On the other hand, the Beam Propagation Method (BPM) has a long history [22–30] in modelling light propagation along both straight and curved waveguides as well as whispering-gallery microcavity eigenmode analyses [31]. Compared to boundary element [20,21,32–35], finite element [36], finite-difference time-domain (FDTD) [37], and free space radiation mode methods [38], BPM remains highly efficient without sacrificing substantial accuracy. By adopting a perfectly matched layer (PML) [39,40] to absorb the light which is otherwise reflected at the computation window and following the procedure formulated in [41] to correct inaccuracies incurred at the refractive index discontinuities in high refractive index contrast waveguide structures, the Finite-Difference Beam Propagation Method (FD-BPM) can achieve high accuracy with a rapid convergence rate. Conventional FD-BPM formulations are based on the Fresnel approximation, where light is assumed to propagate close to the propagation axis [24–26, 28]. To overcome this limitation for bent waveguide modelling, high-order algorithms known as wide-angle BPM [42] or the conformal mapping approach [43] are desirable. Alternatively, BPM may be reformulated in cylindrical coordinates systems to analyze such structures [44–46].

In this work we simulated, within MATLAB, the light propagation in a WGM microcavity by implementing FD-BPM in a cylindrical system as shown in Fig. 1. The field perturbation from a nanobead attached to the microcavity and quality factor degradation arising from cavity deformations were investigated. The computed field distribution correctly includes the radiative field component, which a mode analysis technique would fail to simulate.

## 2. Formulation

From the Helmhöltz equation, the electric field component *E* in a
two-dimensional whispering-gallery microcavity under the scalar approximation satisfies

*n*(

*ρ*,

*ϕ*) is the refractive index of the cavity,

*k*

_{0}= 2

*π/λ*is the wave number in free space, and

*λ*is the vacuum wavelength of the light circulating in the cavity. For a perfect WGM cavity with azimuthal symmetry, the refractive index is independent of

*ϕ*(

*n*(

*ρ*,

*ϕ*) =

*n*(

*ρ*)). The electric field can be approximated as propagating in the form of where

*m*=

*m*+

_{r}*jm*is a complex constant whose real and imaginary parts are

_{i}*m*and

_{r}*m*, respectively. Note that, in general, the real part

_{i}*m*can be any real number for a given wavelength

_{r}*λ*. When a certain wavelength

*λ*

_{mr}yields an integer value for

*m*, resonance occurs. This integer value

_{r}*m*represents the azimuthal order of a cavity resonance mode while the imaginary part

_{r}*m*characterizes the attenuation of the field in the azimuthal direction.

_{i}*ψ̂*

_{mr}(

*ρ*) is the normalized ${m}_{r}^{\text{th}}$ azimuthal modal field distribution such that the squared norm |

*A*|

^{2}represents the circulating power of the mode. In addition, multiple wavelengths may yield an identical integer

*m*where

_{r}*ψ̂*

_{mr}correspond to resonance whispering-gallery modes with the same azimuthal order

*m*but different transverse modes. Both quantities can be obtained from the nonzero solution of the mode equation, in turn described as an eigenequation with the complex eigenvalue

_{r}*m*

^{2}and eigenvector

*ψ̂*

_{mr}derived from Eq. (1):

*n*(

*ρ*,

*ϕ*) =

*n*(

*ρ*) +

*δn*(

*ρ*,

*ϕ*) where the perturbation

*δn*(

*ρ*,

*ϕ*) ≪

*n*(

*ρ*), one may reformulate

*E*as where

*m̄*is a reference value such that

*ψ*(

*ρ*,

*ϕ*) varies slowly along the azimuthal direction or, equivalently, the slowly varying envelope approximation (SVEA) holds. This is mathematically written as

*m̄*is arbitrary as long as SVEA holds; however, if the wavelength of the light is close to the resonance wavelength of the ${m}_{r}^{\text{th}}$ order unperturbed WGM, it is convenient to select

*m̄*=

*m*. We will therefore drop the bar in the rest of the text for convenience. Alternatively, one may treat

*m̄*=

*m*(

*ϕ*) as a

*ϕ*-dependent quantity where

*m*(

*ϕ*) is obtained from solving Eq. (3) at each angle

*ϕ*for higher accuracy. From Eq. (5), we obtain the wave evolution along the azimuthal direction according to

*ρ*,

_{p}*ϕ*) of each grid (

_{l}*p*,

*l*) can be expressed as

*ρ*=

_{p}*p*Δ

*ρ*,

*ϕ*=

_{l}*l*Δ

*ϕ*, and

*ψ*(

*ρ*,

_{p}*ϕ*) =

_{l}*ψ*, one can evolve the field at

_{p,l}*ϕ*from a previous azimuthal angle

_{l}*ϕ*

_{l}_{−1}according to

*ρ*and Δ

*ϕ*are grid spacings along the

*ρ̂*and

*ϕ̂*directions, as illustrated in Fig. 1. Also,

*n*is the refractive index of the waveguide structure at each point. Collecting

_{p,l}*ψ*into a ket form |

_{p,l}*ψ*〉 = (

_{l}*ψ*

_{p0,l},

*ψ*

_{p0+1,l},...

*ψ*

_{p0+N,l})

*and rearranging Eq. (7) into a matrix form, we obtain where*

^{T}*H̃*and

*D̃*are two tridiagonal matrices. By adopting standard FD-BPM procedures [47], one may obtain the field evolution via Eq. (9) from the excitation field at

*l*= 0. In stark contrast to a boundary element method (BEM, as reported in [32–35]) that requires the field of the entire structure for successive azimuthal angle steps, this expression solely requires the field at

*ϕ*to compute that at

_{l}*ϕ*

_{l}_{+1}. Consequently, the efficiency of BPM is orders of magnitude greater than that of BEM and is capable of solving large-scale cavity structures that exceed BEM’s capabilities. Unlike BEM, BPM may also model cavity refractive index profiles that are not necessarily piecewise homogeneous.

## 3. Results and discussion

To characterize the BPM, we first tested it on a perfect silica microring resonator immersed in
water. The refractive index of the silica ring was 1.4508 + *j*(7.11
× 10^{−12}) [48, 49] at a wavelength of 970 nm and the surrounding water had
a refractive index of 1.327 + *j*(3.37 × 10^{−6})
[50].

The resonator had a 45-*μm* major radius and a
10-*μm* minor diameter. To simplify the analysis, we reduced the
three-dimensional waveguide structure to a two-dimensional one through the use of an effective index
method (EIM) [51] along the
*z*-direction. The cavity was excited by a modal field obtained from the mode solver.
To minimize the spurious reflection at the edges of the computation window, a
4-*μm* PML [39] was
placed at the edge of the computation window. The PML is implemented by replacing the radial
derivative with

*σ*(

*ρ*) is defined as

*ρ*

_{0}is the inner edge of the PML layer and

*σ*

_{0}is a wavelength independent constant characterized by the light attenuation strength within the PML. To optimize the PML performance and determine the optimal value for

*σ*

_{0}, a simple experimental simulation was conducted. An optical beam was launched towards a 2-

*μm*PML through a straight waveguide and the power reflectivity was measured for different values of

*σ*

_{0}. The different values of reflectivity in dB versus

*σ*

_{0}are presented in Fig. 2. The insets show the field intensity for three different

*σ*

_{0}values, wherein the white lines are the inner PML edges. Figure 2(a) illustrates an extreme case where the light reaches the outer edge of an undamped PML (i.e.

*σ*

_{0}= 0) and is reflected back. In the case of an overdamped PML (i.e.

*σ*

_{0}= 10

^{20}), as illustrated in Fig. 2(c), light is reflected at the PML’s inner edge. The optimum value of

*σ*

_{0}is found at 2.5 × 10

^{16}in Fig. 2(b), where a minimal reflectivity of −50 dB is attained. Note that over a wide span of

*σ*

_{0}between 10

^{16}and 2 × 10

^{17}the reflectivity at the computation window edge is kept below −30 dB. The optimal value was then utilized for the remaining simulations given that

*σ*

_{0}was insensitive to, for example, the incident angle at the PML edge, the wavelength, and the waveguide structures inside the computation window [40].

In Fig. 3(a), we plot the intensity distribution in
logarithmic scale by setting a 36-*μm* window size, 1601 grids in the
*ρ̂* direction, and 3000 grids in the
*ϕ̂* direction for 2*π* radians, where we
excited the ring with its fundamental mode. The PML thickness at the inside and outside boundaries
is 5 *μm*. The resonance wavelength and the real as well as imaginary parts
of the mode number calculated by the mode solver are respectively 970.25 nm, 458, and 4.5 ×
10^{−5}. The solid green lines in Fig. 3(a)
are showing the resonator edges and the dashed green lines define the PML boundary. As can be seen
in this figure, a small portion of the energy radiates towards the computation window edge yet the
adopted PML efficiently prevents the otherwise spurious reflection from returning to the resonator.
Figure 3(b) provides a portion of the final intensity
distribution frame from the supplementary video clip, Media
1 (included in an external, high-resolution video clip that combines
this article’s multimedia content [52]), elucidating circular field propagation for 1/4 of one round trip. Note that
simulation of the 1-mm diameter microdisk in Fig. 3(b)
spanned 23 seconds on a desktop computer equipped with a 3.10 GHz AMD FX-8120 8-core processor and
32 GB of memory. In contrast, a simulation of such a large cavity has yet to be reported from other
established methods, such as BEM, finite element method (FEM), or FDTD according to the
authors’ knowledge.

We further computed the quality factor *Q* for different grid schemes according to
*Q* = 2*πm _{r}P*(

*ϕ*= 0)/(

*P*(

*ϕ*= 0) −

*P*(

*ϕ*= 2

*π*)).

*P*(

*ϕ*) is the total power at an azimuthal angle

*ϕ*. Figure 4(a) is the plot of quality factor vs. radial and azimuthal grid spacings. The computation window is set to 20

*μm*in the

*ρ̂*direction and

*π*in the

*ϕ̂*direction. As is depicted, the

*Q*converges from 5.4×10

^{6}towards 4.99×10

^{6}by reducing the grid spacing from 200 nm to 3.1 nm in the radial direction and 0.196 radians to 0.0015 radians in the azimuthal direction. We also computed the relative error by adopting a reference

*Q*at infinitesimal grid spacing predicted by the established Richardson extrapolation approach [53]. A set of

*Q*values were computed at different grid spacings, then the expected value at infinitesimal grid spacing was obtained by treating the data set as a function of grid spacing for a least square fit. As seen in Fig. 4(b) the relative error of

*Q*was reduced to 5 × 10

^{−4}from 8 × 10

^{−2}, accordingly. Furthermore, in Fig. 4(c) we plot the

*Q*as well as the relative error as a function of azimuthal grid spacing Δ

*ϕ*by setting Δ

*ρ*= 3.1 nm. In Fig. 4(d) the same quantities are plotted as a function of radial grid spacing Δ

*ρ*. The least square fits to the relative error curves indicate a convergence rate of 0.9 in the azimuthal direction and 2.8 in the radial direction. Clearly, this suggests that faster convergence is attainable by adopting higher-order finite-difference schemes.

In Fig. 5, we plot the resonance wavelength (blue cross
markers) and corresponding relative error (red cross markers) vs. radial grid spacing
Δ*ρ* by setting Δ*ϕ* = 0.003
radians. Here, the resonance wavelength *λ _{res}* is obtained from
the round trip total phase shift ΔΦ of the electrical field computed from the
wavelength

*λ*adopted in the BPM according to ${\lambda}_{\mathit{res}}=\frac{\mathrm{\Delta}\mathrm{\Phi}\lambda}{{m}_{r}}$, where the resonance wavelength

*λ*= 970.21 nm calculated via Richardson extrapolation is used as a reference. A least square fit to the relative error indicates a

_{res}*O*(Δ

*ρ*) convergence rate.

To further demonstrate the validity of the formulation, we launched an arbitrary Gaussian field
at the input of the same structure as illustrated in the insets (red curve) of Fig. 6. For comparison, we also plotted the fundamental mode profile obtained
by the mode solver and displayed it as the blue curve in the insets. As shown in the insets of Fig. 6, after propagating 1, 25, and 125 rounds in the resonator
from left to right, the circulating field distribution gradually evolves into the mode profile. To
quantitatively analyze the field evolution, we de-fine a normalized overlap factor
Γ̂ = 〈*ψ̂ _{o}* |

*ψ̂*〉.

_{i}*ψ̂*is the normalized mode profile and

_{i}*ψ̂*is the normalized output field profile after each round trip of propagation. The magnitude of the overlap factor and its departure from unity are respectively plotted in Fig. 6 as blue and red curves. As is shown, the overlap factor reaches 0.99 at the 200

_{o}^{th}round trip. In Fig. 7(a), we excited the cavity with a continuous wave (CW) mode and plotted the accumulated power normalized to the input power at the cavity’s input cross section

*P/P*as a function of the number of rounds light circulated in the cavity. After circulating more than 25, 000 rounds,

_{in}*P/P*saturates to 1.385 × 10

_{in}^{7}(blue curve) when resonance occurs. The saturation power is in excellent agreement with the theoretical prediction of ${P}_{\mathit{sat}}/{P}_{\mathit{in}}={\left(\frac{Q}{\pi {m}_{r}}\right)}^{2}=1.388\times {10}^{7}$ (with a quality factor caculated by the mode solver of 5.36 × 10

^{6}). The relative deviation of the accumulated power from the saturation power (i.e. the red curve) is depicted as well. We further plotted the total power when the input light wavelength had reached the full wave at half maximum point (

*λ*= (1 + 1/2

_{FWHM}*Q*)

*λ*, i.e. the green curve). As expected, the power saturated to half of that corresponding to resonance. In Fig. 7(b), the normalized saturation power and its relative error are plotted as a function of grid size in the

_{res}*ρ̂*direction. Next, we applied BPM to whispering-gallery mode nanodetection modelling [4, 5]. Here we simulated the induced resonance wavelength shift caused by single polystyrene beads binding to the surface of a silica microtoroid immersed in water [5], wherein the refractive index of the polystyrene bead was 1.576 +

*j*(3.63 × 10

^{−4}) [54]. We modelled the 3D microtoroid and bead structure in 2D using the effective index method (EIM) [51], placing the nanoparticle on the micro-toroid’s equator. The motive for this is to maximize the resonance wavelength shift; however, it is possible to choose other geometries. Note that the 2D refractive index profile is not piecewise homogeneous and thus cannot be modelled with a boundary element method, such as that of [20,21,32–35]. The field evolution across a 400-nm radius bead attached to the cavity is displayed in Fig. 8(a), while a zoomed-in plot of field distribution around the bead is displayed as an inset. We further obtained the resonance wavelength shift from the additional phase shift of the electrical field attributed to the bead (i.e. $\mathrm{\Delta}{\lambda}_{\mathit{res}}=\frac{\mathrm{\Delta}\mathrm{\Phi}}{2\pi m}{\lambda}_{\mathit{res}}$) and plotted them as a function of bead radius in Fig. 8(b) as red cross markers. For comparison, we also plotted the corresponding shift predicted by the perturbation method [55]. As shown, both are in good agreement aside from the fact that there is a identifiable departure due to the 2D simplification with EIM. We believe that a three-dimensional full wave beam propagation method should model the shift with higher accuracy. Finally, we applied the BPM to model the ellipticity effect and hence emulate a nonideal cavity due to fabrication imperfections. Setting the major radius at

*ϕ*= 0 to 45

*μm*and sweeping the other radius from 32.5

*μm*to 55

*μm*, we calculated the quality factor as explained before using BPM. The window size is set to 50

*μm*with 2201 grids in the

*ρ̂*direction and 2

*π*radians with 3000 grids in the

*ϕ̂*direction. We excited the ring with its fundamental mode for a 45-

*μm*major radius and 10-

*μm*minor diameter ring resonator. Figure 9 shows the quality factor as a function of ellipticity, in which the ellipticity is defined as

*R*

_{0°}= 45

*μm*is the fixed radius and

*R*

_{90°}is the variable radius for

*ϕ*= 90°. Inset (a) in Fig. 9 shows the field intensity for

*R*

_{90°}= 37.5

*μm*while inset (b) corresponds to the extreme case of

*R*

_{90°}= 32.5

*μm*. The time duration for Media 2 and Media 3, consisting of two rounds of propagation, is 3.5 minutes.

## 4. Conclusion

In conclusion, we implemented a 2D Finite-Difference Beam Propagation Method for simulating light propagation in whispering-gallery mode microcavities. With this method, we demonstrate the calculation of key optical properties such as resonance wavelength, quality factor in cases where the azimuthal symmetry is absent due to singular perturbations from nano particles, and ellipticity of the cavity. The field scattering that arises from asymmetry is clearly visible from our simulation. The BPM is capable of modelling whispering-gallery mode cavities ranging from micron to millimeter or even meter scale with no additional computational resource requirements. Therefore, cylindrical BPM addresses the need for in-depth study of areas such as biosensing with WGM’s where azimuthal symmetry is perturbed. The presented 2D FD-BPM is capable of simulating other deformations like those for quadrupoles, spirals, stadiums, or quasistadiums [56, 57]. Other deformation effects such as cavity emission and wave-chaos can also be treated by this method. Numerical techniques that are used for WGM cavity deformation simulations have limitations on the maximum cavity size, while such a limitation is of the least concern for BPM. Research regarding full-vectorial three-dimensional BPM, PML optimization in the presented simulation domain, and wave chaos analysis is forthcoming.

## References and links

**1. **M. L. Gorodetsky, A. A. Savchenkov, and V. S. Ilchenko, “Ultimate Q of optical microsphere
resonators,” Opt. Lett. **21**, 453–455 (1996). [CrossRef] [PubMed]

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

**3. **F. Vollmer, D. Braun, A. Libchaber, M. Khoshsima, I. Teraoka, and S. Arnold, “Protein detection by optical shift of a resonant
microcavity,” Appl. Phys. Lett. **80**, 4057–4059 (2002). [CrossRef]

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

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

**6. **J. Dominguez-Juarez, G. Kozyreff, and J. Martorell, “Whispering gallery microresonators for second harmonic
light generation from a low number of small molecules,” Nat.
Commun. **2**, 1–8 (2010).

**7. **J. Knittel, T. G. McRae, K. H. Lee, and W. P. Bowen, “Interferometric detection of mode splitting for whispering
gallery mode biosensors,” Appl. Phys. Lett. **97**, 1–3 (2010). [CrossRef]

**8. **Y. Sun and X. Fan, “Optical ring resonators for biochemical and chemical
sensing,” Anal. Bioanal.Chem. **399**, 205–211 (2011). [CrossRef]

**9. **S. Spillane, T. J. Kippenberg, and K. J. Vahala, “Ultralow-threshold raman laser using a spherical dielectric
mcirocavity,” Nature **415**, 621–623 (2002). [CrossRef] [PubMed]

**10. **B. Min, T. J. Kippenberg, and K. J. Vahala, “Compact, fiber-compatible, cascaded raman
laser,” Opt. Lett. **28**, 1507–1509 (2003). [CrossRef] [PubMed]

**11. **M. Cai and K. J. Vahala, “Highly efficient hybrid fiber taper coupled microsphere
laser,” Opt. Lett. **26**, 884–886 (2001). [CrossRef]

**12. **A. Polman, B. Min, J. Kalkman, T. J. Kippenberg, and K. J. Vahala, “Ultralow-threshold erbium-implanted toroidal microlaser on
silicon,” Appl. Phys. Lett. **84**, 1037–1039 (2004). [CrossRef]

**13. **T. Lu, L. Yang, R. V. A. van Loon, A. Polman, and K. J. Vahala, “On-chip green silica upconversion
microlaser,” Opt. Lett. **34**, 482–484 (2009). [CrossRef] [PubMed]

**14. **T. Lu, L. Yang, T. Carmon, and B. Min, “A narrow-linewidth on-chip toroid raman
laser,” IEEE J. Quantum Electron. **47**, 320–326 (2011). [CrossRef]

**15. **S. M. Spillane, T. J. Kippenberg, K. J. Vahala, K. W. Goh, E. Wilcut, and H. J. Kimble, “Ultrahigh-Q toroidal microresonators for cavity quantum
electrodynamics,” Phys. Rev. A At. Mol. Opt. Phys. **71**, 013817 (2005). [CrossRef]

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

**17. **J. K. S. Poon, J. Scheuer, S. Mookherjea, G. T. Paloczi, Y. Huang, and A. Yariv, “Matrix analysis of microring coupled-resonator optical
waveguides,” Opt. Express **12**, 90–103 (2004). [CrossRef] [PubMed]

**18. **J. Hong, W. P. Huang, and T. Makino, “On the transfer matrix method for distributed-feedback
waveguide devices,” J. Lightwave Technol. **10**, 1860–1868 (1992). [CrossRef]

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

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

**21. **C.-L. Zou, H. G. L. Schwefel, F.-W. Sun, Z.-F. Han, and G.-C. Guo, “Quick root searching method for resonances of dielectric
optical microcavities with the boundary element method,” Opt.
Express **19**, 15669–15678 (2011). [CrossRef] [PubMed]

**22. **M. D. Feit and J. J. A. Fleck, “Light propagation in graded-index optical
fibers,” Appl. Opt. **17**, 3990–3998 (1978). [CrossRef] [PubMed]

**23. **D. Yevick and B. Hermansson, “Efficient beam propagation
techniques,” IEEE J. Quantum Electron. **26**, 109–112 (1990). [CrossRef]

**24. **J. Saijonmaa and D. Yevick, “Beam-propagation analysis of loss in bent optical
waveguides and fibers,” J. Opt. Soc. Am. **73**, 1785–1791 (1983). [CrossRef]

**25. **W. Huang, C. Xu, S.-T. Chu, and S. K. Chaudhuri, “The finite-difference vector beam propagation method:
Analysis and assessment,” J. Lightwave Technol. **10**, 295–305 (1992). [CrossRef]

**26. **J. V. Roey, J. van der Donk, and P. E. Lagasse, “Beam-propagation method: analysis and
assessment,” J. Opt. Soc. Am. **71**, 803–810 (1981). [CrossRef]

**27. **W. Huang, C. Xu, and S. Chaudhuri, “A finite-difference vector beam propagation method for
three-dimensional waveguide structures,” IEEE Photonics Technol.
Lett. **4**, 148–151 (1992). [CrossRef]

**28. **B. Hermansson and D. Yevick, “Propagating-beam-method analysis of two-dimensional
microlenses and three-dimensional taper structures,” Opt. Soc. Am.
A **1**, 663–671 (1984). [CrossRef]

**29. **H. Rao, R. Scarmozzino, and R. M. Osgood, “A bidirectional beam propagation method for multiple
dielectric interfaces,” IEEE Photonics Technol. Lett. **11**, 830–832 (1999). [CrossRef]

**30. **R. Scarmozzino, A. Gopinath, R. Pregla, and S. Helfert, “Numerical techniques for modeling guided-wave photonic
devices,” IEEE J. Sel. Top. Quantum Electron. **6**, 150–162 (2000). [CrossRef]

**31. **A. Ahmed, R. Koya, O. Wada, M. Wang, and R. Koga, “Eigenmode analysis of whispering gallery mode of
pillbox-type optical resonators utilizing the FE-BPM formulation,”
IEICE Trans. Electron. **E78-C**, 1638–1645
(1995).

**32. **W. Yang and A. Gopinath, “A boundary integral method for propagation problems in
integrated optical structures,” IEEE Photonics Technol.
Lett. **7**, 777–779 (1995). [CrossRef]

**33. **T. Lu and D. Yevick, “Boundary element analysis of dielectric
waveguides,” J. Opt. Soc. Am. A **19**, 1197–1206 (2002). [CrossRef]

**34. **T. Lu and D. Yevick, “Comparative evaluation of a novel series approximation for
electromagnetic fields at dielectric corners with boundary element method
applications,” J. Lightwave Technol. **22**, 1426–1432 (2004). [CrossRef]

**35. **T. Lu and D. Yevick, “A vectorial boundary element method analysis of integrated
optical waveguides,” J. Light-wave Technol. **21**, 1793–1807 (2003). [CrossRef]

**36. **H. Deng and D. Yevick, “The nonunitarity of finite-element beam propagation
algorithms,” IEEE Photonics Technol. Lett. **17**, 1429–1431 (2005). [CrossRef]

**37. **S.-T. Chu and S. Chaudhuri, “A finite-difference time-domain method for the design and
analysis of guided-wave optical structures,” J. Lightwave
Technol. **7**, 2033–2038 (1989). [CrossRef]

**38. **M. Reed, T. M. Benson, P. C. Kendall, and P. Sewell, “Antireflection-coated angled facet
design,” Proc. Inst. Electr. Eng. **143**, 214–220
(1996).

**39. **J.-P. Berenger, “A perfectly matched layer for the absorption of
electromagnetic waves,” J. Comput. Phys. **114**, 185–200 (1994). [CrossRef]

**40. **W. P. Huang, C. L. Xu, W. Lui, and K. Yokoyama, “The perfectly matched layer (PML) boundary condition for
the beam propagation method,” IEEE Photonics Technol. Lett. **8**, 649–651 (1996). [CrossRef]

**41. **G. R. Hadley, “High-accuracy finite-difference equations for dielectric
waveguide analysis I: Uniform regions and dielectric interfaces,” J.
Lightwave Technol. **20**, 1210–1218 (2002). [CrossRef]

**42. **G. R. Hadley, “Wide-angle beam propagation using pade approximant
operators,” Opt. Lett. **17**, 1426–1428 (1992). [CrossRef] [PubMed]

**43. **S. Lidgate, P. Sewell, and T. Benson, “Conformal mapping: limitations for waveguide bend
analysis,” IEE Proc. Sci. Meas. Technol. **149**, 262–266 (2002). [CrossRef]

**44. **M. Rivera, “A finite difference BPM analysis of bent dielectric
waveguides,” J. Lightwave Technol. **13**, 233 (1995). [CrossRef]

**45. **H. Deng, G. H. Jin, J. Harari, J. P. Vilcot, and D. Decoster, “Investigation of 3-D semivectorial finite-difference beam
propagation method for bent waveguides,” J. Lightwave
Technol. **16**, 915–922 (1998). [CrossRef]

**46. **M. Krause, “Finite-difference mode solver for curved waveguides with
angled and curved dielectric interfaces,” J. Lightwave
Technol. **29**, 691–699 (2011). [CrossRef]

**47. **K. Kawano and T. Kitoh, *Introduction to Optical Waveguide Analysis* (John
Wiley, 2001).

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

**49. **R. Kitamura, L. Pilon, and M. Jonasz, “Optical constants of silica glass from extreme ultraviolet
to far infrared at near room temperature,” Appl. Opt. **46**, 8118–8133 (2007). [CrossRef] [PubMed]

**50. **G. M. Hale and M. R. Querry, “Optical constants of water in the 200-nm to 200-m
wavelength region,” Appl. Opt. **12**, 555–563 (1973). [CrossRef] [PubMed]

**51. **K. S. Chiang, “Performance of the effective-index method for the analysis
of dielectric waveguides,” Opt. Lett. **16**, 714–716 (1991). [CrossRef] [PubMed]

**52. **M. A. C. Shirazi, W. Yu, S. Vincent, and T. Lu, “Whispering-gallery mode propagation
simulations,” http://youtu.be/SJpEIkmsfMs (2013).

**53. **W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, *Numerical Recipes in C*, 2
(Cambridge University, 1992), chap. 16, pp.
718–725.

**54. **X. Ma, J. Q. Lu, R. S. Brock, K. M. Jacobs, P. Yang, and X.-H. Hu, “Determination of complex refractive index of polystyrene
microspheres from 370 to 1610 nm,” Phys. Med. Biol. **48**, 4165–4172 (2003). [CrossRef]

**55. **M. A. Waldron, “Perturbation theory of resonant
cavities,” Proc. IEE Part C Monogr. **107**, 272–274 (1960). [CrossRef]

**56. **H. G. L. Schwefel, H. E. Tureci, D. A. Stone, and R. K. Chang, “Progress in asymmetric resonant cavities: Using shape as a
design parameter in dielectric microcavity lasers,” in *Optical
Microcavities*, K. Vahala, ed. (World Scientific,
2005).

**57. **Y.-F. Xiao, C.-L. Zou, Y. Li, C.-H. Dong, Z.-F. Han, and Q. Gong, “Asymmetric resonant cavities and their applications in
optics and photonics: a review,” Front. Optoelectron. China **3**, 109–124 (2010). [CrossRef]