COMSOL Simulations of the Main Properties of the Whispering Gallery Mode Resonators for Bio-Sensing Applications

1. Introduction

Whispering gallery mode resonators (WGMRs) are curved shape optical structures, like spheres [1], microdisks [2], microrings [3], and toroids [4] that confine light due to the total internal reflection [5]. WGMRs are interesting objects because they can form optical resonances with ultra-high quality (Q) factors of ${10}^9-{10}^{11}$ [6,7]. Additionally, they concentrate light to such an extent that nonlinear optical phenomena appear. Microresonator frequency comb generation is one of the derived effects [8]. Applications of WGMRs include frequency filtering, lasing, temperature, force, and displacement sensing, electric- and magnetic-field sensing, and bio-sensing [9,10,11,12].

WGMRs are tested experimentally mainly by monitoring light intensity that is transmitted through them via prism or tapered fiber coupling [13,14]. Linewidth scan of laser emission accompanied by time-domain plot of transmitted signal reveal the shape and the behavior of resonance signals of the WGMRs [15,16]. The full width at half maximum (FWHM) of the resonance peaks and their frequency positions are parameters that describe the interrelation between pumping laser light and the surrounding environment through these resonators. The shift in the observable resonance data reflects the changes of the resonator or the external environment, thus becoming a tool of sensing.

Coating the WGMRs with a functionalizing layer to enhance optical proprieties impacts both the Q factor and the resonance peak location describes the ability to store energy inside a microcavity [16,17].

COMSOL Multiphysics with Wave Optics Module is a popular finite element software that is used to model light propagation in WGMRs [18,19,20]. COMSOL provides neat visualization of light field in optical structures as color-map values of its electric and magnetic components and energy density. Numerical simulations allow easier to understand effects of light propagation in WGMRs, provide additional details on the field structure in addition to resonance positions of WGMRs and are easier to comprehend than advanced analytical theory of WGMRs based on Maxwell’s equations [21].

COMSOL simulations provided in this paper show the basic and advanced effects of WGMRs: light can be coupled in them by prism and tapered fiber structures, resonance positions and Q factors of resonances depend on the coupling parameters, "biomolecules" attached to the resonator shift resonance positions, surface functionalization, and surface roughness affect both resonance position and Q factor, and others. Basic theory is provided to be comparable to simulation data.

2. Basic Theory of WGMRs

2.1. Resonances and free spectral range

When the length L of the light path loop in an optical structure with the refraction index n is equal to the positive natural number m of wavelength in this media $(\lambda /n)$, a resonance condition is formed:

$$\begin{array}{c}\hfill L=\left(\frac{\lambda }{n}\right)m,\end{array}$$

(1)

which gives resonance wavelengths:

$$\begin{array}{c}\hfill {\lambda }_m=nL·\frac{1}{m}.\end{array}$$

(2)

Corresponding resonance frequencies are:

$$\begin{array}{c}\hfill {\nu }_m=\frac{c}{nL}·m.\end{array}$$

(3)

If light path loop L is a circle with radius a then

$$\begin{array}{c}\hfill {\lambda }_m=2\pi an·\frac{1}{m},\end{array}$$

(4)

$$\begin{array}{c}\hfill {\nu }_m=\frac{c}{2\pi an}·m.\end{array}$$

(5)

If resonances are searched in a circular resonator with radius r, then (2) and (3) can be used with a condition $a=r$ only as a rough approximation. Advanced theory based on Maxwell’s equations [22,23] show that for a first radial mode of spherical WGMRs placed in air resonances are observed for corresponding TE and TM modes:

$$\begin{array}{cc}\hfill {\lambda }_m\left(TE\right)=& 2\pi rn·\frac{1}{m_1},\end{array}$$

(6)

$$\begin{array}{cc}\hfill m_1\approx & m\left(1+\frac{1.856}{m^{2/3}}+\frac{1}{m}\left(\frac{1}{2}-\frac{n}{\sqrt{n^2-1}}\right)\right),\end{array}$$

(7)

$$\begin{array}{cc}\hfill {\lambda }_m\left(TM\right)=& 2\pi rn·\frac{1}{m_2},\end{array}$$

(8)

$$\begin{array}{cc}\hfill m_2\approx & m\left(1+\frac{1.856}{m^{2/3}}+\frac{1}{m}\left(\frac{1}{2}-\frac{1}{n\sqrt{n^2-1}}\right)\right).\end{array}$$

(9)

