Analytical Approximations of Whispering gallery Modes in Anisotropic Ellipsoidal Resonators

Numerical evolutions of whispering gallery modes of both isotropic and anisotropic spheroidal resonators are presented and used to build analytical approximations of these modes. Such approximations are carried out mainly to have the possibility to have a manageable analytic formulas for the eigenmodes and eigenfrequencies of anisotropic resonators. A qualitative analysis of ellipsoidal anisotropic modes in term of superposition of spherical modes is also presented.

where the major axis is a = 1.9 mm and the rim radius ρ = 0.25 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 b = 0.7 mm [22]. Following Ref. [4], we consider this resonator to be made of Lithium Niobate (LiNbO 3 ), an anisotropic crystal characterized by the following dielectric tensor:ε where ε ⊥ = 5.3 is the dielectric constant in the xy plane and ε = 6.5 is the dielectric constant parallel to the z axis. However, in considering the resonator as isotropic, we assume that its dielectric tensor will be diagonal, with nonzero elements that are equal to ε ⊥ . 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 Fig. 2. Here the external radius of the torus being a = 1.9 mm equal to the major axis of the ellipsoidal resonator, and the radius of the circular cross section being r = ρ = 0.25 mm, equal to the rim radius of the ellipsoidal resonator. A sketch of this approximation is shown in Fig. 3, where both the lateral and top sections of the two resonators are shown.
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 replace the real spheroid with semiaxes a and b with effective semiaxesā = a + σ andb = b + σ, 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 value: where k 0 is the vacuum wavenumber, χ = 1 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 Ref. [13]. For the rest of the paper we will implicitly assume that the resonator is closed with semiaxes a and b, remembering that the solutions we will present can be easily adapted to the open resonator by simply substituting a and b withā andb.

B. B. 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 ε. (color online) Sketch of the 3D geometry of the toroidal (left) resonator that approximates the ellipsoidal resonator and its cross section in the xz plane. The radius of the circular cross section is r = ρ = 0.25 mm, and the external radius of the torus is equal to a = 1.9 mm.
FIG. 3. (Color Online) Lateral (a) and top (b) view of the ellipsoidal resonator (blue solid line) and the approximating torus of circular cross section (red shaded circle). The radius of the torus cross section r matches the rim radius ρ = b 2 /a of the ellipsoid. The external radius of the torus matches the major axes of the ellipsoid (that appear as a circle of radius a in the top view), while the internal radius of the torus is such that its circular cross section has a radius of r.
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 xz plane thanks to the axial rotation symmetry along the z 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 [20] or the magnetic field components [21], depending on what is the best choice with respect to the problem to be solved.
Our results are based on the discussion presented in Ref. [20], 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 [17]. For this reason, in general, if one assumes that the resonator's constituent medium has constant magnetic susceptibility, then the magnetic field H 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 E.
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 either the case of isotropic (Section III.A) or anisotropic (Section III.B) 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 Ref. [18]). For all the simulations presented in this paper, the penalty constant is equal to one.ε −1 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 Eq. (4), and with the aid of Galerkin's method of weighted residuals [19] 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 Eq. (4) in its weak form counterpart, we multiply both sides of Eq. (4) by the complex conjugate of a test magnetic fieldH * 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 Eq.(4), that reads: 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 [19].
The axial symmetry of the problem suggests to describe the resonator in cylindrical coordinates {r, z, ϕ}. With this choice of coordinate system the above equation turns into a 2D equation, because the ϕ-dependent part of the magnetic field H can be factored out in the form of a complex exponential exp (iM ϕ), where M is the azimuthal quantum number, i.e.. the number of nodes of the WGM in the xy plane. The total magnetic field can be then written as H(r) = exp (iM ϕ) H r (r, z)r + H z (r, z)ẑ + iH ϕ (r, z)φ , 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 Eq. (6) and its solution to a form easy to implement in COMSOL; this procedure is clearly presented in detail in Ref. [20], so we address the interested reader to this paper. We just point out that for the isotropic case one has to put ε ⊥ = ε in Eq. (8) and following in Ref. [20] , being ε ⊥ the dielectric constant in the xy plane of the resonator and ε the dielectric constant parallel to the z-axis. In the case of uniaxial resonator, conversely, these two quantities have to be different.

A. A. Isotropic resonator
In Figs. 4 and 5 two examples of the field distribution for WGMs in the cross sectional xz plane for both spheroidal (left) and toroidal (right) resonator are shown. In particular, Fig. 4 shows the fundamental WGM characterized by the three quantum numbers l = m = 1000 and q = 1, being l, m and q the orbital, azimuthal and radial quantum numbers respectively. Fig 5 instead shows the sixth WGM characterized by having l − m = 4 and q = 2.
As can be seen from these figures, and from their x-and z-sections shown in Figs. 6-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 xz, saying that the toroidal modes well approximate the spheroidal mode is equal to say 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. Fig. 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, i.e. on the order of O(l −1 ) for the WGMs presented here. Normally one can find in literature [23] that the most common level of accuracy in determining the eigenvalues of a resonator is of O(l −1 ) 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 x and z direction. This results in a complete equivalence between the ellipsoidal and toroidal boundary condition. As can be seen from Fig. 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, i.e. they are almost identical.
When higher order modes are considered (like for example the one depicted in Fig. 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 x-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.

B. B. Anisotropic resonator
For the case of anisotropic resonators, results are reported in Fig.11 and 12 for the bidimensional xz cross sections of both spheroidal and toroidal modes and in Figs.13-16 for their sections along the x-and z-axes respectively.
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 xy plane and the one parallel to the z direction are assumed to be ε ⊥ = 5.3 and ε = 6.5 respectively, as stated in section II.A. 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 0.4%.

C. C. 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, for example, 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 xz cross sections of the spheroidal solutions. In order to describe also the rotational symmetry of the spheroid, this circumference must be rotated around the z-axis, generating a toroidal resonator like the one depicted in Fig. 2. For this resonator, the ϕ-dependent term in the eigenmodes simply consists (due to the axial symmetry) in the complex exponential of the form exp (iM ϕ), where M is the number of nodes of the eigenmode in the xy-plane, i.e. around the z-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 x-direction and trigonometric functions along the z-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 [10]. 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 is already available [11][12][13]. 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 x = kr. The solutions to this equation are the Bessel functions of the first kind J l (x), 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 Ai(x) [15]. 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 x k is the k-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 other displacement effects (β parameter). The weighting coefficients c k has to be chosen in such a way that the sum can be faithfully truncated to the order N . With this expansion, it is then possible to address the reals 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 Fig. 13 can be very well reproduced by using only one term in the expansion (8) with γ = 25 and β = 4.62. In this case the term x k = −2.3381 will be the first zero of the Airy function of the first kind Ai(x). For the more complicated pattern of Fig. 14, instead, the ellipsoidal WGM can be very well approximated by using three differently weighted Airy functions, so that N = 3, with γ = 23 and β = 4.46. The weighting coefficient in this case are c 1 = 0.15, c 2 = −0.62 and c 3 = 1.8.
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, i.e., the WGM structure along the z 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 β = lα is the fitting parameter (l 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 b k and d k 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 [14], 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, 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 Fig. 16, for which l = 1000 and l−m = 4, a very good analytical approximation can be built with the help of Eq. (10) by considering only the even part of the expansion (i.e., only cosine terms), because the mode shown in Fig. 16   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, than 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.

IV. IV. 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