Statistics of multiple-scatterer induced frequency splitting in whispering gallery microresonators and microlasers

We investigate numerically and experimentally the statistics of the changes in the amount of frequency splitting upon the adsorption of particles one-by-one into the mode volume of whispering gallery mode (WGM) microresonator and microlasers. This multiple-particle induced frequency splitting (MPIFS) statistics carries information on the size and the number of adsorbed particles into the mode volume, and it is strongly affected by two experimental parameters, namely the WGM field distribution and the positions of the particles within the mode volume. We show that the standard deviation and maximum value of the MPIFS are proportional to the polarizability of the particles, and propose a method to estimate particle size from the MPIFS if the only available data from experiments is frequency splitting.


INTRODUCTION
In optical microresonators with circular boundaries, whispering gallery modes (WGMs) traveling along the resonator perimeter are supported. Such resonators have the ability to trap light in a strongly confined area with extremely low losses [1][2][3]. They have attracted great interests in the past few years due to their high quality factors (Q) and small mode volumes (V ), which enable ultra-high light intensity within the resonators [4]. This makes WGM resonators promising platforms to study light-matter interactions and nonlinear optics, and to achieve high-sensitivity label-free sensing of biological and chemical materials as well as nanoparticles [5][6][7][9][10][11][12]. When the evanescent tail of the WGM field interacts with the material of interest adsorbed into the resonator mode volume, the effective optical path length traversed by the WGM field along the periphery of the resonator changes, giving rise to changes in the spectral properties (e.g. linewidth or frequency) of the resonant mode. Such changes are monitored and recorded as the sensing signal and processed to extract information on the adsorbed material.
Eigenmodes of a WGM resonator are twofold degenerate with identical polarization and resonance frequency: clockwise (CW) and counter-clockwise (CCW) travelling wave modes (TWMs) denoted as a cw and a ccw , respectively. When sub-wavelength scatterers enter the mode volume, back-scattering of the WGM field from the scatterers into the cavity mode volume couples the CW and CCW modes, lifting the degeneracy (i.e., mode splitting or frequency splitting). As a result, two new orthogonal eigenstates which are standing wave modes (SWMs) appear as superpositions (i.e., (a cw + a ccw )/ √ 2 and (a cw − a ccw )/ √ 2) of the two counter-propagating TWMs [6,7,[13][14][15][16][17]. Since light-scatterer interaction strength is proportional to the light intensity at the location of the scatterer within the mode volume, two SWMs are affected differently by the scatterers, because the intensity distributions of these two SWMs are π/2 phase-shifted from each other (i.e, phase here refers to the spatial distance between the nodes of the SWMs, and the distance between two adjacent nodes corresponds to π). This is reflected in the different resonance frequencies and linewidths of the two SWMs. When frequency splitting of the two modes is larger than their mean linewidth, the two SWMs are resolved as a doublet (two resonant modes) in the transmission spectrum [18][19][20].
Scatterer induced mode splitting in a WGM resonator has been used to detect and size single nanoparticles and viruses adsorbed onto the resonator [6,7,21]. The two resonance modes forming upon splitting act as a self-referencing system [22] helping to discriminate interfering perturbations from the ones of interest and to minimize, if not eliminate, the effects of some noise sources which affect both modes on the detection and size estimation. Figure 1 shows a series of measured frequency splitting in response to adsorption of Polystyrene (PS) nanoparticles on a microtoroidal resonator. Each discrete change in frequency splitting corresponds to one particle binding event. Since WGM field distribution along the resonator surface is not uniform, the heights of discrete upward or downward steps differ for each binding event. Unlike single-side shift of resonance-shift based detection, frequency splitting either increases or decreases with binding particles [7]. As a PS particle attaches onto a silica resonator, resonances of the split modes shift to lower-frequency side (due to the higher refractive index of PS than that of air), with the shift amount determined by overlap of the particle with the field distributions of the two SWMs. If the higher-frequency mode shows a larger shift than the lower-frequency mode as shown in the lower panel of Fig. 1(b), frequency splitting decreases leading to observed downward discrete step in Fig. 1(a).
In real-life applications, it is not easy to control or measure the positions of individual particles on the resonator. Thus, the change in frequency splitting appears to be random. It would be interesting to investigate the statistics of frequency splitting in response to a large number of particle binding events. We demonstrate that statistical properties of multiple-particle induced frequency splitting (MPIFS) are related to the polarizability of the particle and thus can be used to extract the information of particle polarizability. We aim to answer the following question: Can we estimate the number and average size of particles in the mode volume of a resonator if the statistics is provided? In the following parts, we first introduce a theoretical model to carry out numerical simulations on MPIFS, and then present a comparative discussion on the results of experiments and numerical simulations of MPIFS. In this study, we consider particles much smaller than the optical wavelength and thus can model nanoparticles as dipoles such that the evanescent field is uniform over each particle. For numerical simulations and experiments, we used microtoroidal resonators as the model device but the results are valid for other types of WGM optical resonators.

THEORETICAL MODEL
Theoretical models of MPIFS based on different frameworks have been investigated in previous work [23][24][25]. In this section, we briefly introduce the model we employ for this study.

Polarizability.
In scatterer-induced mode splitting, the amount of mode splitting 2g (spectral distance between the split modes) and the additional dissipation 2Γ (linewidth difference of the split modes) are functions of the polarizability α [6,18]. Their ratio allows one to estimate the polarizability from experimental data (2g and 2Γ can be measured in the experiments) without the need for knowing the location of the particle on the resonator and the mode volume V of the WGM [6]. Therefore, it is worth discussing on how to define or assign polarizability accurately to particles. In (1)-(3), f (r) represents the normalized (i.e., normalized to the maximum value) WGM field magnitude at the particle position r, ω c is the resonance angular frequency, n m is the refractive index of the surrounding medium, and ν is the speed of light in the surrounding medium. The concept of polarizability arises in calculating the response of a particle in an electromagnetic field. In the simplest form where the particle is a small sphere and the field is uniform over the particle, this response can be calculated by assigning an induced dipole whose polarizability satisfies Clausius-Mossotti relation. In the electrostatic limit, α for a homogenous spherical particle of radius R λ (Rayleigh approximation), where λ is the wavelength of the light, is given by with n p and n m denoting the refractive indices of the particle and the surrounding medium, respectively. Clearly the polarizability depends on the volume and the dielectric constant of the particle, as well as the dielectric constant of the environment the particle is embedded in. At the Fröclich frequency which minimizes (n p 2 + 2n m 2 ), the polarizability α sphere and hence the induced homogenous polarization inside the particle experience a resonant enhancement associated with dipolar surface plasmons. For a non-absorbing surrounding, n m 2 is real, and hence the resonant condition reduces to Re[n p 2 ] = −2n m 2 . The non-vanishing imaginary part Im[n p 2 ] limits the magnitude of the polarizability at resonance. The dielectric constant (i.e., refractive index) of metal particles is strongly dependent on the wavelength (i.e., frequency ω) of the incident light (i.e., n p (ω)), and is negative for most frequencies in the visible range. Thus it can satisfy the resonance condition.
For non-spherical particles such as ellipsoid or rod-like particles, the expression in (4) is no longer valid. For example, the dipolar polarizability of an ellipsoid with principal axes a, b and c should be defined along each of its principal axes j, and is given as [26,27] πabcn m 2 n p 2 − n m 2 n p 2 + L j n m 2 (5) where L j denotes the geometrical depolarization factor satisfying L j = 1. For spherical particles, since a = b = c = R, we have L 1 = L 2 = L 3 = 1/3, implying that a spherical particle has the same dipolar polarizability along all directions. For spheroidal particles with a = b, on the other hand, we have L 1 = L 2 and the depolarization effects depend on c/a. For particles where R λ is not fully satisfied (e.g., R > 25 nm in the visible wavelength), the particle does not experience a uniform field over its volume, and the electrostatic limit cannot be justified. In this case, retardation effects must be accounted for [28][29][30]. Retardation effects take place in such particles because opposite charges in the induced dipole mode become largely separated -approximately one particle diameter-such that one end of the particle feels the changes in the other end with a phase delay due to the finite speed of light. Thus the period of the oscillations increases to accommodate this phase delay. Hence, a reduction in the depolarization field and mode energy is seen. This phase difference over the particle volume is especially significant for particles of larger size (e.g., R > 50 nm illuminated with visible light). Consequently, the expression for α given in (4) should be modified to take the retardation effects into consideration as suggested by Kuwata et al. [28] for particles of arbitrary shapes. Indeed, for large spherical particles, retardation effects can be seen in the multipole extension of the fields and the corresponding Mie coefficients. Retardation effects in large particles are manifested as red-shift of the plasmon resonance, appearance of higher multipoles such as quadrupoles in addition to dipoles, and radiative loss as the particle size increases. For small nanoparticles (e.g., R < 25 nm) in the electrostatic limit, on the other hand, plasmon frequency is insensitive to particle size.
Polarization of the light plays an important role in the light-nanoparticle interactions. The knowledge on the polarization properties of the light is crucial, in particular, for non-spherical nanoparticles, because the induced polarization depends on the component of the electric field along each of the principal axes. In contrast to spherical particles where the diameter is the same along different polarizations (e.g., transverse electric (TE) or transverse magnetic (TM)) of the incident light field, spheroidal and rod-like particles have different diameters along different light polarizations [31,32]. Therefore, for such non-spherical particles, different polarizations of light field will experience difference polarizabilities and lead to different dipole moments, even for small particles satisfying the Rayleigh limit [12,26,27,29]. Moreover, higher multipolar orders (i.e., the first order corresponds to small particle limit: dipolar approximation) which should be considered for larger particles are polarization dependent and should be taken into account. Therefore, to accurately assign polarizability to a detected particle of unknown shape, one should look at the response by varying the polarization of the light field. For the specific system we consider here, that is the interaction of the WGM field with a nanoparticle, the interacting field is an evanescent field and is not homogenous over the volume of large particles. Therefore, the electrostatic approximation may not hold anymore and the multipole expansion of the polarizability and the Mie coefficients should be considered [33,34]. It is well-known that higher multipolar orders contribute to the scattering and absorption more in the case of evanescent fields than of plane waves, and more for larger particles than the smaller particles (i.e., Mie coefficients of higher orders decrease very rapidly for small particles, thus their contribution is not significant for smaller particles) [33]. Moreover, these higher orders are polarization dependent resulting in polarization dependent scattering and absorption cross-sections which are larger for TM-polarized light [27].
In WGM resonators, the electric field component of TE modes is axial with no radial component whereas the electric field component of TM modes is predominantly in the radial direction with a relatively small azimuthal component (see Fig. 2). In addition, the electric field of TE modes is continuous along the boundary between the resonator and the surrounding in contrast to a discontinuity experienced by the radial component of the TM electric field. This difference in TE and TM modes will certainly modify the scattering and absorption properties of a particle placed within the WGM field, and the contribution of higher multipolar orders will be significantly higher for the TM mode than the TE mode leading to higher scattering and absorption cross-sections for the TM mode. As pointed out by Chantada et al. [23], scattering of WGM field by particles may lead to resonance broadening which scales with the scattering cross-section of the particle and the polarization of the field for small particles. This type of scaling for linewidth broadening is difficult for large particles because the process is a mixture of scattering from the particle and the propagation of light within the particle that can couple back into the resonator.
In this study, we consider small spherical particles such as polystyrene (PS) nanoparticles R ≤ 50 nm, sodium chloride (NaCl) nanoparticles of R = 25 nm and R = 30 nm, and gold (Au) nanoparticles of R = 15 nm and R = 25 nm. The wavelength λ of the WGM is in the 1550 nm band, which is far from the plasmon resonances of the gold nanoparticles. Thus we do not consider plasmonic enhancement of scattering and absorption. The field is evanescent outside the resonator with a characteristic length of the order of the wavelength, hence the particles are completely within the field, although the field is not uniform over their volumes. A consequence of the polarization dependent scattering cross-section of the particles in the evanescent WGM field is that TM modes are more sensitive to particle-induced perturbations and enable the detection of smaller particles than TE modes. In our experiments, the motivation has been the detection of smaller particles. Therefore, by fiber polarization controllers we choose the polarization of the input light which allows us to detect the smallest particles possible. The polarization is kept fixed throughout the study. Under these conditions, we can also safely neglect retardation effects caused by geometric depolarization and use the polarizability expression in (4).

General Model for Scatterer-Induced Mode Splitting
The simplest case is one single particle (N = 1) in the resonator mode volume. Figure 4(a) depicts the field distributions of the formed orthogonal SWMs. The particle locates at the node of one of the SWMs and the antinode of the other. With more particles adsorbing onto the resonator, the spatial field distributions of the SWMs are perturbed and redistributed to maximize the coupling strength between CW and CCW modes [25]. Spatial distributions of the two SWMs depend on polarizabilities and positions of all attached particles. As a result, both SWMs experience resonance shift and linewidth broadening with the amounts determined by how much their fields overlap with the physical locations of the particles, i.e., the mode whose anti-node is closer to the particles is affected more significantly than the mode which has more particles at its node. Figure 4(b) shows the field distributions of the SWMs when a second particle (N = 2) binds onto the resonator.
To study the effects of mode splitting quantitatively, we consider N particles entering the resonator mode volume. Resonance shift and linewidth broadening of two SWMs are given by the combined effect of all individual particles whose contribution is determined by its position relative to the SWM distributions. This can be explained by the fact that the strength of the interaction between the light and the scatterer, which determines the amount of resonance shift and linewidth broadening, is proportional to the field intensity at the particle location. Here, we denote the resonance frequency and linewidth of the initial (i.e., before the adsorption of the particle) degenerate mode as (ω 0 ,γ 0 ) and of the split modes after adsorption of N particles as (ω N ) with the superscript describing two SWMs and the subscript representing the N -th particle. We define φ N as the spatial phase difference between the first particle and the anti-node of one SWM, and ψ i as spatial phase distance between the i-th particle and the first particle. Frequency shift (∆ω N ) of two SWMs after N particles enter the mode volume are written as [25] ∆ω (1) and ∆γ (1) where 2g i = −α i f 2 (r)ω c /V and 2Γ i = α i 2 f 2 (r)ω c 4 /(3πν 3 V ) correspond to resonance shift and linewidth broadening if the i-th particle is the only adsorbed particle, α i is the polarizability of the i-th particle, f (r) represents the normalized (i.e., normalized to the maximum value) WGM field magnitude at the particle position r, ω c is the resonance angular frequency, ν is the speed of light in the surrounding medium, and V is the WGM mode volume.
Frequency splitting after the adsorption of N particles in the resonator mode volume is found from (6) and (7) as which takes its maximum when one of the split modes is maximally shifted from the degenerate mode while the other split mode experiences a minimal shift. Recall that the resonance condition for the WGM resonator requires that an integer multiple of the light wavelength is equal to the optical path length that light travels within the resonator. A change in the optical path length shifts the resonance wavelength. Thus, the shifts of the split modes are related to the change in the optical path lengths they propagate within the resonator, i.e., increase in the amount of splitting corresponds to increase in their optical path difference. For the split mode which experiences the maximal (minimal) shift with respect to the degenerate mode, the optical path length should change maximally (minimally). This is consistent with the modern interpretation of Fermat's principle that the optical length of the path followed by light between two fixed points is an extremum. Thus, we can analyze the shifts of the split modes and determine the conditions for their maximization (or minimization) by taking the derivative of the expressions given in (6) and (7) with respect to φ N . Taking the first derivatives and arranging the resultant expressions using trigonometric identities, we arrive at Setting ∂∆ω (1) N ∂φN = 0 and ∂∆ω (2) N ∂φN = 0, we find that the functions in (6) and (7) have their extrema (critical points) at the values of φ N satisfying Taking the second derivatives of (6) and (7) and inserting the expression in (12), we arrive at It is then clear that for N i=1 g i cos(2ψ i ) = 0, extrema corresponding to minima (maxima) for ∆ω (1) N (i.e., shift of the lower frequency split mode) are maxima (minima) for ∆ω (2) N (i.e., shift of the higher frequency split mode) and vice versa. Thus, one of the SWMs has maximum resonance shift with respect to the degenerate mode whereas the other one shows minimal shift, giving rise to maximum splitting. For N i=1 g i cos(2ψ i ) = 0, we find from (12) that tan(2φ N ) = ∓∞ implying 2φ N = mπ + π/2 and hence cos(2φ N ) = 0 and sin(2φ N ) = ∓1. Substituting this into (10), we obtain S N = N i=1 2g i sin(2ψ i ) . For each new particle entering the mode volume, the distributions of the two SWMs rotate so as to ensure that φ N satisfies (12). This takes place at each binding event leading to changes in mode-splitting spectrum which enable detection and size measurement of binding particles.
Single-shot size measurement of individual nanoparticles and viruses have been demonstrated by tracking the changes in the resonance frequencies and linewidths of the split modes in a passive microtoroid [6,25]. Since the difference of the linewidths and the resonance frequencies of the split modes both scale as f 2 (r)/V , their ratio g/Γ becomes independent of the particle position and the resonator mode volume, enabling estimation of the size of each detected particle without knowing its position in the mode volume and the size of the resonator mode volume. It is apparent that this size measurement method requires that the changes in the resonance frequencies and the linewidths are measured with high accuracy. However, as the size R of the particle decreases, the changes in the linewidths become too small to be measured with high accuracy as the parameter Γ scales as R 6 . Since the parameter g is proportional to R 3 , there is a region of R below which resonance frequency changes can be measured but the changes in linewidths cannot be resolved. Thus, size measurement either cannot be made or can be erroneous although the detection of the particle is still possible.
In a recent work [7], we have demonstrated that mode splitting in a WGM microlaser [8] enables enhanced detection sensitivity capable of detecting smaller particles beyond the reach of passive resonators employing mode splitting. Each of the split modes is a lasing line which when detected at a photodetector gives rise to a beat signal whose frequency is the same as the amount of mode splitting. This provides an easy way to detect splitting and hence monitor continuously the binding of nanoparticles onto the microlaser. However, it is almost impossible to perform single-shot measurement of linewidth difference of the split lasing lines as the laser linewidths are very small. Thus, although smaller particles can be detected by monitoring the beat frequency, single size measurement cannot be done due to the lack of information on the linewidth difference. In such cases, where changes in frequencies are detected, we can extract an average size for the ensemble of adsorbed particles. It is thereby important to investigate how the splitting spectra are affected as multiple particles are entering the resonator mode volume.
In Sections and , we study numerically multiple-particle induced frequency splitting in microtoroids as a function of particle number and size when taking into account particle positions in the mode volume. In all numerical simulations, we consider the WGM with distribution shown in the lower panel of Fig. 2 and a wavelength of λ = 1550 nm. The maximum field is f max = 0.36 and the mode volume is V = 280 µm 3 . We assume particles enter the mode volume sequentially and uniformly, i.e., ψ has uniform distribution from 0 to π and θ has uniform distribution from 0 to 2π/3. These assumptions are reasonable for our measurements in which particles land on different positions of the resonator with equal possibility. Moreover, in the numerical simulations we assume that all the nanoparticles in the ensemble are spherical and have the same size and polarizability. However, the analysis can be extended to ensembles where the size or the polarizability distribution of the nanoparticles follows a certain statistical distribution such as Gaussian or Poissonian.

MULTIPLE-PARTICLE INDUCED FREQUENCY SPLITTING SN
As indicated by Eq. (10), the frequency splitting S N changes with each particle entering the resonator mode volume. Figure 1 shows the experimentally obtained S N as a function of time as gold nanoparticles are consecutively deposited onto a microtoroid. Measurement details are explained in [7]. In Fig. 1, we see that S N either decreases or increases with each particle and there is no observed trend which can relate S N to the number of adsorbed particles, indicating that S N varies with particle position on the resonator. It is interesting to look at the N -particle induced S N statistically which can be obtained by a large number of repeated trials.
In numerical simulations, N particles with radius R are randomly deposited in the mode volume of a resonator, and the induced frequency splitting is calculated using Eqs. (6), (12) and (10). This process is repeated for 10000 times to obtain a statistically significant distribution of S N . Distributions of S N obtained for various values of N and R are given in Fig. 4, which depicts that the distribution of frequency splitting becomes broader (i.e., standard deviation increases) and mean frequency splitting shifts to higher values as the size of the particles and the number of deposited particles increase. The mean S N µ and the standard deviation S N σ of the distributions of S N are given in Fig. 5 as a function of N and R, from which the curve fittings reveal . These agrees well with the results in Ref. [23]. The coefficients of the linear relations are determined by the WGM field distribution f and the particle positions ψ and θ. It should be pointed out that for each single set of experiment realization, the frequency splitting is random and does not necessarily follow the curves in Fig. 5.

PARTICLE-INDUCED CHANGES IN FREQUENCY SPLITTING ∆SN
In this section, we study the change in frequency splitting in response to adsorption of individual particles. In other words, we study the statistics of the amount of change in the frequency splitting of the WGM of interest upon binding of single particles in order to estimate the size or polarizability of the particles and the number of adsorbed particles. Here, we define the change of frequency splitting induced by the binding of the i-th particle as ∆S i = S i − S i−1 (Fig.  6).
As discussed in Section , for an ideal resonator without particles, the CW and CCW modes are degenerate without frequency splitting, i.e., S 0 = 0. As the first particle enters the resonator mode volume, we have ψ 1 = 0 and φ 1 = 0.
Inserting (15) in (14), we obtain which shows the dependence of S 2 on g 1 , g 2 , and the distance between the two particles (i.e., ψ = |ψ 1 − ψ 2 | = |ψ 2 | since ψ 1 = 0). In the special case of g 1 = g 2 = g, we find φ 2 = ψ 2 /2 indicating that one of the SWMs has its anti-node at the midpoint between the two particles. Frequency splitting in this case is S 2 = |4g cos(ψ)| which is maximized (minimized) for ψ = mπ (ψ = mπ + π/2) for an arbitrary integer m. Moreover, (16) implies that S 2 varies in the range between |2g 1 − 2g 2 | and |2g 1 + 2g 2 | for varying values of ψ. Thus change in the frequency splitting ∆S 2 = S 2 − S 1 in response to the second particle position is in the range from − |2g 2 | to |2g 2 | for a large value of 2 |g 1 |. This is demonstrated by the simulation result given in Fig. 7(a). In practice, 2g 2 is unpredictable because of the unknown θ. Therefore, uncertainty of the second particle position affects the frequency splitting through both 2g 2 and ψ. In Fig. 7(b), we give the distribution of ∆S 2 for different values of 2g 2 and ψ. Statistically, for a large number of repeated tests, ∆S 2 follows some distribution determined by the distributions of 2g 2 and ψ. If the two particles are identical, we find from Eq. (16) that ∆S 2 is proportional to the particle polarizability α. These discussions can be readily extended to more particles, in which case the established frequency splitting is analogous to S 1 and the new coming particle is analogous to the second particle. In general, as the i-th particle enters the resonator mode volume, frequency splitting could either increase or decrease, and the amount of change depends on the position of the i-th particle on the resonator with its maximum possible value equal to 2g i . N particles entering the mode volume of a resonator one-by-one consecutively lead to N discrete changes in frequency splitting ∆S i , i = 1 . . . N . Figure 8 presents histograms of these changes (∆S 1 ∼ ∆S N ) for various N and R. Theoretically, for a large enough number of particle binding events, the histogram of splitting changes is symmetric around zero as shown in the bottom panel of Fig. 8(a). The maximum possible change in the histogram equals to the maximum value of 2g, i.e., ∆S max = αf max 2 ω c /V , which is achieved at θ = 0. The standard deviation ∆S σ of the histogram is proportional to α with the coefficient determined by the distribution of particle positions. These are demonstrated by the negligible impact of particle number N on the distribution of ∆S when N is sufficiently large (Fig. 8(a)), and by the broader distribution of ∆S with increasing R (Fig. 8(b)). The dependence of ∆S σ and ∆S max on α can be used to extract the information of particle polarizability and thus particle size [7]. In each set of particle deposition experiment, ∆S σ and ∆S max vary from their expected values due to the uncontrollable particle positions. We conduct numerical simulations to study quantitatively the dependence of the expectations of ∆S σ and ∆S max on particle size. We perform the simulations as follows. N particles of radius R are randomly deposited into the mode volume, and the splitting changes ∆S 1 ∼ ∆S N are calculated to obtain ∆S σ and ∆S max . This is repeated for 10000 times to calculate their mean values which are plotted as a function of R in Fig. 9. A linear dependence on R 3 is obtained for the expected values of ∆S σ and ∆S max . This linear dependence can be used to estimate the size of identical particles in an ensemble measurement. In real measurements, it is impractical to repeat experiments 10000 times for each ensemble of particles to get the expected values of ∆S σ and ∆S max . Instead, one can perform only one ensemble measurement and use ∆S σ and ∆S max values obtained in that specific measurement to estimate expectations. In such a case, the detected particle number N in the ensemble is a crucial parameter determining the accuracy of the estimation. For small N , ∆S σ and ∆S max in each measurement may vary a lot from their expectations, and thus leads to a large measurement error. However, errors can be reduced by increasing N , as shown in Fig. 10. The standard deviations of ∆S σ and ∆S max are smaller for larger N (Fig. 11). The larger the N is, the closer the estimated values to the real values are, and therefore the more accurate the particle measurement is.
In ensemble measurement of multiple-particle binding events, the particle polarizability can be estimated from the expression ∆S max = αf max 2 ω c /V if the WGM field distribution and the mode volume are known. However, due to the many supported modes in a resonator, it is difficult to decide which mode is measured. In this case, reference measurements with particles of known size can be used to estimate the polarizability or size of unknown particles.
For example, when two groups of particles of the same material but different size are deposited on a resonator, the ratio of the particle polarizability is equal to the ratio of ∆S σ or ∆S max obtained for the two groups of particles. By comparing the ratio, we can eliminate the effects of the field distribution and the particle positions on the estimation. It should be noted that estimation using ∆S max gives only the lower limit for the particle polarizability because there is always a non-zero possibility that all the observed splitting changes are smaller than ∆S max . Moreover, ∆S max is susceptible to perturbations due to contaminants on the resonator surface. In the measurements in Section , we use ∆S σ for particle size estimation.

MEASUREMENTS OF NANOPARTICLE SIZE: EXPERIMENTAL RESULTS
In a previous work, we have demonstrated detection and size measurement of nanoparticles down to 10 nm in radius using a microlaser [7] obtained by optically pumping Erbium (Er 3+ ) or Ytterbium (Yb 3+ ) doped high-Q microtoroidal resonators above the lasing threshold power levels. In the presence of mode splitting, single lasing frequency of the microlaser splits into two, which interfere and lead to a beatnote signal when detected at a photodetector of sufficiently large bandwidth [35,36]. Changes in frequency splitting are thus translated to changes in the beatnote frequency which is processed as the sensing signal. The ultra-narrow laser linewidth allows detection of small particles which will go undetected if a passive resonator was used. However, it is difficult, if not impossible, to measure the linewidth difference between the split lasing lines [39]. This, in turn, makes it difficult to extract the size information of the detected particle at a single-shot measurement. Thus, ensemble measurements of particles with identical size should be performed to assign an average size or polarizability to the detected particle ensemble. In this section, we present size measurement results for sodium chloride (NaCl) and gold nanoparticles using Er 3+ doped silica microtoroid lasers [37,38], and study the dependence of size estimation accuracy on the number of particles in the ensemble. Figure 12(a) depicts the frequency of the beat signal (i.e., frequency splitting) observed as individual NaCl nanoparticles (R = 25 nm and R = 30 nm) are adsorbed into the mode volume of the microlaser. In Fig. 12(a), each discrete change in frequency splitting indicates adsorption of a particle. The statistical distribution of the changes in the splitting is plotted in Fig. 12(b). We have rejected splitting changes which lie within the noise level of beat-note signal. Thus, the histograms in Fig. 12(b) depicts gaps around 0 value. As expected, larger particles lead to broader distributions. By setting one group of particles as reference, size of the other particle group is estimated by taking the ratio of ∆S σ from the measured dataset. Due to the missing information near 0 in the histograms, ∆S σ is obtained at different threshold values to get the particle size. Detailed explanations can be found in the Supplementary information of [7]. For the ensemble of particles with R = 30 nm (R = 25 nm), we estimate the size as R = 30.42 nm (R = 24.81 nm) when the measurement results with particle ensemble of R = 25 nm (R = 30 nm) are used as reference.
The estimation results agree with the size information provided by the manufacturer, but this single estimation result does not tell us how accurate it is as an estimate of the true value. To determine the accuracy of size estimation, we use bootstrap method to obtain the confidence interval of the estimate [40,41]. In the bootstrap approach, the measured original dataset is randomly resampled to form a new dataset having the same length as the original dataset. Each of the resampled data is obtained from random sampling and replacement of the original data points. For each resampled dataset, we calculate the particle size. This resampling process is done 1000 times giving us a collection of 1000 size estimations. The distribution of these estimations approximates the distribution of the actual particle size. The confidence interval can thus be obtained by using the appropriate upper and lower percentages of the distribution. Using this method, 95% confidence intervals for measurements in Fig. 12 are found as (21.38nm, 29.49nm) and (25.41nm, 35.05nm) for two different particle sizes.
One can reduce the estimation error by increasing the number of particles that are detected. We performed experiments with gold nanoparticles to verify it. Using a single microtoroidal laser, we detected around 400 particles of R = 15 nm and 400 reference particles of R = 25 nm. Distribution of the estimated size using bootstrap method is shown in Fig. 13(a) with mean value of 14.97 nm. To study the effect of the number of detected particles on the accuracy of size estimation, we chose the first N (1 ≤ N ≤ 400) data points from the original measured data set as the new data set to estimate the particle size and its 95% confidence interval. The results as a function of N are plotted in Fig. 13(b). Obviously, the confidence interval decreases with increasing particle number, suggesting a higher accuracy. For particle number of N = 250, the mean value of estimated size is R = 14.92 nm with the confidence interval of [13.94 nm − 15.94 nm] which is close to the data provided by the manufacture (mean 15.15 nm with the coefficient of variance of < 8%). These results imply that the proposed size measurement method does not necessarily require that all the particles in the ensemble have exactly the same size, although for the numerical simulations we considered the same size. Ensembles of particles of the same material and shape but not exactly the same size (i.e., size of the particles in the ensemble obeys a statistical distribution with a mean and standard deviation) can also be characterized with the proposed method.
Finally, we would like to discuss the similarities and differences of size measurement method proposed and demonstrated in this work and that used in the reactive-sensing scheme by Vollmer et al. [42], who have provided an analytical expression which relates the resonance shift to the particle size and the size of the microsphere resonator. In the reactive-sensing scheme with passive resonators, with the help of an analytical expression, an average size is assigned to an ensemble of identical particles using the observed maximum frequency shift after a number of particles bind to the resonator. In this method only the maximum frequency shift is used but the distribution of the all observed frequency shifts is ignored. In such estimations based on only one value chosen from the statistical distributions of the events, one may expect to have larger errors. For example, it will be always questionable whether the maximum shift that one observes after, say N binding events, is still the maximum shift after (N + 1)-th binding events, and how one can be sure that the maximum observed shift is due to binding of a single particle of interest but not due to the binding of aggregated particles or some other larger particles present in the solution. On the other hand, making use of the whole set of obtained results and their statistical distribution will help reduce false measurements due to outliers. Therefore, we expect that the method we propose in this study will provide a better size estimation.
The method proposed in this paper can be directly applied to reactive sensing schemes too. We should emphasize here that in particle detection and measurement using mode splitting in passive WGM resonators, we do not need the method proposed here. In this case, as we have shown previously [6,25], we can directly estimate the size of each particle by comparing the mode splitting spectra (changes in the frequencies and linewidths of the split modes) just before and after a single particle binding. We proposed this new method of size estimation for WGM microcavitylaser-based sensing, because in this case the ultranarrow laser linewidths are difficult to measure accurately. This method can also be used for measuring small nanoparticles which do not induce sufficiently high dissipation that can create a measurable amount of changes in linewidths. Since linewidth information is missing, we cannot measure the size of each detected particle but instead we give an average size for the ensemble using the statistics of the changes in mode splitting. In principle, reactive-sensing scheme can be realized using microlasers, too, by monitoring the shifts in the lasing frequency. However, the laser frequency shift induced by small nanoparticles is so small that one cannot use commercially available optical spectrum analyzers to monitor the shift induced by a single nanoparticle binding, and novel measurements techniques are required. On the other hand, mode splitting method is easy to use for microlaser-based sensing because the signal at the detector is a beatnote signal whose frequency corresponds to the amount of mode splitting. Microlaser-based sensing methods will certainly benefit from research targeting the detection of ultranarrow linewidths and the changes in linewidths.

CONCLUSIONS
In this paper, we have studied numerically and experimentally the statistical properties of MPIFS. We have showed by simulations and experiments that the difficulty of precisely controlling the positions of the adsorbed particles within the mode volume leads to decrease or increase in the amount of frequency splitting with varying step heights as the particles are adsorbed one-by-one onto the resonator. Despite this randomness in frequency splitting for each individual binding events, statistical analysis shows that the expected value of frequency splitting increases linearly with the square root of particle number. We have demonstrated that statistics of changes in frequency splitting can be used to extract the information of particle polarizability and hence the particle size if the refractive index of the particles is known. Although our experiments are performed with NaCl and Au nanoparticles and numerical simulations consider PS and Au nanoparticles, our results and the proposed size measurement method are valid for all types of particles. However, as mentioned before, this study does not consider plasmonic effects. This statistical analysis based size and polarizability estimation method can be used in both the resonance-shift based detection (reactive sensing) and laser frequency splitting techniques where the amount of spectral shift and change in splitting frequency are measured. A possible application area of the proposed method is the characterization of nanoparticle generators and sources. Histograms of the changes in frequency splitting induced by NaCl particles of different sizes. Particle radius is labeled in the plot. Total particle number in the histograms is 81 for 25nm-particle and 111 for 30nm-particle. Red curve is Gaussian fitting. In the measurement, we detected around 400 reference particles and 400 measured particles.
(b) Size estimation as a function of detected particle number for both reference and measured particles. Blue circles and red crosses are the mean and 95% confidence interval of size estimation obtained from bootstrap method of 1000 resampling.