Abstract
Numerical evolutions of whispering gallery modes of both isotropic and anisotropic spheroidal resonators are presented and are used to build analytical approximations of these modes. Such approximations are carried out mainly to have the possibility to have manageable analytic formulas for the eigenmodes and eigenfrequencies of anisotropic resonators. A qualitative analysis of ellipsoidal anisotropic modes in terms of superposition of spherical modes is also presented.
1. Introduction
Thanks to their extremely high quality-factor (up to ), whispering gallery mode (WGM) resonators are nowadays very promising devices for applications in optoelectronics [1], bio-sensing [2, 3], nonlinear optics [4], and other fields of applied physics. Although the first analysis of such devices dates back to 1939 [5], this topic remained unconsidered for many years. Recently, however, the need of compact optical devices with great performances (together with the possibility to manufacture not only silica microspheres but WGM resonators of any material and shape [6, 7]) had contributed to make WGMs experience a new renaissance. Unfortunately, while the theory of WGMs in microspheres is well-established, and precise calculations of eigenmodes, radiative losses, field distributions, and other physical quantities have been carefully carried out [1], when approaching the study of anisotropic devices, analytical solutions cannot be found for geometries different from an ideal sphere [8] or cylinder. In this case, it is more convenient to solve the problem with the help of numerical codes. An analytical solution, on the other side, would be preferred, as it makes easier to catch the physics behind a phenomenon and it gives the possibility to easily predict its evolution. But how to validate an analytical approximation of a solution when a real analytical solution does not exist? To answer to this question, we propose in this paper to use a numerical finite element solver (COMSOL Multiphysics [9]) to firstly find an appropriate approximation of the problem in exam and then use that approximation to check the validity of a simple analytical model for the solution.
In this paper we study whispering gallery modes (WGMs) in isotropic and anisotropic dielectric ellipsoidal resonators. Whereas nonspherical WGMs have been already studied in the past, and their analytical expression was given in terms of spherical-like functions [10, 11] or superposition of spheroidal modes [12], to the best of our knowledge no papers have been written on analytical solution for the anisotropic case. We then intend to fill this gap by finding with COMSOL a suitable approximation for WGMs in these resonators and then presenting a correspondent robust and accurate analytical solution. The approximation we present is based on the observation that since WGMs are localized near the resonator surface, an effective approximation of an ellipsoidal resonator near its surface could be represented by a toroidal resonator of circular cross section, the radius of the torus being equal to the major axis of the resonator and the radius of the circular cross section being equal to the rim radius of the resonator at the considered surface. This result gives us the possibility to substitute, under certain conditions, the complete set of spheroidal wave functions [12] that characterize the ellipsoidal resonator with the complete set of spherical wave functions [8, 13]. This permits us to manage a simple and analytical set of eigenmodes that can be easily used for theoretical predictions of nonlinear or quantum optical effects in these resonators.
This paper is organized as follows: in Section 2 the model used to simulate with COMSOL such a resonator is presented, together with a brief recall of the weak form expression for Maxwell’s equations in a rotationally symmetric resonator that is used as a model for the calculations. Section 3 is then devoted to present the results of the COMSOL simulations we made both for the isotropic and for the anisotropic cases, with a direct comparison between the field distributions and the eigenvalues of both the ellipsoidal and the toroidal resonators. Finally, in Section 4 conclusions are presented.
2. The Model
2.1. Geometry of the Resonator
Let us consider an ellipsoid of revolution filled by a dielectric medium, rotationally invariant around the -axis and characterized by the in-plane major axis , that is, the radius of the circumference in the plane , and the out-of-plane minor axis , that is, the minor axis of the ellipse in the plane . Such an ellipsoid is depicted in Figure 1, and its cartesian equation is the following: where the major axis is mm and the rim radius is mm; the value of the minor axis is then determined from these two numerical values with easy geometrical considerations and it turns out to be mm. According to [4], in fact, it is possible to obtain the value of in the following way. Let us firstly define as the angle between the vector describing the position of point on the surface of the ellipse and the major axis of the ellipse itself. Therefore, at , the rim radius of the ellipse is simply given by . This result is used here to infer, from the data presented in [4], the correct value of the minor axis . Following [4], we consider this resonator to be made of lithium niobate (LiNbO3), an anisotropic crystal characterized by the following dielectric tensor: where is the dielectric constant in the -plane and is the dielectric constant parallel to the -axis. However, in considering the resonator as isotropic, we assume that its dielectric tensor will be diagonal, with nonzero elements that are equal to .
(a)
(b)
In order to make a suitable approximation, we compare the results obtained from this resonator with the ones obtained from a toroidal resonator with circular section like the one depicted in Figure 2. Here the external radius of the torus being mm is equal to the major axis of the ellipsoidal resonator, and the radius of the circular cross section being mm is equal to the rim radius of the ellipsoidal resonator. A sketch of this approximation is shown in Figure 3, where both the lateral and top sections of the two resonators are shown.
(a)
(b)
(a)
(b)
Even if an analytical solution for the open ellipsoid (with the fields continuous at the resonator surface) does not exist, it is possible to transform the open resonator into a closed one (with ideal metallic boundaries, i.e., the electric field is zero at the resonator surface) by simply considering that the field of a WGM outside the resonator is evanescent (i.e., exponentially decaying from the resonator surface) and replacing the real spheroid with semiaxes and with effective semiaxes and , where is the depth upon which the optical field penetrates in the surrounding medium under total internal reflection. For grazing angles, does not depend on the angle and has the following value: where is the vacuum wavenumber, for quasi-TE modes, and for quasi-TM mode, being the dielectric constant of the isotropic resonator. A detailed analysis on the origin of this term can be found in [11]. For the rest of the paper we will implicitly assume that the resonator is closed with semiaxes and , remembering that the solutions we will present can be easily adapted to the open resonator by simply substituting and with and .
2.2. Weak Formulation of the Electromagnetic Problem
Before entering deep in the subject of this paper, we briefly intend to recall to the mind of the reader the weak formulation of Maxwell’s equations with Galerkin’s method of the weighted residuals, that we use to simulate with COMSOL whispering gallery modes in a general axisymmetric dielectric medium with permittivity tensor .
We choose to implement such a model rather than using a standard mode solver because the latter cannot be easily configured to fully exploit the axial symmetry of the problem and they experience some problems when dealing with WGMs. Then, in order to obtain from them an accurate solution, a fully 3D model must be implemented, making the calculations very time-consuming and complicated.
Galerkin’s method, on the other hand, has the advantage that can easily take into account the symmetry of the problem, giving in this case the possibility to reduce the dimension of the problem from 3D to 2D, by solving the problem only in the -plane thanks to the axial rotation symmetry along the -axis of the resonator.
This method, compared with mode solvers, allow us to save a lot of computational time and memory and gives us the possibility to calculate all the fields in the post processing phase of the simulation and only if we are interested in those quantities.
This method can be applied either for solving directly the electric field components [14] or the magnetic field components [15], depending on what is the best choice with respect to the problem to be solved.
Our results are based on the discussion presented in [14], where the calculation of the electromagnetic field distribution on an axisymmetric resonator is generalized to an arbitrary number of dielectric media embedded in an outer dielectric resonator, and they will be particularized in this work to study the simpler case of a resonator made by only one dielectric medium.
The electromagnetic field inside the resonator obeys Maxwell’s equations in continuous macroscopic media [16]. For this reason, in general, if one assumes that the resonator’s constituent medium has constant magnetic susceptibility, then the magnetic field will be continuous across the resonator volume, and the problem can be easily solved in terms of the magnetic field rather than the electric field .
In the general case, however, the medium is not isotropic, and the dielectric constant became a tensor. For the sake of clarity, we will approach the problem from a general point of view, leaving the dielectric constant as a tensor in order to obtain a general expression for the magnetic field inside such a resonator.
After obtaining the general expression, we will specialize it to the case of either isotropic (Section 3.1) or anisotropic (Section 3.2) resonator.
By eliminating the electric field in favor of the magnetic field in the Maxwell’s equations, after some simple algebra the problem reduces to solve the following equation: where is the so-called penalty constant that controls the presence of spurious solutions (i.e., solutions with nonzero divergence inside the resonator) that may arise during the computation. The presence of this term is quite common when a weak form differential problem has to be solved with the help of finite element codes (see, for instance, [17]). For all the simulations presented in this paper, the penalty constant is equal to one. is the inverse relative permittivity tensor (sometimes named in literature as the impermeability tensor), assumed to be independent of field strength.
In order to have the complete solution of this equation, suitable boundary conditions must be applied: in our case we will assume that the resonator is closed and so we will use the so-called electric wall boundary conditions, namely, These two equations imply that the magnetic field must be tangent to the resonator surface, while the electric field is normal to the resonator surface.
From (4), and with the aid of Galerkin’s method of weighted residuals [18], it is possible to obtain the desired weak form of the electromagnetic problem, that will be implemented in COMSOL to be solved.
In order to reduce (4) in its weak form counterpart, we multiply both sides of (4) by the complex conjugate of a test magnetic field and then integrate over the resonator’s dielectric volume.
Then we integrate by parts over the spatial coordinates and, after the disposal of the surface integrals with the help of boundary conditions, we finally arrive to the desired weak form of (4), that reads as follows: where the integral is taken over the whole resonator volume, and the first term under the integral stands for . The three terms appearing in the integrand correspond exactly to the weak form terms required to define an appropriate model within a partial differential equation solver [18].
The axial symmetry of the problem suggests to describe the resonator in cylindrical coordinates . With this choice of coordinate system the above equation turns into a 2D equation, because the -dependent part of the magnetic field can be factored out in the form of a complex exponential , where is the azimuthal quantum number, that is, the number of nodes of the WGM in the -plane. The total magnetic field can be then written as , where the imaginary unit in the -component is added to allow writing all the three components of the magnetic field in terms of an amplitude multiplied by a common phase factor. As it can be seen, the initially full 3D problem has been reduced to a 2D problem, thus resulting in a simplification of the problem and a speeding up of the computational time needed to execute the calculations.
The next step is to transform (6) and its solution to a form easy to implement in COMSOL; this procedure is clearly presented in detail in [14], so we address the interested reader to this paper. We just point out that for the isotropic case one has to put in (8) and following in [14], being the dielectric constant in the -plane of the resonator and the dielectric constant parallel to the -axis. In the case of uniaxial resonator, conversely, these two quantities have to be different.
3. Results
3.1. Isotropic Resonator
In Figures 4 and 5, two examples of the field distribution for WGMs in the cross sectional -plane for both spheroidal (a) and toroidal (b) resonators are shown. In particular, Figure 4 shows the fundamental WGM characterized by the three quantum numbers and , , , and being the orbital, azimuthal, and radial quantum numbers, respectively. Figure 5 instead shows the sixth WGM characterized by having and .
(a)
(b)
(a)
(b)
As can be seen from Figures 4 and 5, and from their - and -sections shown in Figures 6, 7, 8, and 9, the toroidal approximates the isotropic spheroidal resonator very well, resulting in the possibility to describe the set of spheroidal modes with the equivalent set of toroidal modes. Since here we are analyzing a 2D problem, and all that matters are the shapes of the modes in the cross section plane , saying that the toroidal modes well approximate the spheroidal mode is equal to saying that these modes are very well approximated by the eigensolutions of the torus.
In order to make the approximation complete, we calculated the eigenvalues corresponding to the resonator’s eigenmodes. Figure 10 shows a comparison between the eigenvalues of the first few WGMs of both spheroidal (red) and toroidal (blue) resonators. As can be seen, also in this case the approximation holds very well. For the analysis and comparison of the eigenvalues of these two cavities, we adopted as a merit parameter for this approximation to hold the mean error on the estimation of the spheroidal eigenvalues with the toroidal ones (i.e., the numerical average of the error on every single eigenvalue). This mean error takes a value roughly around 0.3% for the isotropic resonators, that is, on the order of for the WGMs presented here. Normally one can find in literature [19] that the most common level of accuracy in determining the eigenvalues of a resonator is of for WGMs. In our case this means that the precision needed for the eigenvalues must be of the order of 0.1%, that is precisely the accuracy lever of the toroidal approximation.
The approximation is very good for the first few WGMs, where the mode is highly confined near the resonator surface both in - and -directions. This results in a complete equivalence between the ellipsoidal and toroidal boundary conditions. As can be seen from Figure 10, the eigenfrequencies are very close to each other, resembling the fact that the modes have a very good overlap. This can be justified by also saying that while the modes of both resonators are highly confined near the resonator surface, they experience the same boundary conditions; that is, they are almost identical.
When higher order modes are considered (like, e.g., the one depicted in Figure 9), some discrepancies start to appear, since the mode starts to penetrate deep in the resonator. Formally, this approximation breaks when the considered mode is too extended in the resonator (namely, along the -direction) in such a way that its confining region is greater with respect to the region in which there is a good overlap between the curvature of the ellipse and the curvature of the circular cross section of the torus. In this case, the modes of the two cavities feel a different boundary (one feels ellipsoidal boundaries while the other one feels circular boundaries), and the discrepancy between these two modes starts to be significative.
3.2. Anisotropic Resonator
For the case of anisotropic resonators, results are reported in Figures 11 and 12 for the bidimensional -cross-sections of both spheroidal and toroidal modes and in Figures 13, 14, 15, and 16 for their sections along the - and -axes, respectively.
(a)
(b)
(a)
(b)
The geometry of the anisotropic resonators is the same as the one described before for the isotropic case. The value of the dielectric function in the -plane and the one parallel to the -direction are assumed to be and , respectively, as stated in Section 2.1.
Figure 17 shows instead the eigenvalue comparison for the anisotropic case for the first seven WGMs. As can be seen, even in the anisotropic case the toroidal approximation very well reproduces the eigenvalue pattern of the spheroidal resonator, and in this case the mean error was estimated to be roughly around .
3.3. Analytic Formulas
The approximations shown in the previous sections are of great interest, since they give us the possibility to approximate with a very good level of accuracy the WGM solutions of Helmholtz equations for an ellipsoid, by means of the WGM solutions of the Helmholtz equation in a toroidal cavity with circular cross section. This is very important from the theoretical point of view, because, as the WGM eigenmodes of the torus approximate well the WGM eigenmodes of the ellipsoid, then it is possible to use simple analytic expressions (as, e.g., the solutions of the Helmholtz equation for the circle) to describe both isotropic and anisotropic ellipsoidal resonators.
The results presented here show that if we take a circumference whose radius is equal to the rim radius of the actual spheroid near its boundary, the solutions for this geometry very well approximate the -cross-sections of the spheroidal solutions. In order to describe also the rotational symmetry of the spheroid, this circumference must be rotated around the -axis, generating a toroidal resonator like the one depicted in Figure 2. For this resonator, the -dependent term in the eigenmodes simply consists (due to the axial symmetry) of the complex exponential of the form , where is the number of nodes of the eigenmode in the -plane, that is, around the -axis. This ansatz highly simplifies the analytical managing of such modes in these resonators, giving the possibility to work with the complete set of circular solutions (Bessel functions along the -direction and trigonometric functions along the -direction) rather than the more complicated spheroidal wavefunctions, concurring in a huge simplification of analytical calculations in such a system. Such calculations usually represent very hard tasks, as they require to manage infinite series of spherical solutions [12]. It has to be pointed out, however, that this kind of approximation is not really needed when dealing with isotropic resonators, since a variety of analytical formulas are already available [10, 11, 20]. The present approximation, on the other hand, is of great importance for anisotropic resonators, for which an analytical expression for the eigenmodes does not exist.
Having this goal in mind, let us first consider the radial wave equation for the circle, that has the following form: where . The solutions to this equation are the Bessel functions of the first kind , as they are the only ones that are finite at the origin. Nevertheless, as usual, it is possible to approximate the Bessel functions of large argument with the correspondent Airy functions [13]. By exploiting these large-limit solutions, it is possible to build an analytical approximation to the ellipsoidal eigenmodes in the form of suitable combinations of Airy function as follows: where is the th zero of the Airy function and the quantities and are two scaling and fitting parameters, respectively, that are used to take into account both the effects of anisotropy ( parameter) and the other displacement effects ( parameter). The weighting coefficients has to be chosen in such a way that the sum can be faithfully truncated to the order . With this expansion, it is then possible to address the real WGMs of the ellipsoidal resonator with a finite number of solutions of the circle with large argument.
With these expressions for the WGMs it is then possible to fit the ellipsoidal modes. For example, the fundamental WGM described in Figure 13 can be very well reproduced by using only one term in the expansion (8) with and . In this case the term will be the first zero of the Airy function of the first kind . For the more complicated pattern of Figure 14, instead, the ellipsoidal WGM can be very well approximated by using three differently weighted Airy functions, so that , with and . The weighting coefficients in this case are , , and .
Proceeding along this way, we have then the possibility to approximate, with a large degree of precision, all the WGMs of the ellipsoidal resonator (either isotropic or anisotropic) with a finite number of eigensolutions of the Helmholtz equation in the circle, rather than using infinite series.
This is of great advantage because it guarantees more manageable expressions and allows one to easily understand the physics that lives behind these modes, as one can interpret them as just a linear combination of few modes from a circle.
In the same manner it is possible to analytically approximate the angular distribution of such modes, that is, the WGM structure along the -direction. For this case we have the possibility to employ the solution to the angular equation for the circle: by writing again the angular part of the mode of the ellipsoid as a finite sum of trigonometric terms: where is the fitting parameter ( being the orbital quantum number that appearing in the angular equation) that, again, takes into account a possible anisotropy even along the angular distribution and and are weighting coefficients.
It has to be pointed out, however, that while for an uniaxial anisotropy in a sphere the anisotropy will affect (at the first order) just the radial part of the mode [8], in the case of an ellipsoid, we must account for an anisotropy fitting parameter even in the angular part of the mode since, as the mode index grows, the mode starts to extend itself into the resonator, breaking the toroidal approximation and then starting to feel ellipsoidal boundary conditions rather than the circular ones. This can be considered and modeled as an effective anisotropy that acts prevalently along the angular direction.
For example, for the sixth mode of the anisotropic ellipsoidal resonator as the one shown in Figure 16, for which and , a very good analytical approximation can be built with the help of (10) by considering only the even part of the expansion (i.e., only cosine terms), because the mode shown in Figure 16 possesses even symmetry with respect to the vertical axis. We then have and the expansion coefficients are , , , and ; the effective anisotropy parameter takes the values and since we consider only angular function with even parity.
As can be seen, the value of the effective anisotropy parameter in the argument of the angular functions is very low but still different from zero. This means that while for a spherical anisotropic resonator (if the anisotropy is small) the angular function is not affected by the anisotropy, in the case of the anisotropic ellipsoidal resonator, this is also true for the low order WGMs. As the order of the WGM grows, then an effective anisotropy (arising from the different shapes of the spherical and ellipsoidal boundaries that the mode sees) starts to play an important role even in the angular distribution of the mode.
4. Conclusions
In this work we have used a numerical finite element solver like COMSOL to find the field distribution of WGMs in a spheroidal resonator. We have proposed to approximate the spheroidal resonator near its boundaries with a toroidal resonator with circular cross section, having the radius equal to the major axis of the spheroid and the radius of the circular cross section equal to the rim radius of the spheroid neat its surface. We presented results for both isotropic and anisotropic resonators and we pointed out that, especially for anisotropic resonators, thanks to the numerical verification of our approximation, it is possible to approximate WGMs of a spheroid with a suitable superposition of modes of the circle. The advantage of our proposed analytical approximation is two-fold: from one side it gives the possibility to easily represent the ellipsoidal WGMs as superposition of few circular eigenmodes rather than infinite series of spherical harmonics. On the other side, for the case of anisotropic resonator, it gives easy manageable formulas to describe modes that otherwise will not have an analytical representation. This could be of great help when challenging complicated calculations on such resonators because at first it allows us to use analytical (i.e., simple and practical) methods to have a first educated guess about how things should go inside the considered system. Second, a posteriori, availability of simple approximate formulas permits very powerful qualitative analysis of the quantitative numerical results. Last, but not least, WGMs found, in the last decade, many different applications, ranging from sensing [21] to nonlinear optics [22, 23], cavity QED [24], and lasing science [25]. We therefore believe that the analytical approximations described here can be used as an important guideline in each of these fields of research to have a better understanding of the physical processes that rule such applications, as well as the fact that they can constitute an educated guess for numerical simulations of WGMs, thus contributing to make them easier and faster.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.