Cylindrical Beam Propagation Modelling of Perturbed Whispering-Gallery Mode Microcavities

We simulate light propagation in perturbed whispering-gallery mode microcavities using a two-dimensional finite-difference beam prop- agation 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.


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][2][3][4][5][6][7][8][9][10][11][12][13][14][15] In this work, we simulated 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.

Formulation
From the Helmholtz equation, the field component E z of the TE wave in a two-dimensional whispering-gallery microcavity satisfies where 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 r (ρ) is the normalized m th r azimuthal modal field distribution such that the squared norm |A| 2 represents the circulating power of the mode and m = m r + jm i is a complex constant whose real part m r represents the azimuthal mode order when the cavity is in resonance with the circulating light and its imaginary part m i characterizes the attenuation of the field in the azimuthal direction. Note that, in general, the real part m r can be any real number for a given wavelength λ . When a certain wavelength λ m yields an integer value for m r , resonance occurs.
In addition, multiple wavelengths may yield an identical integer m r where eigensolutionsψ m r correspond to resonance whispering-gallery modes with the same azimuthal order m r but different transverse modes. Both quantities can be obtained from the nonzero solution of the mode equation, in turn described as an eigenequation with eigenvalue m 2 derived from Eq. (1): If the aforementioned symmetry is broken due to an azimuthal angle dependent perturbation of the refractive index n(ρ, φ ) = n(ρ) + δ n(ρ, φ ) where the perturbation δ n(ρ, φ ) n(ρ), one may reformulate E z as wherem 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 It is necessary to point out that the choice ofm is arbitrary as long as SVEA holds; however, if the wavelength of the light is close to the resonance wavelength of the m th r order unperturbed WGM, it is convenient to selectm = m. We will therefore drop the bar in the rest of the text for convenience. Alternatively, one may treatm = 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 Discretizing the computation window uniformly so that the coordinates (ρ p , φ l ) of each grid (p, l) can be expressed as ρ p = p∆ρ, φ l = l∆φ , and ψ(ρ p , φ l ) = ψ p,l , one can evolve the field at φ l from a previous azimuthal angle φ l−1 according to Here ∆ρ and ∆φ are grid spacings along theρ andφ directions, as illustrated in Fig. 1. Also, n p,l is the refractive index of the waveguide structure at each point. Collecting ψ p,l into a ket form | ψ l = (ψ p 0 ,l , ψ p 0 +1,l , . . . ψ p 0 +N,l ) T and rearranging Eq. (7) into a matrix form, we obtaiñ whereH andD are two tridiagonal matrices. By adopting standard FD-BPM procedures, one may obtain the field evolution via Eq. (9) from the excitation field at l = 0.

Results and Discussions
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 ) [45, 46] at a wavelength of 970 nm and the surrounding water had a refractive index of 1.327 + j(3.37 × 10 −6 ) [47]. 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) [48] 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 [37] was placed at the edge of the computation window. The PML is implemented by replacing the radial derivative with in which σ (ρ) is defined as where ρ 0 is the inner edge of the PML layer and σ 0 is a constant. 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 reflected power was measured for different values of σ 0 . The different values of reflected power 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. The σ 0 value is 0 in Fig. 2(a), 2.5 × 10 16 (i.e. the optimal value) in Fig. 2(b), and 10 20 in Fig. 2(c). The optimal value was then utilized for the remaining simulations. In Fig. 3, 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 field distribution is plotted in Fig. 3. The resonance wavelength and the mode number calculated by the mode solver are 970.25 nm and 458 nm, respectively. The solid blue lines in Fig. 3 are showing the resonator edges and the dashed blue 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.
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 φ . Fig. 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, 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 the Q calculated by Richardson extrapolation as a reference value. As seen in Fig. 4(b) this 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 λ res = ∆Φλ m r , 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 define 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 200 th round trip.
In Fig. 7, we excited the cavity in continuous wave (CW) mode and plot the accumulated power normalized to the input power at the input of the cavity P/P in as a function of the number of rounds light circulated in the cavity. As illustrated, P/P in saturates to 1.95 × 10 6 (blue curve) when resonance occurs. The saturation power is in excellent agreement with the theoretic prediction of P sat /P in = Q 2 πm r that can be obtained from a mode solver with a relative error around 10 −6 after circulating more than 25, 000 rounds, as indicated in the red curve. We further plotted the total power when the input light wavelength had reached the full wave at half maximum point (λ FW HM = (1 + 1/2Q)λ res , i.e. the green curve). As expected, the power saturated to half of that corresponding to resonance.    Next, we applied the 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]. We modeled the 3D microtoroid and bead structure in 2D using the effective index method (EIM) [48]. 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. ∆λ res = ∆Φ 2πm λ 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 [49]. 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. 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. Fig. 9 shows the quality factor as a function of ellipticity, in which the ellipticity is defined as where R 0 = 45 µm is the fixed radius and R 90 is the variable radius for φ = 90 • . The inset in Fig. 9 shows the field intensity for the extreme case of R 90 = 32.5 µm.

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. Such properties are not readily obtainable through first-order perturbations on whispering-gallery modes as traditionally treated by other groups. Therefore, cylindrical BPM addresses the need for in-depth study of areas such as biosensing with WGM's where azimuthal symmetry is perturbed. A full vectorial three-dimensional BPM approach will be investigated in future research.