1. Introduction

Whispering-gallery-mode (WGM) cavities have been adopted in a wide range of applications [1–8]. By confining light energy within a mode volume on the order of hundreds to thousands of cubic microns [9, 10], WGM resonators [1, 11] can be used to sense minute local environment changes. In particular, nanoparticles adsorbed to the cavity surface will be polarized by the WGM and cause a shift of cavity resonance wavelength [12, 13]. In addition, it has been recently reported that a localized surface plasmon resonance (LSPR) can largely enhance the sensitivity of a microcavity [14–17] by increasing the field intensity around bound plasmonic particles.

To guide research in the aforementioned fields, numerical techniques are required. Currently, a perturbative approach [18, 19] that employs two-dimensional mode solvers such as those implemented in the finite element method [20] can be adapted to predict particle induced resonance shifts of an axisymmetric WGM. This technique involves the discretization of the two-dimensional cross section of the unperturbed WGM along the transverse direction to obtain a modal field distribution while the resulting resonance wavelength shift, e.g. from nanobead adsorption events, is computed as a perturbation to the mode. This yields highly accurate results for plain nanoparticles that do not trigger significant distortion to the cavity modal field. In the event that gold nanoparticles land on the surface of a WGM cavity, the associated plasmonic effects focus the light to the proximity of the bead, hence yielding larger inaccuracy. Such inaccuracy invalidates the perturbation approach in the case of large metallic beads or at a wavelength close to the plasmon resonance [19]. To circumvent this problem, a refinement scheme for the first order perturbation approach was proposed [15, 21] to model the WGM cavities with bound plasmonic particles. By approximating the WGM using a polarized plane wave, such a scheme employs an additional optical scattering solver to refine the field surrounding the bound particles. While such a refinement scheme provides an effective solution to nanoscale plasmonic and non-plasmonic binding problems, the plane wave approximation approaches its limit for larger particles when the field across the particle is no longer uniform. Further, perturbative approaches cannot simulate the resonance wavelength and quality factor of a WGM cavity that does not exhibit axisymmetry under the current formalism [22, 23]. Alternatively, one may choose first-principle techniques that directly solves Maxwell’s equations, Helmholtz equations or Green functions with numerical techniques such as finite difference, finite element or boundary element [24–30] methods to perform a full wave analysis over the three-dimensional WGM microcavities. Unfortunately, such techniques require the discretization of the WGM along all three spatial dimensions. This necessitates extensive computational resources compared to 2-D mode solver based techniques due to orders of magnitude increase in grid points. To accommodate such difficulties, first principle techniques are commonly used in conjunction with the effective index technique to approximate the three-dimensional WGM as a two-dimensional one. In [25], a new approach is proposed to simulate a WGM with a localized perturbation by treating the perturbed section of the WGM as a three-dimensional closed cavity, where artificial boundary conditions such as a perfect electrical wall are imposed at the input and output ends of the section. A three-dimensional mode of this perturbed section is then computed with the finite element method to approximate the distorted WGM field and other corresponding parameters. Such a method can yield better accuracy compared to a perturbation approach by applying a correction to the WGM field distorted by the bead. In this case, the number of grid points per computation step that are involved is larger than the perturbation approach yet fewer than first principle techniques. On the other hand the boundary condition imposed in [25], as pointed out by the authors of that article, is based on the assumption that the light scattered by the adsorbed bead will quickly ”heal” within the perturbed section and that the model ”represents a resonator with many particles [located periodically] along the circumference.” Such an assumption no longer holds for the binding of larger beads or for beads that are located at random positions, thereupon introducing significant phase shift and amplitude distortion of the otherwise unperturbed WGM field. The field obtained from this method is a standing wave pattern of the cavity, which departs from the realistic situation where a travelling wave is involved. Furthermore, the multi-scale nature of a WGM-nanoparticle system introduces additional challenges with these methods. For a large size perturbation such as a fiber taper placed close to the cavity, a large section is needed for the model. This causes the resulting mesh to become too large for a conventional computer to handle. Here, we modify the mode-matching method (MMM) to model WGM cavities and demonstrate its capability to address all the challenges mentioned above.