Similarity between equations (4), (6), (8) is seen. Let us explore the situation when a resonator with $n=1.4537$ (fused silica, 780 nm region) is used. Then m = 5 gives $m_1\approx 7.3$ and $m_2\approx 8.0$, and m = 500 gives $m_1\approx 513.9$ and $m_2\approx 514.8$. This shows that if $a=r$ is used in (4) then the number m of wavelengths in some resonance differs from actual numbers $m_1,m_2$ significantly. Resonances for fused silica WGMR sphere can be determined near 780 nm spectral position. According to equation (6) - (9) the parameters $\lambda$, $\nu$ and m for resonators with radius r = 2 µm, r = 20 µm and r = 200 µm are given in Table 1. These values give orientation for the search for resonances in COMSOL models.

The distance between two closest resonances is called the Free Spectral Range (FSR). For wavelength scale, FSR is:

$$\begin{array}{c}\hfill FSR\left({\lambda }_m\right)={\lambda }_{m+1}-{\lambda }_m=-\frac{{\lambda }_m^2}{(nL+{\lambda }_m)}\approx -\frac{{\lambda }_m^2}{nL},\end{array}$$

(10)

which means that peaks are not equidistant. We assumed that $nL>>{\lambda }_m$, which is valid for large resonators.

In frequency scale, FSR is:

$$\begin{array}{c}\hfill FSR\left({\nu }_m\right)={\nu }_{m+1}-{\nu }_m=\frac{c}{nL},\end{array}$$

(11)

which means that resonances are equidistant.

2.2. Quality factor, finesse and intensity distribution of resonances

The quality factor of the resonator system is defined as the ability to store the energy inside the cavity [24,25]:

$$\begin{array}{c}\hfill Q=2\pi \frac{storedenergy}{energylossperoscillationperiod}.\end{array}$$

(12)

If amplitude of the optical signal oscillates and decays as $e^{i\omega t}{e}^{-t/\left(2\tau \right)}$, then

$$\begin{array}{c}\hfill Q=\tau \omega =2\pi \nu \tau =2\pi c\tau /\lambda .\end{array}$$

(13)

The full width at half maximum (FWHM) of resonances for $\omega$, $\nu$, and $\lambda$ are:

$$\begin{array}{ccc}\hfill \Delta \omega & =& \frac{1}{\tau },\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill \Delta \nu & =& \frac{1}{2\pi \tau },\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill \Delta \lambda & =& \frac{{\lambda }^2}{2\pi c\tau }.\hfill \end{array}$$

(16)

And Q factor can be expressed as:

$$\begin{array}{c}\hfill Q=\frac{\omega }{\Delta \omega }=\frac{\nu }{\Delta \nu }=\frac{\lambda }{\Delta \lambda }.\end{array}$$

(17)

In case of various processes with corresponding Q-factors $Q_i$ that influence the system, the final Q-factor of the system is calculated as:

$$\begin{array}{c}\hfill \frac{1}{Q}=\frac{1}{Q_{rad}}+\frac{1}{Q_{ss}}+\frac{1}{Q_{cont}}+\frac{1}{Q_{mat}}+\cdots .\end{array}$$

(18)

For typical WGMRs these sub-processes are [6]: $Q_{rad}$ is formed by intrinsic radiative (curvature) losses, $Q_{ss}$ – scattering losses on residual surface inhomogeneities, $Q_{cont}$ – losses introduced by surface contaminants, $Q_{mat}$ – material losses.

Finesse describes the resonator and is defined as FSR divided by FWHM:

$$\begin{array}{c}\hfill \mathcal{F}=\frac{FSR}{\left(FWHM\right)}.\end{array}$$

(19)

Taking into account formulas (11), (13), (15), and (19), we get:

$$\begin{array}{c}\hfill \mathcal{F}=\frac{c}{nL\Delta \nu }=\frac{2\pi \tau c}{nL}=\frac{Qc}{nL\nu }.\end{array}$$

(20)

Here we will describe a model of a waveguide coupled to a circular resonator taking into consideration the field decay in the system. This system is called an optical all-pass filter. The coupling of waveguide and resonator is described by reflection coefficient u, transmission coefficient $t=\sqrt{1-u^2}$, and one loop path length in the resonator L. This gives phase shift per loop $\varphi =2\pi nL/\lambda$, and the field intensity decay rate per loop is $e^{-\beta }$ with $\beta =t_0/\left(2\tau \right)=nL/\left(2c\tau \right)$, where $\tau$ is the decay period, $t_0$ is the time needed for the light to travel one loop in the resonator, and c is the speed of light. Incident light intensity is $I_0$.

Transmitted field intensity of all-pass filter is [26]:

