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
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 . Compared to boundary element [20,21,32–35], finite element , finite-difference time-domain (FDTD) , and free space radiation mode methods , 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  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  or the conformal mapping approach  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.
From the Helmhöltz equation, the electric field component E in a two-dimensional whispering-gallery microcavity under the scalar approximation satisfiesEq. (1): Eq. (3) at each angle ϕ for higher accuracy. From Eq. (5), we obtain the wave evolution along the azimuthal direction according to Fig. 1. Also, np,l is the refractive index of the waveguide structure at each point. Collecting ψp,l into a ket form |ψl〉 = (ψp0,l, ψp0+1,l,... ψp0+N,l)T and rearranging Eq. (7) into a matrix form, we obtain 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 ϕl to compute that at ϕ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) .
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)  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  was placed at the edge of the computation window. The PML is implemented by replacing the radial derivative withFig. 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 = 1020), 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 × 1016 in Fig. 2(b), where a minimal reflectivity of −50 dB is attained. Note that over a wide span of σ0 between 1016 and 2 × 1017 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 .
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 ), 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πmrP(ϕ = 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×106 towards 4.99×106 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 . 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 , where the resonance wavelength λres = 970.21 nm calculated via Richardson extrapolation is used as a reference. A least square fit to the relative error indicates a 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〉. ψ̂i is the normalized mode profile and ψ̂o 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 200th 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/Pin as a function of the number of rounds light circulated in the cavity. After circulating more than 25, 000 rounds, P/Pin saturates to 1.385 × 107 (blue curve) when resonance occurs. The saturation power is in excellent agreement with the theoretical prediction of (with a quality factor caculated by the mode solver of 5.36 × 106). 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 (λFWHM = (1 + 1/2Q)λres, 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 ρ̂ 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 , wherein the refractive index of the polystyrene bead was 1.576 + j(3.63 × 10−4) . We modelled the 3D microtoroid and bead structure in 2D using the effective index method (EIM) , 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. ) 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 . 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 asFig. 9 shows the field intensity for R90° = 37.5 μm while inset (b) corresponds to the extreme case of R90° = 32.5 μm. The time duration for Media 2 and Media 3, consisting of two rounds of propagation, is 3.5 minutes.
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
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]
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]
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]
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]
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]
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]
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]
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]