MMM has been widely used in modeling and designing straight waveguide discontinuities, junctions, and filters in Cartesian coordinates [31–34]. In MMM, the optical structure is divided into a number of slices along the propagation axis. In each slice a complete orthonormal set of eigen modes derived from an infinitely long waveguide, whose cross section along the propagation direction has an identical refractive index profile to the slice, are employed to form the basis of the solution space where an arbitrary electromagnetic (EM) field distribution can be expanded upon. The field will then be projected to the next slice along the propagation direction as if each modal field propagates with its distinct propagation constant. The new field will then be re-expanded on the eigen modes native to the new slice. The wave propagation behaviour can be simulated by progressively projecting fields via the above-mentioned approach from one slice to the next along the propagation direction. By analogy, we partition a WGM system along the azimuthal direction and obtain full vector WGM modes at each slice. Mode matching between neighbouring slices is then performed to emulate wave propagation in the WGM cavity, while optical properties such as resonance wavelength and quality factor are calculated accordingly.

2. Theoretical formulations

2.1. Ideal whispering-gallery-mode microcavities

An ideal or axisymmetric whispering-gallery-mode microcavity has a complex refractive profile ñ(ρ, z, ϕ) independent of the azimuthal angle ϕ (ie. ñ(ρ, z, ϕ) = ñ(ρ, z)) in a concentric cylindrical coordinate system, as shown in Fig. 1(a). The electric field Ẽ(ρ, z, ϕ) of light at wavelength λ₀ that propagates inside the cavity has the form

 

Fig. 1: (a) Light circulating along the azimuthal direction in a whispering-gallery-mode microcavity (e.g. a silica microtoroid). A cylindrical coordinate system is used for modelling purposes. (b) Light propagating from ϕ to ϕ + δϕ as it passes by a bound particle.

Download Full Size | PDF

2.2. Non-ideal whispering-gallery-mode microcavities

In general, a whispering-gallery-mode microcavity may have a ϕ-dependent refractive index profile n(ρ, z, ϕ) if particle binding occurs, surface imperfections exist, a waveguide taper is placed close to the cavity, or other geometrical disturbances are introduced. For simplicity, we drop ∼ symbols for all the physical quantities to distinguish them from those defined in an ideal WGM.

At an azimuthal angle ϕ₀, the electrical field distribution at the cross section E(ρ, z, ϕ₀) can be expanded onto the normalized WGM mode ê(ρ, z, ϕ₀) at ϕ₀ according to

By equating the right hand sides of (6) and (8), multiplying both sides by $\frac{1}{2{\eta }_0}{n}_r\left(\rho ,z,{\varphi }_0+\delta \varphi \right)\widehat{\mathbf{e}}\left(\rho ,z,{\varphi }_0+\delta \varphi \right)$, and integrating over the cross section at ϕ₀ + δϕ, we obtain the evolution of A(ϕ₀ + δϕ) according to

We may then obtain an additional loss term m_m(ϕ₀) characterized by

To satisfy the resonance condition, the overall phase change of the field after propagating through a 2π azimuthal angle should be 2Mπ:

Together with (7) we obtain the resonance wavelength of the whispering-gallery-mode microcavity λ_res by simply taking the arithmetic mean of the resonance wavelengths λ_r(ϕ) at each cross section of the cavity, where λ_r(ϕ₀) at each ϕ₀ can be obtained from the WGM mode solver for the case of an ideal WGM whose cross-sectional refractive index profile is identical to that at ϕ₀.

It is worth emphasizing that the simple expression in (13) is valid for all whispering-gallery-mode microcavities, including those with or without axisymmetric properties, those that are subject to dielectric or metallic particle binding, or those in the presence of a tapered fiber or any other structures. The equation above can be further simplified if the perturbation or non-axisymmetric parts of the microcavity are localized to a small region, as is demonstrated in the application section below. Similarly, one may derive the overall quality factor of an asymmetric or perturbed WGM cavity from the algebraic average of the imaginary part of δm as

As can be seen, in the absence of the axisymmetry, an additional loss term from the mode mismatch between neighbouring cross sections occurs. This yields the decrease of overall Q_tot factor from the inclusion of the mode mismatch Q_m.

3. Application and discussion

3.1. The ideal WGM cavity

To validate the formulated algorithm, we first examine the case of a WGM with perfect axisymmetry. It’s evident that the WGM mode at each cross section is identical and that the mode mismatch loss vanishes. In addition, λ_r is identical at any azimuthal angle ϕ and the cavity resonance wavelength λ_res = λ_r. Therefore, an ideal WGM cavity can be treated as a special case under the current MMM formulation.