$$\begin{array}{ccc}\hfill I& =& I_0-\frac{I_{1max}}{1+{(2\mathcal{F}/\pi )}^2{sin}^2(\varphi /2)},\hfill \end{array}$$

(21)

where

$$\begin{array}{ccc}\hfill I_{1max}& =& \frac{(1-u^2)(1-e^{-2\beta })}{{(1-u{e}^{-\beta })}^2},\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill \mathcal{F}& =& \frac{\pi \sqrt{u}{e}^{-\beta /2}}{1-u{e}^{-\beta }},\hfill \end{array}$$

(23)

Using (17) and (20) we get analytical values for:

$$\begin{array}{ccc}\hfill \Delta \nu & =& \frac{c}{\pi nL}\frac{1-u{e}^{-\beta }}{\sqrt{u}{e}^{-\beta /2}},\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill Q& =& \frac{L}{(\lambda /n)}\frac{\pi \sqrt{u}{e}^{-\beta /2}}{1-u{e}^{-\beta }},\hfill \end{array}$$

(25)

Minimal and maximal value of a transmitted field (21) is:

$$\begin{array}{ccc}\hfill I_{min}& =& I_0{\left(\frac{u-e^{-\beta }}{1-u{e}^{-\beta }}\right)}^2,\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill I_{max}& =& I_0{\left(\frac{u+e^{-\beta }}{1+u{e}^{-\beta }}\right)}^2,\hfill \end{array}$$

(27)

giving a resonance depth:

$$\begin{array}{cc}\hfill I_{res}& =I_{max}-I_{min}=I_0\frac{4u{e}^{-\beta }(1-u^2)(1-e^{-2\beta })}{{(1-u^2{e}^{-2\beta })}^2}.\end{array}$$

(28)

Close to resonances ${sin}^2(\varphi /2)\approx {((\varphi -{\varphi }_m)/2)}^2$ and intensity (21) can be approximated as:

$$\begin{array}{c}\hfill I=I_0-\frac{I_{1max}}{1+K{(\nu -{\nu }_m)}^2},\end{array}$$

(29)

with K as a coefficient. Absorption resonance according to (29) is a Lorentzian function.

Equations (21 - 29) derive the main characteristics of resonances of WGMRs - resonance depth $I_{res}$, resonance width $\Delta \nu$, and resonance shape - that can be visualized using mathematical simulations.

3. Modelling of WGMRs

COMSOL Multiphysics is finite element numerical calculation software that allows one to model various physical processes like light propagation in optical systems, heat propagation, mechanical stress distribution, etc. Wave Optics module has to be activated to describe light in optical systems as waves. WGMRs can be modeled in 2D, 2D axisymmetric, or 3D mode. Resonance frequencies are obtained in "Eigenfrequencies" study selection or by scanning the frequency of the incident field and monitoring the outcome port or value of the Energy density time average of the resonator.

Finite element numerical calculation software derives values of parameters in various spatial structures and computing time for this depends on the number of mesh points. For exploration of wave optics processes mesh size has to be smaller than the fraction of wavelength of the field. This limits the application of these calculations to 2D and 3D structures that are much larger than the wavelength of the field. Large structures can be modeled relatively fast only in 2D axisymmetric modes. Therefore 3D mode was used to model WGMRs radius of up to 2 µm, 2D mode was used to model WGMRs radius of up to 20 µm, and 2D axisymmetric mode was used to model WGMRs radius up to 200 µm. When circle resonators are explored in 2D mode they represent cylindrical structures. Spherical resonators can be modeled only in 2D axisymmetric and 3D modes.

Refraction coefficients for the resonator, prism, and fibers were selected to be equal (n=1.4537, fused silica) and resonances were explored in the spectral region around 780 nm.

3.1. Free space, prism and tapered fiber coupling of cylindrical resonator

When a wide plane wave is irradiated on a cylindrical resonator (2D mode) then no resonances are observed and we find that the resonator works only as a focusing lens. When we irradiate spatially limited plane wave on a cylindrical resonator and it touches only tangentially its perimeter then WGM resonances are formed (Figure 1a). Plane wave moves from left to right. Colors in the Figure represent values of $Ez$ component of the intensity of the electric field of a light wave. This is known as free space coupling, which is an inefficient way of coupling light. The efficiency is less than ${(\lambda /a)}^2$ [27], where a is the beam cross-section radius (comparable with sphere radius). This approach is used experimentally to couple light inside a liquid spherical resonators, for example, glycerol droplet [28].

