Due to ultra high quality factor (106 − 109), axisymmetric optical microcavities are popular platforms for biosensing applications. It has been recently demonstrated that a microcavity biosensor can track a biodetection event as a function of its quality factor by using phase shift cavity ring down spectroscopy (PS-CRDS). However, to achieve maximum sensitivity, it is necessary to optimize the microcavity parameters for a given sensing application. Here, we introduce an improved finite element model which allows us to determine the optimized geometry for the PS-CRDS sensor. The improved model not only provides fast and accurate determination of quality factors but also determines the tunneling distance of axisymmetric resonators. The improved model is validated numerically, analytically, and experimentally.
© 2013 Optical Society of America
Microcavities are widely used in many applications such as biosensing, nonlinear processes, lasers and optomechanics. Most previously published work with microcavity biosensors has employed measurement of the change in resonant frequency of the sensor as a function of a biodetection event [1, 2] with only a few instances of using the quality factor as a sensing metric [3–5]. It has been shown that the quality factor of microcavities can be measured in the time domain using phase shift shift cavity ring down spectroscopy (PS-CRDS) . This approach has recently been applied to biosensing applications . In order to extend this work to ultra sensitive biosensing applications (e.g. determination of binding kinetics, detection of pM − fM concentrations, and proteins, viruses etc), it is important to design a cavity with the optimum parameters for high sensitivity. One of the goals of the current work, is to rapidly estimate the quality factors of the microcavities with high accuracy to address these applications.
Whispering gallery mode (WGM) radiation occurs due to tunneling of photons through a potential barrier from the microcavity to a allowed wave propagation zone (see Fig. 1). The tunneling phenomenon has been utilized in various photonic devices such as photonic band gaps , super lattices , and asymmetric cavities . Recently Tomes et. al.  have performed remarkable experiments to image the tunneling process of axisymmetric microcavities. These experiments indicate that the tunneling phenomenon can be utilized for potential applications such as biosensing. However, to design the appropriate devices for a particular biosensing application, it is again necessary to accurately calculate the tunneling distance of axisymmetric microcavities. Another goal of our work is to address this issue.
In any modeling technique, reflections from boundaries of the computation domain are induced due to radiation produced by the whispering gallery modes (WGM) of a microcavity. For an accurate electromagnetic model, absorbing boundary conditions or perfectly matched layers (PML) are required to reduce these unwanted reflections. PML act as artificial boundaries that truncate the computation domain of open region scattering problems in the finite element method. There have been previous attempts to develop finite element models (FEM) of axisymmetric cavities which incorporate the PML. In 2005, Chinellato et al.  introduced a FEM model which was implemented in MATLAB. Whilst their model was capable of simulating the behavior of very small resonators (< 3μm) it was not adequate for resonators possessing the dimensions typically used in experiments. In 2009, Karl et al.  developed a 3D FEM model in JCMsuite for studying a micro-pillar cavity. However, this model does not consider the suppression of false solutions, a well known problem in finite element formulations [14, 15]. Moreover, the quality factor of the modes was estimated by fitting the Lorentzian peak to the calculated spectrum of the cavity, thus introducing an extra approximation.
In 2007, Oxborrow  developed a FEM for open axisymmetric resonators in COMSOL without invoking any transverse mode approximation to Maxwell’s equations, representing an advance on previous work [17, 18]. In his work, he showed that the model could simulate resonators of arbitrary cross section in optical and microwave regimes, thus removing the size limitation in . Another difference from  is in terms of suppression of false modes; Oxborrow used a simple penalty term in his master equation whereas Chinellato et al. used Nédélec edge and modified Lagrange nodal element functions to avoid the spurious modes. With Oxborrow’s formulation in COMSOL, 3D rotationally symmetric problems are reduced to 2D and are solvable in seconds which is vast improvement over [12,13]. (For further discussion and further comparison to other works see references [16, 19].)
The advantages of Oxborrow’s model have made it a popular choice among researchers and it has been widely used in numerous research works e.g. [20–23]. However in Oxborrow’s model no PML was implemented and as a result the WGM quality factor could not be determined accurately. The quality factor due to the WGM radiation was estimated by placing a bound on its minimum and maximum possible values. These maximum and minimum values were determined by executing the model multiple times with different boundary conditions. Moreover, the model lacks the capability to estimate the quality factor of multiple modes simultaneously.
In order to provide accurate determination of the WGM quality factor, we have improved Oxborrow’s model by modifying its master equation and implementing the PML along the boundaries of the computation domain. In the present work, geared towards sensing applications, we expand and refine the model presented in our earlier work . Other researchers  have subsequently reported the inclusion of PML in Oxborrow’s model for determining resonant frequencies and corresponding mode profiles of microring resonators but mathematical details have not been provided.
Our modified model does not have any of the drawbacks of Oxborrow’s model. Moreover, we have computed the quality factors of all the modes without using any fitting algorithms as opposed to the approach taken in . Furthermore, with the modified Oxborrow’s method, tunneling distances can also be accurately extracted for microcavities of various shapes. In the present work, a simple expression for computing tunneling distance of microtoroidal cavities is also provided.
In our model, we treat the PML as an anisotropic absorber and implement it in the cylindrical coordinate system. Our model is applicable to any axisymmetric resonator geometry but due to the availability of analytical expressions for spherical resonators, we have validated the model by determining the quality factors and tunneling distances of a silica microsphere in air. We have found that our simulation results are in excellent agreement with the analytical results. We also apply our model to microtoroidal cavities immersed in liquid and show that our results are consistent with those obtained by experiments. The model is then used to determine the optimum parameters of a microtoroidal cavity sensor based upon phase shift cavity ring down spectroscopy.
The paper is organized as follows. In section 2, an improved FEM model is introduced by incorporating the PML. The optimal parameters of the PML are discussed in section 2.1. The analytical expressions of the quality factor and the tunneling distance of microspheres are presented in section 2.3. These expressions are then used in section 3.1 for comparing the modeling results. In the same section, an empirical expression for tunneling distance of fundamental modes of a toroidal cavity is also provided. In section 3.2, the model is validated experimentally by measuring the quality factors of the microtoroidal cavities in liquid. In section 3.3, the model is validated numerically by running convergence test. In section 4, we show that the model can be applied for sensing applications. In this section, we show that for a sensor, whose sensing metric is change in the quality factor, an optimum geometry exists for achieving maximum sensitivity.
2. Mathematical description
Applying Galerkin’s method to the wave equation and after using the boundary conditions for open resonators, one can arrive at the FEM equation in the weak form :Eq. (1) represents a penalty term to suppress false solutions. None of the field components will depend upon the azimuthal coordinate ϕ in the axisymmetric resonators, resulting in reduction of the 3D problem to a 2D problem.
2.1. Perfectly matched layer formulationEq. (3)
where trpml, tupml, tlpml are the PML thicknesses in the radial, +z and −z directions respectively and rpml, zupml, zlpml are the locations of the start of PML in the radial, +z and −z directions respectively. nmedium is refractive index of the medium, N is order of the PML, and G is a positive integer.
In the PML expressions (sr, sz, r̃), the imaginary component contributes to the attenuation of waves in the PML but at the same time, due to the discrete nature of the FEM mesh, a large imaginary component will introduce reflections at the interface between the PML and the medium. In order to determine the optimal value for the imaginary component, we have investigated linear, quadratic, and cubic PML of different thicknesses for various values of G by running many simulations for various sphere diameters. To deduce the optimum values of the parameters, we then compared the simulation results for QWGM of spherical cavities with the analytical ones. The simulation results show that a linear (i.e. N = 1), and λ/4 thick PML with a G value of 5 is optimum. We have also used these optimum values for the simulations of the microtoroidal cavities. The location of PML is also important to obtain accurate results. The radial PML should be greater than the tunneling distance (t) of the microcavity and the z PML should be greater than FWHM (wz) of WGM along the z axis. After running series of simulations we find the following as optimum values (w.r.t. (r,z) = (0,0), see Fig. 1):
2.2. Finite element method equation with perfectly matched layer
In order to incorporate the PML, we have reformulated Eq. (1) in the following way:Eq. (9) into the FEM software COMSOL, a full vectorial finite element model of a silica sphere in air can be obtained. By using the eigenvalue solver in COMSOL, resonant frequencies (fr) of all the modes can easily be determined. Quality factor due to the WGM radiation can be calculated as :
2.3. Analytical expressions of the spherical resonator
2.3.1. Quality factor
The quality factor due to WGM radiation losses in a spherical microcavity can be written as :
m can be calculated using the characteristic equation for WGM frequencies :
It should be noted that Eq. (11) is an asymptotic solution for the Qwgm of the spherical resonator which requires m ≫ 1, however the error is less than 1% for m ≥ 19 . In our comparison for the analytical results, we have applied Eq. (11) to silica spheres with mode numbers (m) ranging from 25 – 80.
2.3.2. Tunneling distance
The Schrodinger equation for a microsphere for a fundamental mode can be written in radial coordinates (r) as :30]:
3.1. Comparison between simulations and analytical results
Figure 1 shows a fundamental TE mode of a silica spherical cavity in air.
We have plotted the quality factor due to the fundamental TE whispering gallery mode radiation for various sphere diameters at 1550nm. Similar results are obtained for the TM mode. Figure 2 shows the comparison of the FEM simulation results and results obtained by using analytical expressions presented in section 2.3. We have also calculated the minimum and maximum Qwgm values using Oxborrow’s model  for each sphere diameter and those are also shown in Fig. 2. It can be seen that the new model provides a more accurate estimate than Oxborrow’s model, and returns a single value of Q rather than a range of values. The slight difference between the analytical and the FEM values is attributed to discretization of the computation domain. The accuracy will improve with finer mesh elements (see section 3.3). It can also be seen that Oxborrow’s bounds do not always straddle the analytical solutions. One of the possible reasons for this discrepancy may be that while deriving these bounds, in order to simplify the derivation, Oxborrow assumes that the modes are transverse but in reality the modes of these axisymmetric resonators are not perfectly transverse .
Figure 3 shows the results for tunneling distance for both microsphere and microtoroidal cavities. An empirical relation for the tunneling distance for fundamental modes of the micro-toroidal cavities can be extracted from the simulation results and is given as:Fig. 3(b). The tunneling distance results are also in agreement with experimental results presented in .
3.2. Comparison between simulations and experimental results
Our model is also appropriate for other axisymmetric resonators where analytical solutions are not available. In order to test this, we have used our FEM model to determine the total quality factor of microtoroidal cavities  and have compared the simulation results with experimental measurement of the quality factor. Mathematically the total quality factor (Qtotal) can be represented by:
In order to validate the model with reasonable range of quality factors, it is necessary to ensure suitable experimental conditions in which to observe these quality factors. We have performed the experiments with the microtoroids immersed in ethanol (refractive index: 1.3538 at 1550nm) by using the fluidic cell described in . We have preferred ethanol over water as ethanol not only provides low refractive index contrast between the silica cavity and the surroundings but also has low absorption as compared to water at 1550nm. This will allow us to select a wide range of microtoroidal cavities whose Qtotal is mainly limited by the Qwgm.
A broadband source (peak wavelength:1530nm) is used to couple the light into microtoroids via a tapered optical fiber (taper waist: ≤ 1μm). The resonant peaks are observed on a high resolution optical spectrum analyzer (Apex AP 2443B, resolution: 0.16pm) and quality factor is determined by Lorentz curve fitting of the peaks (Q = λ/Δλ). The tapered fiber is positioned (using 10nm resolution nanostages) along the equator of the cavity to couple light into fundamental transverse modes. It is ensured that the tapered fiber does not touch the cavity during the measurements. In order to minimize the effect of Qcoupling on Qtotal, the power at the peak wavelength of the source is set around −48dbm (Optical spectrum analyzer sensitivity: −70dbm) and all the measurements are taken in highly undercoupled regime. The comparison between simulation and experimental results is shown in Fig. 4.
3.3. Computational speed and numerical accuracy
For the results presented in Figs.1–4, we have used the quadratic Lagrange elements of triangular shape with average element size (longest edge) of 0.3μm and the number of degrees of freedom of order 105. As an example, a microtoroidal cavity model (Computation domain size:25μm × 16μm, D = 30μm, d = 5μm) with the aforementioned statistics took only 35s to find the first 20 eigenvalues on a quad core 64 bit operating system PC.
To check numerical accuracy of the model, we have run a convergence test. Figure 5 shows the plot of relative error (Er = (QFEM − Qexact)/Qexact) as a function of the degrees of freedom (DOF) of our model. The DOF is related to the number of mesh elements, and the basis functions used in the FEM solver. Convergence of the solution is of order 2 which is expected for the Galerkin method for quadratic Lagrange elements . It can be seen in the Fig. 5 that for small DOF, convergence rate is slower than the expected; the reason for this is the poor approximation of the curved boundary of the cavity by the coarsely meshed triangular elements.
4. Discussion and application to PS-CRDS microcavity sensor
The results presented in section 3 show that our finite element model is both physically and numerically accurate. The quality factors for all the modes (fundamental and higher order modes) are also obtained by one single simulation rather than multiple simulations, which was the case for the original Oxborrow model. Moreover, no prior knowledge of any of the mode frequencies is required to obtain the quality factors. The model also gives fast results and without any fitting algorithms.
It should be noted that we have neglected dispersion for estimating the quality factors and tunneling distances. Whilst dispersion must be considered for very high Q applications where an equidistant modal spectrum is required , it has negligible effect for biosensing applications.
We have also applied our model to the PS-CRDS microtoroidal cavity sensor  in order to determine the optimum parameters for an application of refractometric sensing. In such a scheme, the microcavity is immersed into water and a small refractive index change (δn) is introduced. This change in refractive index will influence the Qtotal of the cavity and can easily be measured via PS-CRDS microtoroidal cavity sensor . Figure 6 shows the modeling results for change in Qtotal of the cavity as a function of D. It can be clearly seen that to achieve maximum sensitivity there exists an optimum geometry for both λ = 633nm and λ = 1530nm. Due to high water absorption at 1530nm, ΔQ is lower than the one at 633nm. Since Qtotal varies with the wavelength, the optimum geometry will be different for each of the wavelength.
In order to understand that why an optimum geometry exists for a ΔQ measurement, we have plotted the individual terms of Eq. (22) as a function of the cavity geometry (Fig. 7). We have assumed Qcoupling = 0 which is a reasonable assumption as experimentally, contribution of Qcoupling can be minimized by taking measurements at multiple input power levels . Figure 7 shows that the optimum region exists close to the point of inflection of the Qtotal curve. Results in Fig. 7 show that as diameter increases, the WGM is better confined within the microcavity (see Fig. 1(a)) and so Qwgm is less influenced by the external environment, with the result that the total quality factor is limited by silica absorption.
Currently the standard experimental approach to determine the total quality factor of a microcavity is either by cavity ring down spectroscopy or Lorentzian fitting of a resonant peak. However, tunneling distance is also related to WGM radiative quality factor i.e. large tunneling distance means high Qwgm and vice versa. This suggests that with model provided here coupled with the experimental procedure outlined in , one can also extract Qwgm even when Qtotal is greater than 108. Such a measurement can possibly open up the way to a new sensing method for microcavity sensors.
We have formulated an empirical expression (Eq. (21)) for the tunneling distance of the fundamental modes of microtoroids. This expression is inspired by the analytical expression for the tunneling distances of microspheres by treating the microtoroid as a sphere of diameter (D + d). This formulation is not surprising as these equations are true only for fundamental modes which lie along the equatorial region of a microcavity. However, this treatment will not be true for higher order modes (non equatorial modes) as the curvature for a microsphere of diameter D + d will be quite different from a microtoroid with D and d as its dimensions (see Fig. 3(b)).
In summary, our finite element model, without any approximation in its master equation and coupled with a PML, not only gives accurate quality factors and tunneling distances, but can also determine the other important parameters (e.g. mode volumes) accurately for a wide range of applications based on axisymmetric microcavities such as biosensing, non-linear processes, and lasers.
We thank Prof. Jonathan Webb at McGill University for useful discussions which assisted in our understanding of PMLs. We are grateful to Prof. Andrea Armani at University of Southern California for her invaluable suggestions for improving the manuscript. We appreciate the efforts of Ashley Maker at University of Southern California, for fabricating the microtoroidal cavities used in this work. We are also grateful to Prof. Tal Carmon at University of Michigan for pointing us to tunneling distance calculations. This work is supported by the NSERC-CREATE training program in Integrated Sensor Systems, McGill Institute of Advanced Materials, Montreal, Canada.
References and links
1. F. Vollmer and L. Yang, “Label-free detection with high-Q microcavities: a review of biosensing mechanisms for integrated devices,” Nanophotonics 267–291 (2012).
2. X. Fan, I. M. White, S. I. Shopoua, H. Zhu, J. D. Suter, and Y. Sun, “Sensitive optical biosensors for unlabeled targets: A review,” Analytica Chimica Acta 620, 8–26 (2008) [CrossRef]
3. L. Maleki and V. S. Ilchenko, “Techniques and devices for sensing a sample by using a whispering gallery mode resonator,” California Institute of Technology, US Patent 6,490,039, (2002).
4. J. L. Nadeau, V. S. Ilchenko, D. Kossakovski, G. H. Bearman, and L. Maleki, “High-Q whispering-gallery mode sensor in liquids,” Proc. SPIE 4629, 72 (2002).
6. J. Barnes, B. Carver, J. M. Fraser, G. Gagliardi, H. P. Loock, Z. Tian, M. W. B. Wilson, S. Yam, and O. Yastrubshak, “Loss determination in microsphere resonators by phase-shift cavity ring-down measurements,” Opt. Express 16, 13158–13167 (2008) [CrossRef] [PubMed] .
7. M. I. Cheema, S. Mehrabani, A. A. Hayat, Y. A. Peter, A. M. Armani, and A. G. Kirk, “Simultaneous measurement of quality factor and wavelength shift by phase shift microcavity ring down spectroscopy,” Opt. Express 20, 9090–9098 (2012) [CrossRef] [PubMed] .
11. M. Tomes, K. J. Vahala, and T. Carmon, “Direct imaging of tunneling from a potential well,” Opt. Express 17, 19160 (2009) [CrossRef] .
12. O. Chinellato, P. Arbenz, M. Streiff, and A. Witzig, “Computation of optical modes in axisymmetric open cavity resonators,” Future Gener. Comp. Sy. 21, 1263–1274 (2005) [CrossRef] .
13. M. Karl, B. Kettner, S. Burger, F. Schmidt, H. Kalt, and M. Hetterich, “Dependencies of micropillar cavity quality factors calculated with finite element methods,” Opt. Express 2, (2009).
14. M. Koshiba, K. Hayata, and M. Suzuki, “Improved finite element formulation in terms of the magnetic field vector for dielectric waveguides,” IEEE Trans. on Microwave Theory Tech. 33, 227–233, (1985) [CrossRef] .
15. R. A. Osegueda, J. H. Pierluissi, L. M. Gil, A. Revilla, G. J. Villalva, G. J. Dick, D. G. Santiago, and R. T. Wang, “Azimuthally dependent finite element solution to the cylindrical resonator,”10th annual review of progress in applied computational electromagnetics 1, 159–170, (1994).
16. M. Oxborrow, “Traceable 2-D finite element simulation of the whispering gallery modes of axisymmetric electromagnetic resonators,” IEEE Trans. on Microwave Theory Tech. 55, 1209–1218 (2007) [CrossRef] .
17. S. M. Spillane, “Fiber-coupled ultra-high-Q microresonators for nonlinear and quantum optics,” PhD. thesis, CALTECH (2004)
21. C. Shi, H. S. Choi, and A. M. Armani, “Optical microcavities with a thiol-functionalized gold nanoparticle polymer thin film coating,” Appl. Phys. Lett. 100, 013305 (2012) [CrossRef] .
22. T. Lu, H. Lee, T. Chen, S. Herchak, J. H. Kim, S. E. Fraser, R. C. Flagan, and K. Vahala, “High sensitivity nanoparticle detection using optical microcavities,” Proc. Natl. Acad. Sci. (2011) [CrossRef] .
23. Y. F. Xiao, C. L. Zou, B. B. Li, Y. Li, C. H. Dong, Z. F. Han, and Q. Gong, “High-Q exterior whispering-gallery modes in a metal-coated microresonator,” Phys. Rev. Lett. 105, 153902, (2010) [CrossRef] .
24. M. I. Cheema and A. G. Kirk, “Implementation of the perfectly matched layer to determine the quality factor of axisymmetric resonators in COMSOL,” in COMSOL conference,Boston, Oct 8 2010.
25. A. Arbabi, Y. M. Kang, and L. L. Goddard, “Analysis and Design of a Microring Inline Single Wavelength Reflector,” in FIO , October 24, 2010.
26. F. L. Teixeira and W. C. Chew, “Systematic derivation of anisotropic PML absorbing media in cylindrical and spherical coordinates,” IEEE Microwave Guided Wave Lett. 7, 371–373, (1997) [CrossRef] .
27. V. V. Datsyuk, “Some characteristics of resonant electromagnetic modes in a dielectric sphere,”Appl. Phys. B54, 184–187, (1992).
28. L. A. Weinstein, Open Resonators and Open Waveguides (The Golem Press, Boulder, Colorado, USA, 1969) pp. 298.
29. J. R. Buck and H. J. Kimble, “Optimal sizes of dielectric microspheres for cavity QED with strong coupling,” Phys. Rev. A 67033806 (2003) [CrossRef] .
30. B. R. Johnson, “Theory of morphology-dependent resonances: shape resonances and width formulas,”J. Opt. Soc. Am. A 10, 2, (1993) [CrossRef] .
32. S. Scheurich, S. Belle, R. Hellmann, S. So, I.J.G. Sparrow, and G. Emmerson, “Application of a silica-on-silicon planar optical waveguide Bragg grating sensor for organic liquid compound detection,” Proc. SPIE 7356 (2009) [CrossRef]
33. M. Gallignani, S. Garrigues, and M. Guardia, “Direct determination of ethanol in all types of alcoholic beverages by near-infrared derivative spectrometry,”Analyst 118, (1993) [CrossRef] .
34. A. Quarteroni and A. Valli, Numerical Approximation of Partial Differential Equations (Springer, 1997).
35. V. S. Ilchenko, A. A. Savchenkov, A. B. Matsko, and L. Maleki, “Dispersion compensation in whispering-gallery modes,” J. Opt. Soc. Am. A 20, 157–162, (2003) [CrossRef] .
36. A. M. Armani, D. K. Armani, B. Min, K. J. Vahala, and S. M. Spillane, “Ultra-high-Q microcavity operation in H2O and D2O, Appl. Phys. Lett. 87, (2005) [CrossRef] .