In the presence of a tapered waveguide that couples the light in and out of the cavity, mode mismatch occurs within the cavity-taper interaction regime. The coupling loss can be viewed as a special case of mode mismatch loss triggered by the waveguide taper. Assuming the taper interacts with the cavity within an azimuthal angle between ϕ₀ and ϕ₀ + Δϕ, we obtain a new formulation of coupling Q as

Surface roughness of the cavity often contributes to the loss as an individual term [11] in the context of MMM, wherein the loss also originates from the mode mismatch between cavity cross sections. Using the same approach, one can obtain the surface roughness associated Q using the equations obtained above, although an integral over a 2π azimuthal angle must be performed. Evidently, the surface roughness could perturb optical path length and yield an addition shift of cavity resonance wavelength. This shift, however, would usually be small in magnitude.

3.2. Nanoparticles on a WGM cavity

In this test case we examine the resonant wavelength shift and decrease of quality factor when N nanoparticles land on the surface of a WGM cavity. Here we assume each nanoparticle i occupies a space lying between ϕ_i and ϕ_i + Δϕ_i. From (13) and (14) we derive that the particle induced resonance wavelength shift Δλ and quality factor Q change follows the expression

Equations (17) and (18) are numerically verified at an operational wavelength around 633 nm on a silica microtoroid platform [2]. The silica microtoroid under investigation has a major radius of 40 μm and a minor radius of 5 μm. The surrounding environment is filled with water. Polystyrene (PS) beads and gold (Au) beads are individually placed on the toroid equator and the refractive indices of silica (1.457+(6.95 × 10⁻¹¹) j), water (1.33168+(1.47×10⁻⁸) j), PS (1.583+(5.29×10⁻⁴) j), and Au (0.1834+3.433 j) are taken from literature [35–38]. We employed a full vector finite element whispering-gallery-mode solver according to [20] to obtain the mode at each cross section. The computation domain was meshed adaptively in order to improve the accuracy of the calculated eigen value. The trapezoidal rule is implemented to perform the integral (17) and (18) numerically to get the total resonant wavelength shift and Q factor degradation of the hybrid cavity mode. Note that for the same mesh density, [25] requires approximately 90 times more grid points than MMM in each computation step assuming a minimum-size slice of a half-wavelength in length along the cavity edge was adopted in [25]. It is demonstrated that one typical MMM simulation of a hybrid WGM-NP system for a 20-nm radius gold bead takes 33 seconds for a relative error of the calculated resonance shift around 0.3 when only two slices along the azimuthal angle are involved. Higher accuracy can be achieved with more densely packed slices along the azimuthal direction and a more accurate mode solver. Note that the calculation is not based on the same mesh scheme as in [25] as our algorithm only requires a two-dimensional mesh at each computation step.

Figure 2 shows the intensity of the normalized toroid-bead hybrid modes at the azimuthal cross section, where the center of the beads are located. As shown in Fig. 2(a), where a polystyrene bead was adsorbed at the equator of the toroid, the intensity around the bead (c.f. Fig. 2(a) inset) displays little distortion from the unperturbed toroid mode. In this case, the resulting resonance shift of the cavity can be extracted with sufficient accuracy using first order perturbation theory [18]. Fig. 2(b) shows the WGM-SP hybrid mode excited at the gold-dielectric interface. As shown in Fig. 2(b)’s inset, there is strong field intensity stored in the surface plasmon resulting in two hot spots aligned along the same direction as that of the electric field (along the ẑ direction displayed for the TE mode). Due to the coupling of the WGM and the SP, the field intensity outside the cavity is stronger than that inside the cavity. Such significant distortion to the WGM precludes the treatment of such problems with first order perturbation theory. On the other hand, MMM intrinsically incorporates the change in mode field distribution and yields accurate evaluations of the perturbed system.

 

Fig. 2: The fundamental mode intensity distribution of a silica microtoroid with (a) a polystyrene bead and (b) a gold bead bound to the equator at a 633-nm wavelength. The modes are plotted at the azimuthal cross section where the center of the bead is located. The insets provide a zoomed-in view of the intensity distribution around the beads.

Download Full Size | PDF

Figure 3 illustrates the calculated terms m_r, m_m, and m_abs arising from the PS (Fig. 3(a)) and gold (Fig. 3(b)) nanoparticles. It is observed that the mode mismatch loss dominates the degradation of the Q factor over material absorption in both cases.

 

Fig. 3: The real and imaginary part of the mode order m along the propagation direction when (a) a 50-nm radius PS bead and (b) a 50-nm radius Au bead are placed at ϕ = 0.