Prism-type structure can be used to introduce light in the resonator (Figure 1b). When scanning incident wavelength and measuring energy density time average integrated across resonator, we get resonance spectra in the resonator. Experimental studies have shown that prism coupling has approximately 80% efficiency [29,30]. This coupling method is robust. Additionally, the surface of the prism can be cleaned, which is advantageous for different bio-sensing applications in liquids because the molecules may attach to both the surface of the WGMR and the prism.

Experimentally it is easier to couple light inside the WGMR using a tapered optical fiber – a single-mode optical fiber that is significantly thinned with a waist less than 10 µm in diameter (Figure 1c). It has 98-99% coupling efficiency [30]. Optical transmittance of tapered fiber degrades due to particle adhesion on the surface. However, cleaning tapered fiber is more challenging than a prism because it is very fragile. Therefore, this coupling method has severe limitations in bio-sensing applications.

Besides the coupling method, various coupling conditions can be modeled. When changing the gap between the resonator and prism as gap $=g·{\lambda }_0$ with ${\lambda }_0$= 780 nm and g changing in the region 0 – 0.5 we get that resonance peak position almost does not change (Figure 2a), but the intensity of resonance peak at first increases and then decreases with largest intensity to be achieved when parameter g ≈ 0.1 (Figure 2b). The simulation results conform with the experimental results when varying the gap between prism and WGMR [31]. An optimal gap allows critical coupling when intrinsic and extrinsic losses match.

Resonance intensities also change when the incident angle $\alpha$ of the light field changes (Figure 2c). Maximum intensity in this case is dependent on the gap between the resonator and prism and is achieved when the incident angle of the light is slightly larger than the angle of total internal reflection $\gamma =arcsin(1/n)$ which is 0.76 rad for fused silica. In case of a larger gap, the coupling appears to be achieved just before reaching the boundary $\gamma$, which could be explored further. Otherwise, maximums in Figure 2c below $\gamma$ do not belong to WGM resonances.

For tapered fiber coupling, the fiber width (${\lambda }_0/n·g$) was changed to monitor intensities of resonance peaks (Figure 2d). It was found that the largest peak of energy density time average in the resonator appears when the width of the fiber is below the wavelength. This explains the need for tapered fibers (< 2µm in diameter) for effective coupling of light into the WGMRs [32,33].

3.2. Resonances and bio-sensing using WGMRs

WGMRs can be used as bio-sensors [10,34]. The sensing principle explored in this paper is the shift of the resonances due to changes on the surface of the resonator.

The process of bio-molecules attaching to the surface can be modeled by small elements that attach to the surface of the resonator. Figure 3a–c gives modeling examples of this process when multiple small "items" touch the surface. In Figure 3a five small items (R = 0.1 µm) are attached to the surface of the resonator (R = 2.0 µm). WGMR is coupled using the coupling prism. In Figure 3b four particles (R = 0.5 µm) are attached to the surface of the resonator (R = 6.0 µm). WGMR is coupled using the tapered fiber. In Figure 3c four different irregularly shaped particles are attached to the surface of the resonator. In all three cases, particles act as surface defects and scatter light. For small particles, when irradiation wavelength $\lambda$ is comparable with the particle size R scattering dominates over absorption [35]. For large particles, where $R\gg \lambda$, the absorption dominates over scattering [35].

Figure 3d–f shows the shift of WGM resonance and degradation of the Q factor due to attachment of spherical (as shown in Figure 3b) and irregularly shaped (as shown in Figure 3c) particles. For Figure 3d, round dust like particles are modelled using R = 0.5 µm, $n_d$ = 1.5 and $k_{id}$ = -0.0007, and, for Figure 3f round gold particles are modelled using R = 0.5 µm, $n_{Au}$ = 0.14891, $k_{iAu}$ = 4.783. WGMR size in all three cases is R = 6.0 µm. As each addition particle attaches to the surface (in simulations particles are attached from right to left as illustrated in Figure 3b-c), the electric field norm of the WGM resonance decreases. The shift of the resonance wavelength depends on the parameters n and k. For round dust particles, the resonance wavelength increases with the number of particles attached to the surface (see Figure 3d inset). Particle shape and size have some impact as well. In the case of four round identical dust particles WGM resonance position shifted 0.72 nm. In the case of four irregularly shaped and sized particles WGM position shifted 0.63 nm. For gold particles, the resonance wavelength decreases with the number of particles attached to the surface. WGM resonance shifted -0.50 nm. In all cases, the Q factor of the WGMR decreases with attached particles, which acts as a surface defect and increases surface scattering losses.

