Method for predicting whispering gallery mode spectra of spherical microresonators

A full three-dimensional Finite-Difference Time-Domain (FDTD)-based toolkit is developed to simulate the whispering gallery modes of a microsphere in the vicinity of a dipole source. This provides a guide for experiments that rely on efficient coupling to the modes of microspheres. The resultant spectra are compared to those of analytic models used in the field. In contrast to the analytic models, the FDTD method is able to collect flux from a variety of possible collection regions, such as a disk-shaped region. The customizability of the technique allows one to consider a variety of mode excitation scenarios, which are particularly useful for investigating novel properties of optical resonators, and are valuable in assessing the viability of a resonator for biosensing.


Introduction
Whispering gallery modes (WGMs) are produced in microresonators by electromagnetic waves that travel along the surface of the material interface. Though the WGMs are bound waves, they produce an evanescent field that extends beyond the surface of the resonator. WGMs represent an important optical phenomenon for sensing due to the sensitivity of their evanescent field to nearby entities, such as biomolecules, which break the symmetry of the field. Optical microresonators have been the subject of much recent investigation since their discovery as an important tool for biological sensing [1]. It has been shown that microresonator WGMs can be sensitive to the presence of virions, animal cells, bacteria [2] and macromolecules such as proteins [1,[3][4][5] and DNA [6,7]. This technique thus represents an important step in the development of labelfree detection technologies [8].
The extremely high Q-factors (> 10 8 ) that are possible by exciting the WGMs of a resonator often utilise a prism, or a tapered fibre coupled to a resonator [9][10][11][12][13]. An alternative approach is to use active microresonators, where the resonator contains a gain medium. In this case, WGMs may be excited using a variety of methods, such as fluorescent dyes [14], surface plasmon effects [8,15] or nanoparticle coatings such as quantum dots [16]. Here, the focus is on coupling strategies that involve active microspheres, and the Finite-Difference Time-Domain (FDTD) computational tool developed in this work will address this coupling scenario in particular. It should be noted, however, that the tool is easily extended to include fibre-coupled methods.
The use of dipole sources in the vicinity of a microresonator allows one to investigate the dependence of the optical modes on the size, shape and refractive index contrast. Efforts to characterise the optical properties of WGMs in resonators using modelling techniques include both excitation from a plane-wave beam [17][18][19], and excitation from an electric dipole source [20][21][22][23][24][25][26][27][28]. The dipole sources can also serve as an effective analogue for fluorescent dyes or embedded nanoparticles that excite WGMs in active microspheres [14].
In this article, the full three-dimensional FDTD method is simulated using the free software package, MEEP [29]. FDTD is chosen specifically for its ability to incorporate effects such as material inhomogeneities, and a variety of resonator shapes, such as shells, ellipsoids or shape deformations to a microsphere. FDTD is also suited to analysing transient optical phenomena, as each time-step is recorded, and intermediate steps are accessible.
The use of two-dimensional FDTD simulations to describe the resonance peak positions has been investigated in the case of microdisks [30]. However, the accurate prediction of the Qfactor of the resonances represents a principal challenge [31], due to the significant dependence of the Q-factor on minute characteristics of the resonator, such as surface roughness, material inhomogeneities and microscopic deformations in the shape of the resonator [15]. As a result, theoretical Q-factors evaluated through analytic models can be difficult to obtain experimentally [15,32]. The flexibility of the FDTD method is able to address this issue by incorporating geometric, material and refractive index features in a way that is intuitive and easy to implement.
The FDTD method will be used to consider a variety of flux collection scenarios, source distribution, and resonator properties. In Sec. 4, an example scenario of a polystyrene microsphere in a surrounding medium of air or water is considered. The Q-factors obtained from the simulation are then compared to those of the Chew analytic model [25]. The angular distribution of the flux density is also investigated.
Analytic models describing WGMs in the literature generally fall into two categories: models that provide only the positions of the TE and TM modes and are unable to predict the profile of the output spectrum [33,34], and models that consider the full behaviour of the electric and magnetic fields both inside and outside the resonator, for an incident plane-wave (Mie scattering) [35] or dipole source [20][21][22][23][24][25][26][27][28]. Although the simpler analytic models used in the literature are comparatively efficient to calculate, the profile of the transverse electric and magnetic (TE and TM) modes are obtained independently, for example, through geometrical arguments [33,34]. This means neighbouring modes do not interact with each other in the model. Both types of analytic model rely on the assumptions that the sphere is perfectly round with no surface roughness, it consists of a medium that is homogeneous in refractive index, and that the total radiated power is collected in the limit of long collection times. Therefore, the formulae for the power output are cumbersome, and must be re-derived for scenarios that consider different resonator shapes or inhomogeneities.
On the other hand, FDTD is easy to customise, as changes to the resonator geometry (spheroids, shells, etc.), refractive index scenario, and introducing material inhomogeneities can be accommodated without re-deriving the boundary conditions of the resonator. Since the field equations are solved at each discrete time-step, FDTD also provides access to transient optical phenomena that may not be evident in the power spectrum at long collection times. In addition, FDTD is able to consider arbitrary positions and alignments of dipole sources, including inhomogeneous distributions of dipoles placed on the surface or throughout the medium. The freedom to define a specific flux collection region at any point in space, or for any length of time, is an important feature for ongoing investigations into a range of coupling scenarios. Measuring the flux from a particular direction and aperture automatically biases the shape of the collected power spectrum. FDTD is therefore especially suitable for simulating the effect of collecting radiation through a restricted region, such as a fibre.
The principal drawbacks in the FDTD method are the computational intensity, and systematic effects due to the discretisation of space. The finite resolution leads to an intrinsic 'surface roughness' in resonator geometries. In the following simulations, the effective roughness is 20-30 nm. This effect can be incorporated into the systematic uncertainty of the calculation. An unanticipated benefit of this effective roughness is a more realistic prediction of the Q-factor of the resonator. Note also that the simulated roughness is comparable in size to that of the typical surface roughness of a physical polystyrene microsphere, which is approximately 15 nm or greater [36].
This method is able to achieve robust predictions for preselecting specific optical properties of a given resonator prior to fabrication. For example, specific modes or wavelength regions can be reliably identified, and tailored for specific tasks by altering the initial configuration of the resonator. Furthermore, the ability to scan over a wide range of parameters may lead to valuable design solutions for biosensing that would not otherwise have been found.

FDTD Method
The Finite-Difference Time-Domain technique simulates the evolution of electromagnetic fields by discretising a volume into a three-dimensional spatial lattice [37]. Maxwell's Equa- Table 1. Computing resources required for a three-dimensional FDTD simulation of a 6 µm diameter sphere excited by a dipole source with a central wavelength of 0.6 µm. The Tizard machine at eRe-searchSA (https://www.ersa.edu.au/tizard) is used in these simulations, which uses AMD 6238, 2.6 GHz CPUs. The number of CPUs, the memory (RAM), virtual memory (VM) and wall-time are listed for a variety of FDTD grid resolution values, ∆x. The resolution in the frequency domain is quoted after being converted to wavelength, ∆λ . The wall-time increases linearly with the flux collection time, which is held fixed at 0. 6  The geometry of the resonator is defined for a dielectric medium, and placed in the discretized volume, with the edges padded with a field-absorbing perfectly-matched layer.
Since the simulation includes the complete set of radiation and bound modes occurring in the system, the only approximations in an FDTD simulation involve discretisation artefacts, finitevolume effects and any assumptions of ideal material properties one might make. By making the grid size finer, the simulation approaches a converged result, and this can be quantified in the form of a systematic uncertainty in the positions of the WGMs (see Sec. 4).
The FDTD method is eminently suitable for accessing transient or emergent optical effects, as each time step is evaluated separately. However, FDTD is computationally expensive. For example, time steps totalling a few hundred wave periods frequently require up to 100 hours of simulation time, for ∼ 30 CPUs on modern supercomputers. Table 1 summarises typical computing resources for a range of resolutions, using the Tizard machine at eResearchSA (https://www.ersa.edu.au/tizard). Nevertheless, for each simulation run, a large amount of data can be extracted by accessing the electric, magnetic and Poynting vector fields at each grid point. This allows one to map out the angular distribution of the energy density for a typical flux plane (see Sec. 4).
To generate an electromagnetic current to excite the WGMs, one or several electric dipole sources are placed in the vicinity of the sphere. Consider the example of a microsphere with a diameter of 6 microns, as illustrated in Fig. 1. An electric dipole is placed on the surface of the sphere and oriented tangentially to the surface, while a circular flux collection region is defined a distance L flux , normal to the sphere, with diameter D flux . This region aggregates the flux in the frequency domain. One may vary the orientation and position of both source and collection region to build up a map of the coupling efficiencies to different WGMs.
The power output spectrum represents an important quantity for assessing the Q-factors and the wavelength positions of the excited WGMs. The total output power (P) in terms of wavelength (λ ) is obtained by integrating the projection of the Poynting vector (S ≡ E × H * ) onto a flux region of area A: for a unit vectorn normal to the collection region. The profile of the power spectrum, plotted as sphere in the x-direction, is placed so that its normal is aligned radially outward along the same axis.
In this illustration, the optical modes are excited from a tangentially-oriented electric dipole source.
a function of λ , will show sharp peaks that correspond to the positions of TE and TM WGMs. The relative heights of the peaks indicate the coupling efficiencies to different modes, which is highly dependent on the geometry of the dielectric medium, the refractive index contrast, and also the method of mode excitation, such as dipole sources of different alignments, coupling to optical fibres or fluorescent dye coatings. The power spectrum can also be normalized to the power output P 0 of the sources that would occur in the absence of the dielectric sphere. A value of P/P 0 > 1 indicates a Purcell enhancement of the dipole by the dielectric medium. The corresponding cavity Q-factors can be obtained by treating the simulation output as that of an experiment, and measuring the central wavelength position λ 0 , and the full-width at halfmaximum, δ λ [32] This formula will be used in estimating the simulation Q-factors in Table 2.

Analytic models
The simplest analytic WGM model that provides the positions of the TE and TM modes is particularly useful for identifying the mode numbers of the WGM peaks that occur in the FDTD spectrum. While this model identifies the resonance positions, which exactly match the results of Mie scattering theory, [35], this approach does not given the total electric and magnetic fields including the dipole source, and is unable to generate the profile of the spectrum. The model is defined by obtaining the resonance condition as a function of wavelength for a dielectric sphere, and follows the formalism described in Refs. [33,34]. Starting from a very general construction, a continuous electric field, oscillating with frequency ω, takes the form E(r) = E 0 (r) e iωt , for position r. This satisfies the wave equation: for a refractive index n(r), and wave-number k = ω/c, where c is the speed of light in a vacuum.
The TE modes are obtained by separation of variables in spherical polar coordinates where the angular vector function X lm (θ ) is defined in terms of the associated Legendre polynomial P m l , and contains contributions in directionsê θ andê ϕ . The function S m (r) satisfies the second order differential equation Up to this point, the electric field simply takes the general form of a field in spherical polar coordinates, embedded in a dielectric medium. Now, the resonance condition may be obtained for a dielectric sphere, by specifying continuity across the boundary in the radial component of the electric field, S m (r). For a refractive index of n 1 inside a sphere of radius R (r ≡ |r| < R), and n 2 outside, the radial function S m (r) takes a Riccati-Bessel form, which can be expressed in terms of the spherical Bessel functions of the first kind, j m , and the spherical Hankel function, h The arguments z r 1,2 ≡ n 1,2 k r are the so-called size parameters, and A m are the coefficients defined through continuity at the boundary. The resonance condition for the TE modes of a sphere of radius R is thus [34] For TM modes, the form of E 0 contains both angular and radial vector components (forê r as the unit vector in the radial direction) where the function T m (r) obeys the equation The solutions to Eq. (12) for T m (r) take a similar form to those of S m (r) above where the coefficient B m , analogously to A m , may be defined through enforcing continuity at the boundary. This leads to the resonance condition By solving for the wave number k, a spectrum of independent TE and TM modes can be obtained, and compared to the FDTD simulation. A more general type of analytic model, such as that developed by Chance, Prock and Silbey [21] and separately by Chew [20,24,25], is able to simulate the full electric and magnetic field behaviour both inside and outside the resonator. This type of model will serve as an excellent benchmark comparison for the FDTD profile spectrum. Both models were shown to be equivalent formulations of spherical resonators excited by a dipole source [24], and obtain the same resonance positions as Mie scattering [35]. The simpler analytic models, however, are still useful at identifying the mode numbers of the peaks corresponding to WGMs in a microsphere. For a tangentially (E y ) or radially (E x ) oriented dipole source on the surface of a sphere of radius R and refractive index n 1 = √ ε 1 µ 1 , the power output equations take the following form (when normalized to the surrounding medium of index n 2 ): Although this analytic model is better suited to simulating scenarios that involve dipole sources at a variety of positions and alignments, the FDTD method is able to collect flux from an arbitrary region (in the near or far field) for any length of time. As a result, FDTD may simulate the experimental condition of collecting power from a nearby optical fibre probe. This is an important point, as this collection scenario will affect the measured power spectrum profile. It is useful to be able to estimate the size of the effect due to collection region, and this will be investigated in the specific case of the angular distribution of the flux, in Sec. 4.

Results
An FDTD simulation of a polystyrene (n 1 = 1.59) microsphere, with a diameter of 6 µm is carried out for a tangentially-oriented electric dipole source, emitting a Gaussian pulse with a central wavelength of 600 nm, and a width of 5 fs 2 . A circular flux collection region with a diameter of 2.58 µm is placed a distance of 240 nm from the surface of the sphere in the x-direction, with its normal aligned along the same axis. Spectral information is then collected for wavelengths in the range 500-750 nm. Ihe finite grid resolution entails an asymmetry in the Gaussian peak, which diminishes as the resolution increases. A comparison of a variety of grid resolutions is shown in the left panel of Fig. 2, using a total flux collection time of 0.6 ps (120 times the width of the pulse). For resolutions in each spatial direction of 22 − 33 nm, one finds that as the resolution decreases, the profile features of the spectrum do not alter significantly; there is only a small offset in the positions of the peaks, and the peak heights. The positions of the WGMs are determined from a flux collection of the frequencies, which are then converted to wavelengths. The temporal resolution is further interpolated to yield a value < 1 nm.
The systematic uncertainty in the resonance positions due to resolution can be quantified by tracking the positions of the WGM peaks. By examining the wavelength positions of the most prominent peaks, one can obtain a converged result to a chosen tolerance, as shown in the right panel of Fig. 2. An error estimate in the peak positions can be obtained by comparing the results for different grid resolutions. For example, using a resolution of ∆x = 25 nm, the most central peak has a wavelength position of 602.11 nm. For a resolution of ∆x = 22 nm, the position is 601.35 nm. This yields an error estimate in the position of the central peak of (602.11 − 601.35) = 0.76 nm, or 0.13%. Note that an even greater spectral resolution (∼ 0.03 nm) is required for the detection of nearby macromolecules by a microsphere [41]. This presents an issue in simulating the effect of detection using FDTD. However, FDTD is suitable for providing a realistic estimate of the WGM Q-factors under a variety of resonator scenarios, and for assessing their feasibility for new designs of biosensors.
By changing the length of time allowed for flux collection, one can also obtain important insights into the distribution of modes in the power spectrum. The left panel of Fig. 3 shows the resultant normalized power spectrum for a variety of flux collection times. A spatial resolution of ∆x = 25 nm is chosen, and also a density of sampled wavelengths corresponding to ∆λ = 0.62 nm. As the collection time increases, the WGM peaks become more prominent. Furthermore, double-peak structures seen at small collection times are no longer apparent for larger collection times, as the sampling of the full mode structure of the radiating cavity does not have sufficient time to be measured by the flux region. For the longer collection times, a Purcell effect is observed (P/P 0 > 1), at wavelengths coinciding with the dominant WGMs. This indicates that the power output of the source is enhanced at these wavelengths by the presence of the dielectric microsphere.
The behaviour of the Q-factor as a function of collection time is shown in the right panel of Fig. 3 for three central WGM peaks. For the peak at λ ∼ 0.61 µm, convergence is achieved beyond a collection time of 1.0 ps, with a variation of 3.4%. For the other two peaks, λ ∼ 0.6 µm and 0.59 µm, variations of 16.9% and 15.3% are measured, respectively. This systematic effect will be considered when reporting the results in Table 2.
A comparison of the analytic model and the FDTD simulated spectrum for a microsphere is shown in Fig. 4 for tangentially and radially oriented dipole sources in both air and water media. The WGM positions from the analytic model are shown as vertical bands, with TE m,n modes in green and TM m,n modes in red, for azimuthal and radial mode numbers m and n, respectively. The width of each band indicates the estimate of the systematic uncertainty due to the finite grid spacing of FDTD. Specifically, the Mie scattering analytic model is used to estimate the shift in the WGM positions due to uncertainty in the sphere diameter, 6 µm ±∆x/2. In this case, the spatial resolution is held fixed at the value ∆x = 22 nm. For the temporal resolution, the spectral density of 800 points yields an uncertainty of (750 − 500)nm/800 = 0.31 nm. The two uncertainties are added in quadrature. The small offset of each peak from the position expected from the analytic model may be due to mixing of nearby modes, which interact with each other to shift the peaks away from their central positions. The TE and TM modes cannot be completely decoupled due to the spherical symmetry, and contributions from both TE and TM modes are expected in the spectrum, in addition to non-whispering gallery radiation modes.
For a tangential electric dipole (top left panel in Fig. 4), one expects that the dominant WGMs excited are the lowest-order radial TE modes since the electric field of each of the TM modes contains a radial component (see Sec  source. It is found that the dominant peaks have a free spectral range (FSR) consistent with a radial mode number of n = 1. Note that the fundamental radial modes (n = 0) cannot be resolved for this index contrast at this finite grid size, since they are known to exhibit extremely large Q-factors experimentally [42] and in the Chew model [25]. In FDTD, the discretisation of space leads to an effective surface roughness, which broadens the peaks, leading naturally to lower Q-factor estimates (see Table 2). In the case of the radial electric dipole (top right panel in Fig. 4), the dominant peaks exhibit an FSR that matches the n = 1 and n = 2 TM modes. In this case, there appears to be very little contribution from the TE modes.
For a lower index contrast scenario, such as polystyrene microspheres in a surrounding medium of water (bottom panels of Fig. 4), the WGM peaks are broadened, and there is a reduced density of modes at wavelengths in the range 500-750 nm. As a result, the peaks observed in the FDTD simulation correspond to the fundamental radial TE/TM modes, for a tangential/radial dipole source, respectively.
Comparisons of the Chew model, the Mie scattering analytic model and the FDTD simulations are shown in Fig. 5, for a surrounding medium of water. In the left panel, a tangential dipole source is used, and the fundamental radial TE modes from the Mie scattering model  (green vertical lines) exactly match the peaks of the Chew model (dashed blue line), as expected by construction [25]. The peaks of the FDTD simulation also correspond to these TE modes, and the FSR matches that of the Chew model, with a small peak offset due to the finite grid resolution. The right panel of Fig. 5 shows the result for a radial dipole. The peaks of both the Chew model and the FDTD simulation correspond to the fundamental radial TM modes.
The simulation can also provide information about the angular distribution of the flux, both directly and indirectly. Since the FDTD simulation records the electromagnetic field values, the flux density, S(r), may be projected onto the circular region of the z-y plane (with normal vector n) for any wavelength value. An analysis of this type is helpful for visualising the distribution, distinguishing the azimuthal and radial modes at different time slices, and identifying transient resonant features. By integrating the flux over a collection time of 1.0 ps, the distribution over the collection plane indicates the angular dependence of the mode at a particular wavelength. As an example, Fig. 6 shows the flux distribution over the collection region, for a tangential dipole source, and a surrounding medium of air. The four panels display the four most prominent WGM peaks occurring in that wavelength region, corresponding to values of 0.576, 0.588, 0.601 and 0.615 µm, respectively. It is apparent that the modes at wavelengths of 0.601 µm (TE 38,1 ) and 0.615 µm (TE 37,1 ) exhibit a flux distribution with peaks spread out in the flux measuring region over a comparatively wide angle. The total power output for these modes Table 2. A summary of the Q-factors and wavelength positions, λ (µm), for the four most prominent WGM peaks, for each plot displayed in Fig. 4. The scenarios considered are: a surrounding medium of air or water with a tangential (E y ) or radial (E x ) dipole source. Due to finite collection time, a systematic uncertainty of up to 14% is expected in the Q-factors.
Scenario peak 1: λ (µm), Q peak 2: λ (µm), Q peak 3: λ (µm), Q peak 4: λ (µm), Q   is likely to vary significantly if the diameter of the region is reduced. In contrast, the mode occurring at 0.588 µm (TE 39,1 ) has a more focused concentration of flux in the centre of the measuring region, corresponding to a narrower angular distribution of flux. The power output for this mode will be consistent over a wider range of flux region diameters. As a consequence of the nontrivial angular distributions of the modes, changes in the power spectrum are also detected. Figure 7 shows the spectrum for the 6 µm sphere for flux planes at different distances from the sphere (left panel), and different diameters (right panel). In each case, a tangential electric dipole is used, and the collection time is held fixed at 1.0 ps. By comparing the spectra from several flux collection plane sizes and positions, changes in the peak heights are observed, indicating the different angular behaviour of the modes. Changes in the distance of the flux region, L flux , result in a similar mode structure, but the overall power output, especially for the dominant WGM peaks, changes as the flux is sampled differently in each region. Alterations in the diameter of the flux region results in more noticeable changes to the mode heights, as the restriction of the measuring region demonstrates the highly angulardependent nature of the relative contribution of each mode to the total flux.
These scenarios indicate that the FDTD simulation package represents a promising development toward a robust, customisable, predictive toolkit for simulating the whispering gallery    modes of microspherical resonators. In this investigation, the FDTD tool has already provided guidance for future designs of resonators, which can be adapted for biosensing applications. For example, in Fig. 4, dominance of TE or TM modes is dependent on the alignment of the dipole source on the surface of the sphere. Microsphere sensors may be able to detect the orientation of an external biomolecule by recording the relative coupling strengths of TE to TM modes.

Conclusion
This work presents the first easily customisable approach for simulating Whispering Gallery Modes of microspheres, which is based on the FDTD method. The tool is easily customisable, and especially suited to investigating a wide variety of resonator scenarios. The FDTD approach is able to simulate changes to the power spectrum by measuring flux in the near-field region for any length of time, using an arbitrary flux collection plane that can resemble measurement from an optical fibre. The comprehensive inclusion of all modes, and access to transient optical properties for a range of time slices, allows for more realistic simulations than is feasible for typical analytic models. In addition, different methods of mode excitation of an active microsphere can be investigated, such as dipole source locations and alignments. The intrinsic surface roughness associated with the finite grid spacing leads to more realistic Q-factors, and changes in the peak positions of the modes due to discretisation effects can be quantified, and absorbed into the systematic uncertainty.
The tool may be extended in a straightforward manner to include different coupling scenarios, such as a tapered fibre, or a distribution of dipole sources placed on the surface of the microsphere, which can act as an effective analogue for the behaviour of fluorescent nanoparticle coatings used to excite optical modes in experiments. Different resonator shapes or material properties, such as index contrast or material inhomogeneities, may also be included without requiring re-derivation of the formulae, making the tool easily generalisable.
In practice, one can scan over a variety of parameters to gauge the effect on the profile of the collected power spectrum. Consequently, a wavelength and resonator combination can be sought systematically, and optical attributes of particular interest may be emphasized or suppressed. The optimisation of the coupling efficiency to certain mode channels allows one to preselect desired resonator properties in a cost-effective way, leading to new design solutions that might not otherwise have been found.
As a first example, the angular behaviour of the energy flux was investigated for several sizes and distances of the collection region, leading to alterations in the coupling efficiencies to different modes. This is observed in the profile of the power spectrum. The magnitude of the flux density projected onto the collection disk was also mapped out for a range of collection times, allowing one to distinguish the distribution of the modes at different collection times.
The new simulation results were compared and contrasted to two different analytic models. The analytic models rely on assumptions of a perfectly smooth dielectric, with the total radiation output measured in the far-field region. It was found that the peak positions from the models matched the FDTD simulations well.
The computational toolkit presented in this paper represents the first step in establishing a realistic and easy-to-use simulator, which is important for facilitating a cost-effective approach to designing tailored optical resonators for biosensing.