Download Full Size | PDF

In Fig. 4 we further investigate the convergence rate as a function of azimuthal angle discretization in the case of a 50-nm radius gold bead binding to the toroid. The percent error was calculated by comparing the results of each δϕ to the expected result obtained through the Richardson extrapolation procedure. One can observe that the simulated resonance wavelength shift converges on the order of O(δϕ²) (black line in Fig. 4(a)) while the convergence rate of Q is O(δϕ^0.9)(black line in Fig. 4(b)). A higher convergence rate is possible by implementing higher order integration procedures to evaluate (17) and (18).

 

Fig. 4: Shift and Q factor vs. grid spacing δϕ along the ϕ̂ direction for a 50-nm radius gold bead. The last point is omitted for the creation of the line of best fit.

Download Full Size | PDF

In Fig. 5, we display the resonance wavelength shifts and the degraded Q factors due to the binding of a polystyrene (Fig. 5(a)) and gold (Fig. 5(b)) bead of different radii, where an unperturbed toroid has a resonant wavelength of 632.747 nm and a theoretical Q factor of 1.65 × 10⁹. The resonance shifts predicted by MMM (red solid square markers) are compared with those predicted by a first order perturbation method (black plus symbol markers). For PS beads, MMM predicts wavelength shifts which are in good agreement with those predicted by the first order perturbation method. The quality factor (blue solid circle markers) of the toroid is unaffected when a bead with radius smaller than 5 nm is attached to the surface, yet there is a drop to approximately 10⁷ if a 50-nm radius bead is attached instead. The calculated wavelength shift and quality factor degradation is in line with the experimental observation reported in [13].

 

Fig. 5: Binding shift and Q factor degradation due to (a) a bound PS sphere and (b) a bound Au sphere for different bead radii.

Download Full Size | PDF

For the case of a gold bead simulation, the excess polarizability is given by α = 4πR³ε_m(ε_p − ε_m)/(ε_p + 2ε_m) [39] where ε_p and ε_m are the permittivities of the particle and the medium, respectively. MMM predicts resonance shifts from 2×10⁻² to 359 fm and Q factor degradation from 3 × 10⁸ to 7.5 × 10⁴ when the attached gold bead radius increases from 2.5 nm to 50 nm. As expected, at an off-plasmon resonance wavelength of 633 nm, the perturbation method matches the MMM results when the bead size is below 20 nm yet for large beads yields greater errors due to the non-negligible field distortion from the bead.

In Fig. 6 we also compared the shifts calculated by the two methods as a function of wavelengths. As illustrated, near the surface plasmon resonance of a gold nanobead (540 nm, as shown in the insert), the accuracy of the perturbation method drops below 50% due to the strong surface plasmon effect that was similarly observed in [19].

 

Fig. 6: Cavity resonance shifts as a function of wavelengths for a 25-nm radius gold bead. The insert shows the excess polarizability of the bead at different wavelengths.

Download Full Size | PDF

To verify our formulation we further compared our results with reported experiments. We simulated the binding of 27.5-nm gold beads to a 50-μm radius microsphere with a 4 nm inter-particle gap as reported in [16, 17]. In this simulation, we obtained a resonance wavelength shift of 32.6 fm per bead adsorption. Compared to the total resonance shift in Fig. 1(d) of [16], we predict that a total of 10 beads are adsorbed to the microsphere surface. Knowing that an unperturbed cavity quality factor of 3 × 10⁶ was reported in the same paper, we estimated that the quality factor will drop to 5.8 × 10⁵ by the adsorption of 10 beads. This value is in good agreement with the reported experimental value on the order of 10⁵.

It is worth noting that the MMM developed in this paper assumes that the cavity operates in a single whispering-gallery-mode. In practice, a WGM may contain multiple propagating modes and so omitting the coupling between different modes could decrease the accuracy of the model. A MMM formulation involving multiple WGM modes is possible, albeit at the expense of higher computational resource requirements.

4. Conclusion

A FEM based mode-matching method (FEM-MMM) has been developed for solving three-dimensional whispering-gallery-mode cavity problems. To demonstrate the validity of the method, numerical results for a microtoroid bound with a plasmonic and non-plasmonic bead were presented. Results obtained in this paper are consistent with those obtained through perturbation theory or experimental demonstration. The investigation of MMM with the inclusion of higher order modes will be carried out in forthcoming publications.