For bio-sensing applications, WGMRs need to be functionalized with a bio-sensitive layer by introducing a binding material to the surface that attracts only specific molecules. The surface functionalization layer can be modeled by adding a thin layer of material with different optical parameters [36,37]. We demonstrate examples with Polymethyl methacrylate (PMMA) ($n_{PMMA}$=1.4860, k=1.4000$·{10}^{-7}$ [38]) and ZnO ($n_{ZnO}$=1.7588, k=0.030390 [39]) layers, which influence the degree of WGM confinement inside the resonator. PMMA layer has a lower refraction index than silica and it makes the WGMR more sensitive to the surrounding environment. ZnO layer has a higher refraction index than silica and it ensures the stability of the whispering gallery mode.

A more realistic rough surface can be generated by using built-in random functions generating a random polar curve [40]:

$$\begin{array}{ccc}\hfill x& =& cos\left(2\pi S\right)·\left(R_0+dr·\left|\sum _{\begin{array}{c}m=-F\\ m\ne 0\end{array}}^F{\left(m^2\right)}^{-\frac{b}{2}}·g_1\left(m\right)·cos\left(2\pi mS+u_1\left(m\right)\right)\right|\right)\hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill y& =& sin\left(2\pi S\right)·\left(R_0+dr·\left|\sum _{\begin{array}{c}m=-F\\ m\ne 0\end{array}}^F{\left(m^2\right)}^{-\frac{b}{2}}·g_1\left(m\right)·cos\left(2\pi mS+u_1\left(m\right)\right)\right|\right)\hfill \end{array}$$

(31)

where $S\in [0;1]$, $R_0$ is the radius of the random polar curve, $dr$ (roughness parameter) is an amplitude multiplier for the random part of the equation, b is a spectral exponent, $g_1$ is a normal distribution random function, $u_1$ is a uniform distribution random function, and F is a spatial frequency resolution. We take the absolute value of the random function so the curve wouldn’t cross the smooth resonator line and would only go outwards of the resonator as a functionalized surface would. The resulting jagged line better approximates a real resonator with an imperfect surface. Figure 4a–c shows the distribution of electric field norm in fused silica spherical WGMR with radius 5 µm coated with 100 nm PMMA and roughness parameter $dr$ 3.7 nm, 9.1 nm and 17.3 nm respectively. A rough surface contributes to the surface scattering losses. The shapes of the electric field norm of the resonance become less regular. Additionally, the rough surface affects the effective refraction index of the WGMR and WGM resonances shift as can be seen in Figure 4d,e for both fused silica WGMR and fused silica WGMR coated with 100 nm PMMA layer. Initially, WGM resonance is shifted by 0.91 nm with the PMMA layer (see Figure 4f). When increasing the roughness the WGM resonance shifts more for WGMR with the PMMA layer than for pure silica WGMR. For 20 nm roughness, the shift difference increased to 1.96 nm. As the roughness increases, the WGMR Q factor gradually decreases for resonators both with and without PMMA layer (see Figure 4g).

For a ZnO layer, the surface roughness plays a more important role as shown in Figure 4h. As the WGM resonance moves into the coating [17] it is more sensitive to surface roughness and affected by surface scattering losses.

COMSOL 2D axisymmetric mode can be used to find resonances of WGMRs, it is suitable for large resonators. The advantage of 2D axisymmetric mode is that we can investigate various side modes the effective mode area and the effective mode radius. It is also used to calculate the dispersion of the resonator. However 2D axisymmetric mode can’t be used to investigate the coupling efficiency with a light source and effects of attached bio-molecules or rough surface.

4. Conclusions

The description of the main effects of WGMRs was given with corresponding COMSOL simulations. Although plane wave irradiating on a resonator does not create resonances, it was shown that they can be obtained when focused tangentially on the perimeter of the resonator. It was shown that resonances can be excited in resonators by prism coupling when the incident angle of the light is slightly larger than the angle of total internal reflection. There was the optimal distance of the fraction of wavelength between the prism and resonator to get the strongest resonances. The width of fiber coupled to the resonator has to be comparable to the wavelength of the light to achieve the strongest resonances thus providing theoretical ground for fiber tapering. Bio-sensing principle of WGMRs was shown with the demonstration that the resonance shift increases with a larger number of "bio-molecules" and particles attached to the surface of WGMR. The roughness surface functionalization layer was more important for layers with a higher index of refraction. The provided results form an understanding of the main effects of WGMRs and surface functionalization. In the future modeling results could be used to explain and accompany experimental results of WGMR surface functionalization or molecule attachment. Some aspects of coupling parameters and surface roughness also could benefit from a more detailed modeling and a more thorough results analysis.