Perturbation theory of nearly spherical dielectric optical resonators

Dielectric spheres of various sizes may sustain electromagnetic whispering-gallery modes resonating at optical frequencies with very narrow linewidths. Arbitrary small deviations from the spherical shape typically shift and broaden such resonances. Our goal is to determine these shifted and broadened resonances. A boundary-condition perturbation theory for the acoustic vibrations of nearly circular membranes was developed by Rayleigh more than a century ago. We extend this theory to describe the electromagnetic excitations of nearly spherical dielectric cavities. This approach permits us to avoid dealing with decaying quasinormal modes. We explicitly find the frequencies and the linewidths of the optical resonances for arbitrarily deformed nearly spherical dielectric cavities, as power series expansions by a small parameter, up to and including second-order terms. We thoroughly discuss the physical conditions for the applicability of perturbation theory.


I. INTRODUCTION
In this work we aim at determining frequencies and linewidths of electromagnetic resonances of nearly spherical dielectric cavities of arbitrary size and (small) deformation. In a spherical dielectric cavity, light can excite whispering gallery modes (WGMs) and circulate about any great circle with small attenuation [1]. In fact, propagation of light along a curved interface between two different dielectric media, is an intrinsically lossy process [2]. This implies that the resonant optical frequencies associated with the WGMs, have small but finite linewidths. The ratio between the frequency of a mode and its linewidth is proportional to the optical quality factor Q of the mode. This quantifies the number of optical cycles the light in the mode will stay confined within in the cavity. Values of Q around 10 11 have been achieved for silica microspheres [3]. Deviations from the spherical shape change frequencies and linewidths of the modes, thus modifying their quality factors by an amount depending on the size and the shape of the deformation. This may be either detrimental or, conversely, very useful for many applications, ranging from WGMs lasers [4], to dielectric microcavities [5]. Therefore, it is highly desirable to have a theory predicting, at least with a certain level of approximation, the frequencies and the linewidths of the electromagnetic resonances of nearly spherical dielectric cavities [6].
In principle, determining such resonances is a conceptually simple boundary-value problem: One must solve Maxwell's equations for the fields inside (medium 1) and outside (medium 2) the cavity, and match these fields at the interface between the two media. However, satisfy- * andrea.aiello@mpl.mpg.de ing boundary conditions on interfaces of arbitrarily complicated shape, is typically a formidable algebraic task. The literature about techniques and methods developed for solving this problem is, without exaggeration, enormous. Among the books we found particularly useful the Stratton's and Jackson's classical texts [7,8], and the perhaps less known but not less valuable works of Grandy [9] and Kristensson [10]. One of the first perturbation approaches to the scattering of electromagnetic waves by dielectric media of arbitrary shape was given by Yeh [11]. This study was further developed and improved by Erma [12]. Later, an important contribution to the perturbation theory of quasinormal modes in open systems was given by Lai et al., [13]. However, a serious problem with perturbation theory, based upon the analogy between the refractive index in electromagnetism and the potential energy in quantum mechanics, arises from the discontinuity of both the refractive index and the normal component of the electric field, occurring at the interface between the resonator and the surrounding medium [14]. Several methods have been proposed to deal with this issue; see, e.g., [15,16].
A different approach that avoids this problem, is the so-called boundary-condition perturbation theory. Recently, Dubertrand et al., presented a boundarycondition perturbation theory for two-dimensional disk resonators [17], which may be seen as an extension to electromagnetic waves of the classical work by Rayleigh for acoustic membranes [18]. Such theory was further developed by Wiersig and coworkers, but still limited to two-dimensional resonators [19,20].
The purpose of our work is to develop a perturbation theory, for the electromagnetic resonances of threedimensional nearly spherical dielectric resonators. As we will see, this requires us to fully account for the vector nature of the electromagnetic field in three dimensions, which is a nontrivial technical challenge. How-ever, although we deal with an effectively open system, our method permits us to avoid the use of quasinormal modes [21,22], and similar techniques [23]. We note the importance of including second-order terms in the theory. In fact, under certain conditions, first-order perturbation theory does not account for the effects of random deformations, which typically averages to zero. However, such deformations manifest nonzero correlations, the effects of which are always disclosed by second-order perturbation theory [24,25]. Further details on the theory, omitted here for brevity, can be found in [26].
The work is organized as follows. In Sec. II we establish the notation and we review the classical Mie solution [27] for the scattering of electromagnetic waves by dielectric spheres. This is functional to the perturbation theory to be developed because the Mie solution will be taken as the zeroth-order approximation. In Sec. III we establish the exact equations for our boundary-condition problem. From Sec. IV to Sec. VIII we thoroughly develop our degenerate second-order perturbation theory, including the case of highly symmetric problems where the degeneracy is not lifted to first order. In Sec. IX we apply our theory, to the simple case of an oblate spheroid resonator. In Sec. X we summarize our work and draw some conclusions. Three appendixes provide some detailed calculations.

II. NOTATION AND SCENARIO
In this section we show how to calculate the optical resonances of a dielectric sphere using the method of Debye potentials [28]. The sphere has radius a and refractive index n 1 and it is surrounded by a medium of refractive index n 2 < n 1 (typically air or vacuum). Both the sphere and the surrounding medium are nonmagnetic, homogeneous and isotropic. In the remainder we will benefit from the following definitions: • c 0 is the speed of light in vacuum.
• k 0 is the (real-or complex-valued) wave number of light in vacuum. • c α = c 0 /n α is the speed of light in a medium of refractive index n α , with α = 1, 2. • k α = k 0 n α is the wave number in a medium of refractive index n α , with α = 1, 2. • The time-independent Debye potentials u E α = u E α (r) and u M α = u M α (r) are scalar fields that, in a medium of refractive index n α , satisfy the Helmholtz equation with α = 1, 2 [29]. • The orbital angular momentum differential operator is defined as Deformation Geometry of a section of the spherical (dark blue) and of the nearly spherical (light blue) resonators. Hereêr is the radial unit vector, and n(θ, φ) is the vector normal to the deformed surface (42). From (49) it follows that cos γ =êr · n/ |n|. Figure 1 illustrates our working scenario. Without loss of generality, in this work we consider observable monochromatic electric and magnetic fields, denoted by E α (r, t) and B α (r, t), respectively, defined by where ω = k 0 c 0 = k 1 c 1 = k 2 c 2 . In a nonmagnetic medium of refractive index n α , the time-independent electric and magnetic fields E α (r) and B α (r), can be written in terms of the two Debye potentials u E α (r) and u M α (r), as [28,29] In the standard jargon, u E α and u M α yield Transverse Electric (TE), and Transverse Magnetic (TM) waves, respectively [29]. Note that from (5) it follows that where the square brackets denote functional dependence. Using the completeness and orthogonality of the spherical harmonics Y lm (θ, φ) [8], and the spherical coordinates (r, θ, φ) with r ∈ [0, ∞), θ ∈ [0, π] and φ ∈ [0, 2π), we can write, with σ = E, M and α = 1, 2. Here and hereafter l,m is shorthand for and for any smooth function s(θ, φ), The radial dependence of the form k α r, is a direct consequence of (5). Substituting (7) into (5), we obtain where the prime denotes the derivative with respect to the argument k α r, and the three vector spherical harmonics Y lm (θ, φ), Ψ lm (θ, φ) and Φ lm (θ, φ) are defined as [30]: withê r = r/r. In general, the functions U E αlm (k α r) and U M αlm (k α r) are expressible as linear combinations of spherical Bessel functions j l (k α r), h (1) l (k α r), and h (2) l (k α r), where j l (k α r) is finite at r = 0, and h (1) l (k α r) and h (2) l (k α r) describe outgoing and ingoing spherical waves, respectively, for r → ∞ (see, e.g., Appendix A of [31]). However, the electric and magnetic fields inside the sphere must be finite everywhere for 0 ≤ r ≤ a. Moreover, we assume that the field outside the sphere is made of outgoing waves only. This implies that we can write the radial parts of the four Debye potentials u E 1 (r), u M 1 (r) and u E 2 (r), u M 2 (r), in the two media as where the radial functions have been defined in terms of the spherical Bessel functions for the fields in media 1 and 2, renamed as The choice of the denominators in (13) just fixes an arbitrary normalization which could be absorbed into the definition of the coefficients a σ αlm . Substituting (12) into (10) we obtain, after a straightforward calculation, where we have defined The numerical coefficients a E 1lm , a E 2lm , and a M 1lm , a M 2lm , are determined by imposing the electromagnetic boundary conditions on the surface of the sphere [8]: It is not difficult to see that using the relationŝ with respect to x to find the resonant wave numbers for both TE and TM waves. Thus, we obtain two countably infinite sets of solutions denoted by and with σ = E, M . In the remainder, we will refer to (30) and (31), as the unperturbed spectrum. A portion of the spectrum of TE and TM resonances of a dielectric sphere with refractive index n 1 = 1.5, surrounded by vacuum with n 2 = 1, is shown in Fig. 2. FIG. 2. Spectrum of the TE and TM resonances of a dielectric sphere of radius a and refractive index n1 = 1.5, surrounded by vacuum with n2 = 1. The values of x ln = k ln a for 1 ≤ n ≤ 10 and 1 ≤ l ≤ 10 are shown as blue bands. The vertical position of the center of each band is equal to Re(k ln a), and the thickness is equal to Im(k ln a). For the first radial mode n = 1 (lighter blue bands) the imaginary part of k nl a quickly decreases as l increases from left to right. Each resonance characterized by the pair of azimuthal and radial numbers l and n is 2l + 1 times degenerate.
By construction, from the (2l + 1)-fold degeneracy of x σ ln it follows that any linear combination of the form l m=−l a σ lm u σ α|lmn (r) where a σ lm are arbitrary numerical coefficients, is still an admissible Debye potential associated with the same eigenvalue x σ ln . Evidently, for a given x σ ln , it is possible to build 2l +1 linearly independent of such combinations. We will make extensive use of this property when developing a degenerate perturbation theory.

III. PERTURBATION OF THE BOUNDARIES
Now we turn to the more general problem of a dielectric electromagnetic resonator of refractive index n 1 , surrounded by a medium of refractive index n 2 < n 1 . Again, both the resonator and the surrounding media are nonmagnetic, homogeneous, and isotropic. The expressions (5)-(7) and (12)- (16) for the fields and Debye potentials are perfectly general, so they remain valid also in the present case. What will change are the boundary conditions (17) that will be replaced by (41), defined later.

A. Describing the deformation
Consider a nearly spherical dielectric resonator, the surface of which can be described by the equation where with h(θ, φ) an arbitrary, smooth, single-valued function of θ and φ defined on the unit sphere S, which describes the deformation of the resonator. The slight departure from the spherical shape is guaranteed by any function h(θ, φ), the maximum value of which is much less than 1 on S: However, as we will see soon, this is not the only condition required for the applicability of the theory. As usual in perturbation theory, it is useful to introduce a formal parameter 0 ≤ ε 1 defined by where |f (θ, φ)| ≤ 1. This parameter is just a mathematical device used to rewrite h(θ, φ) in a more convenient form for later developments. At the end of the calculations we will restore the physical deviation h(θ, φ) by replacing everywhere εf (θ, φ) with h(θ, φ). With this definition (37) becomes The standard electromagnetic boundary conditions on the surface of the dielectric body can now be written as: where f = f (θ, φ), and the vector n = n(θ, φ) = ∇F (r, θ, φ) normal to the surface of the dielectric at r = a 1 + εf (θ, φ) , is given by where Let ∆ = ∆(r, θ, φ) denote either E 1 − E 2 or B 1 − B 2 . Then, we can rewrite the boundary conditions (41) in the suggestive form where f = f (θ, φ). This expression is exact;as no approximations have been done up to now. However, we have written it in such a way to isolate the unperturbed first termê r × ∆ (a, θ, φ), coincident with (17), from the perturbed second term delimited by curly brackets. To develop a meaningful perturbation theory we require this second term to be O(ε) with respect tô e r × ∆ (a, θ, φ). This is certainly true for the tangential partê r × ∆ (a + a εf, θ, φ) − ∆ (a, θ, φ) because by definition However, for the radial part we have and |n | is potentially unbounded. To show this, suppose that, for example, in the neighborhood of the direction (θ, φ) the deformation of the resonator could be described by with p > 0. This implies that Clearly, ε p can be of the order of unity or bigger if p ≥ 1/ε. In this case the angle γ(θ, φ) between the vectorê r normal to the unit sphere S along the direction (θ, φ), and the vector n(θ, φ) normal to the deformed sphere in the same direction, defined by can be arbitrarily close to π/2. When this occurs, we have the so-called strongly winding boundaries [20]. If this is not the case, then we have weakly winding boundaries. Through this work, we assume that the latter condition is always verified. Thus, Eqs. (38) and (46) imply that the physical conditions for the applicability of the perturbation theory, require that the magnitude of the deformation function h(θ, φ) and its gradient ∇h(θ, φ), both evaluated on the unit sphere S, must be O(ε), namely

B. Developments
To proceed further, we must rewrite (41) in a more convenient form. Note, however, that in this section we will not assume the smallness of the deformation as in (38). Following Sec. II of [12], we want to show that of the six equations (41), only four are independent. Let us write ∆ = ∆ rêr + ∆ θêθ + ∆ φêφ ≡ ∆ rêr + ∆ , and n =ê r − n θêθ − n φêφ ≡ê r − n , where, ∆ denotes either E 1 − E 2 or B 1 − B 2 , evaluated at r = a(1 + εf ). Calculating n × ∆ = 0 we obtain the three equations It is easy to see that when (51b) and (51c) hold true, then (51a) is automatically satisfied, because Therefore, choosing (51b) and (51c) as independent equations, and multiplying (51b) byê φ and (51c) byê θ , we can rewrite the four independent boundary conditions (51b) and (51c) as where all the fields are evaluated at r = a 1 + εf (θ, φ) . We can expand Eqs. (53) in terms of Ψ l m (θ, φ) and Φ l m (θ, φ) solely, the radial components being absent, to obtain with X = Ψ, Φ. For reasons that will soon be clear, on the left-hand sides of (55), we have made explicit the dependence on x = k 0 a. This arises from the radial dependence of the fields evaluated on the surface of the dielectric resonator: After integration with respect to the angular variables θ and φ in (55), we are left with the dependence on x only.
Formally, at this point our problem is perfectly posed: All what we have to do is determine the values of x (namely, the resonant wave numbers) such that the four equations possess nontrivial solutions for the coefficients a E 1lm , a E 2lm and a M 1lm , a M 2lm . These coefficients enter in (55), via the expressions of the electric and magnetic fields written in terms of the four Debye potentials u E 1 (r), u M 1 (r) and u E 2 (r), u M 2 (r), defined by (7) and (12). Needless to say, solving the system of nonlinear algebraic equations (57), is a formidable task. However, in principle, the way to proceed is direct: Substituting (15a) and (15b) into (55), after some long but straightforward calculations we can write explicitly the four equations (57) as where we have defined the matrix elements of types A and B as, respectively, with α = 1, 2 and X = Ψ, Φ. In (62) we have defined where e (θ, φ) is defined by (43), and Note that in (64) and (65) we have used (16) to write with α = 1, 2 and W = Ψ, Φ, Y . It should be noticed that in Eqs. (58) and (61) the pair of indices l , m comes from (57) and the sum with respect to l, m, originates from the expressions of the fields (15a) and (15b). Moreover, for reasons that will soon be clear, it is instructive to rewrite these four equations in the suggestive matrix form where we have defined the 4 × 4 matrix M l m lm and the 4 × 1 vector ψ lm as and respectively. Note that the incognita x in this system is the same for all the (infinite) terms of the sum with respect to l, m. Therefore, it is not possible to solve each matrix equation of the sum independently. This is why in the remainder we will develop a perturbation scheme to solve (67) with respect to x. From a physical point of view, the dependence of x on all indices (l, m) denotes the coupling between all the modes of the resonator, due to the departure from the spherical shape. We remark that the homogeneous linear system (67) is exact, and is valid irrespective of the shape and the magnitude of the deformation, and of the size of the resonator. In principle, it contains all the information about the resonances of the deformed resonator. Were we able to solve it numerically, we would not need to develop a perturbation theory. However, this is not the case.

C. The unperturbed problem
As a first step towards a perturbation theory, we must verify that the system (67) reduces to Eqs. (25) and (26) for ε → 0, that is, when the resonator is perfectly spherical. In this case, from (42) it follows that n = 0 and Eqs. (63) become where Eqs. (16) have been used and we have defined with α = 1, 2. Substituting (70) into (62) we readily find Inserting these values into (67), we obtain the algebraic system where g αl = g αl (n α x) with α = 1, 2. The block-diagonal form of this matrix equation reveals that for a spherical dielectric resonator the TE and TM waves are uncoupled. Therefore, the 4 × 4 system (73) naturally splits into two independent 2 × 2 systems, which are for TE waves and   g 1l (n 1 x) for TM waves. The first system (74), possesses the nontrivial solution a E 1lm = a E 2lm if and only if −n 1 g 1l (n 1 x) + n 2 g 2l (n 2 x) = 0. Using (71) it is easy to see that the last condition is equivalent to (27). Similarly, the second system (75), admits the solution a M 1lm = a M 2lm , provided that g 1l (n 1 x)/n 1 − g 2l (n 2 x)/n 2 = 0. Again, from (71) it follows that this condition is equivalent to (28).
We have thus demonstrated that our system of equations (67), correctly reproduces the well-known set of equations for the electromagnetic resonances of a dielectric sphere.

A remark
The matrix (68) looks formidable. However, it has in fact a quite simple structure and admits a clear physical picture. To show this, let us introduce the shorthand M = M l m lm , and omit the indices l , m and l, m everywhere. Then, it is not difficult to see that we can rewrite (68) as a block matrix: where the matrix-valued function T is defined by and Hence, M has only eight different elements of four different types, four types per each of the two media labeled by α = 1, 2.
From (73) it follows that at ε = 0 (spherical resonator), This means that at ε = 0, M 0 describes the TE resonances of a perfect sphere. The rest matrix M 0 − M 0 | ε=0 gives their corrections due to self-coupling between TE modes. The same reasoning remains valid if we replace M 0 with T (M 0 ), and TE waves with TM waves. The off-diagonal matrices V and T (V ), evidently yield the coupling between TE and TM waves due to the departure from the spherical shape, for they vanishes at ε = 0:

IV. QUANTUM-LIKE PERTURBATION THEORY
The main goal of this work is to study how the electromagnetic vibrations of a dielectric resonator are affected by a slight departure from the exact spherical form. Such a departure is quantified by the small parameter 0 ≤ ε 1 defined by (39). Equations (29) define the resonances of the unperturbed physical system, which is a dielectric sphere of radius r = a. Let us denote by x (0) any solution of either f TE l (x) = 0 or f TM l (x) = 0, the type of wave being irrelevant for the following discussion. We assume the existence of a neighborhood of ε = 0 where the algebraic system of equations (67) possesses a nontrivial solution for x = x(ε), such that Following the classical Rayleigh's scheme of perturbation theory [18], we would like to determine x (1) from a set of first-order equations in ε, x (2) from a set of second-order equations in ε, and so on. To achieve this goal in a systematic and direct manner, we find it convenient at this stage to adopt a quantumlike notation to represent the linear system of (nonlinear) equations (67). This is possible because we can always associate a linear operator with a matrix and vice versa. However, we remark that in this work the quantum formalism is just a useful notational tool that permits us to solve an entirely classical problem.
A. Linear algebra in quantum-like notation To begin with, let us introduce the fictitious vector states |l, m with l = 0, 1, · · · , ∞, and m = −l, −l + 1, · · · , l. By hypothesis, they are orthonormal, and form a complete basis in an infinite-dimensional whereÎ ∞ is the identity operator in E ∞ and here and hereafter the circumflex will mark operators in infinitedimensional Hilbert spaces. We remark that the vector states |l, m are artificial in the sense they do not represent either the scalar spherical harmonics Y lm (θ, φ) or the vector spherical harmonics (11). They are a mathematical tool that we use to solve our problem in an efficient way. Next, we define the four basis vectors |i with i = 1, 2, 3, 4. We assume that they are orthonormal, and span a four-dimensional Hilbert space, denoted by E 4 , where they form a complete basis: with I 4 the 4 × 4 identity matrix in E 4 . Then, the tensor product Hilbert space is by construction spanned by the vectors By definition, the completeness relation for E reads where here and hereafter Equipped with this paraphernalia, we can rewrite (67) as follows. First, we introduce the vector state |ψ ∈ E , such that where, using (69), Second, we define the operatorM =M(x) via the matrix elements l , m , i|M|l, m, j = i| l , m |M|l, m |j where, according to (68), Finally, using this notation we can rewrite (67) aŝ This can be easily proven by multiplying this equation by l , m , i| from the left and using the closure relation (90), where (93) and (95) have been used. We make an important remark: Unlike the case of quantum mechanics, here we have no guarantee that the operatorM is Hermitian. As a matter of fact, in general, it is not. This is why, in the remainder, we will make extensive use of biorthogonal bases generated by the right and left eigenvectors of non-Hermitian operators. The ultimate reason for the presence of non-Hermitian operators in our theory, is that dielectric resonators are intrinsically leaky systems.

B. Formal expansion
Before diving into the development of a rigorous perturbation theory, in this subsection we provide a general outline of the theory, irrespective of the precise form of the resonances spectrum.
The goal is to solve (97), here rewritten aŝ whereM (ε) is defined by (95), and has been defined in (83), with x (0) = x(0). Moreover, we assume that also the operatorM (ε) and the vector |ψ(ε) can be expanded as power of ε aŝ and where, by definition, Substituting (101) and (103) into (99), we obtain M (ε) |ψ(ε) All the terms proportional to the same power of ε must sum to zero. Thus, we obtain the chain of equations, etc. To solve iteratively these equations, we must first choose the initial state |ψ (0) (actually, the initial set of states) associated with the unperturbed eigenvalue x (0) . We will take for x (0) one of the 2l + 1 times degenerate solutions of (32), that is x (0) = x σ ln . We will see that such a solution is associated with a degenerate subspace of dimension 2l + 1. However, before starting to solve (106), it is useful to illustrate some general properties of the operatorM.
possesses some general properties which are key to the development of the perturbation theory. These properties are proven in Appendix A. In this section we will present the plain results, which are summarized bŷ where we have defined with n ≥ 1. We remark again that all the operators in (108) are not necessarily Hermitian. Note that the operatorD is independent of the order index n. BothD (0) andD are diagonal with respect to the basis |l, m , that is, l, m, i|D|l , m , j = δ ll δ mm i|D l |j , where the operators D with g αl = g αl n α x (0) , (α = 1, 2), given by (71). Moreover, Eq. (A19) gives where the prime denotes the derivative with respect to the argument: [g αl ] = dg αl (u)/du| u=nαx (0) . Note that it is possible to rewrite D l as In practice,V (n) may (and, in general, it will) depend on x (0) , x (1) , · · · , x (n−1) , but not on x (n) .

V. ZEROTH-ORDER EQUATION
In this section we will start a systematical analysis of (106), solving the chained equations order by order.
Using twice the resolution of the identity (90) and Eqs. According to (93), in the remainder we will also occasionally use the more compact notation where Since the vectors {|l, m, i } form a complete basis in E , then (114) is satisfied when all the coefficients of the expansion (114) are identically zero, that is when In matrix form this equation reads where g αl = g αl n α x (0) , (α = 1, 2) and (111) has been used.
We have already solved this system of equations in Sec. III C, and we have found two different results for TE and TM waves. Therefore, also now we will consider these two cases separately.

A. TE waves
Let us choose a pair of values l 0 , x In this case (118) becomes where and Equation (120) turns into an identity for where ψ l0m are, at this stage, arbitrary numbers. By definition, this solution (123) is valid only for l = l 0 . However, Eq. (117) must be zero for all values of l. Therefore, the solutions of (117) must be Then, we can write |ψ (0) as where we have defined At this stage, the 2l 0 + 1 coefficients ψ (0) l0m in (126a) are still undetermined. However, irrespective of their values, we always havê because from (121) and (123), it follows that This equation can be interpreted as an eigenvector equation with eigenvalue equal to 0. A direct calculation actually shows that |α E 0 ∈ E 4 belongs to the biorthogonal and where λ 0 = 0 and λ 1 = z E +1. Note that throughout this paper, we use dotless letters ı and  , as indices running from 0 to 3, while ordinary letters i and j are indices running from 1 to 4. The left eigenvectors where If, additionally, we define and we can build a complete and biorthogonal set of bases for E 4 , denoted by |α E ı , α E ı | [32]. It is a simple exercise to verify that these basis vectors satisfy the standard normalization condition for bi-orthogonal vectors, and that they form a complete basis for E 4 , where I 4 is the 4×4 identity matrix. Such a biorthogonal basis will be very useful for the next steps in perturbation theory.

B. TM waves
In this case we choose a pair of values l 0 , x We proceed as for the TE case and we write again (118) as with Next, (138) turns into an identity for where ψ Then, we can write |ψ (0) as where |ϕ (0) is defined by (126a) and By definition, from (139) and (143), it follows that Now the complete and biorthogonal set of bases for E 4 is |α M ı , α M ı | and it is defined by and where and with λ 0 = 0 and λ 1 = z M − 1. By definition,

A. Some preparatory remarks
In this section we will focus on the degenerate subspace of dimension N l0 = 2l 0 +1, denoted by D 0 ⊆ E , generated by a solution x and the vectors (140).

B. Solving the equations
We consider now the change of the degenerate vectors |ψ (0) µ0 , (µ = 1, · · · , N l0 ), when the sphere is deformed. Proceeding as in Sec. IV B, we write This set of equations (165) must hold for all µ = 1, · · · , N l0 . Note that for each value of µ the corrections x (n) µ (ε) might be different, because the departure from the spherical shape typically removes the degeneracy.
As it is customary in quantum perturbation theory, we normalize the vector |ψ µ (ε) as As in quantum mechanics, this normalization is particularly convenient for the later developments of the theory, and does not affect any physical quantity. Equation (166) implies that |ψ (n) µ for n ≥ 1, has no component along ψ Note, however, that |ψ (n) µ may have components along From (97) it follows that the perturbed vector |ψ µ0 (ε) must satisfyM (ε)|ψ µ (ε) = 0.
Substituting (165) into this equation and proceeding as in Sec. IV B, we obtain at first order in ε, where (108) has been used. To solve this equation we must project it on the three orthogonal subspaces D 0 , D I , and C .

Projecting along D0
Multiplying (169) from the left by ψ (0) ν0 | and recalling (110), we obtain where (160b) has been used to cancel the leftmost term in (169). This equation implies that which can be suggestively rewritten as The denominator α 0 |D l0 |α 0 is just a number, as shown in Appendix A. Conversely, the numerator α 0 |V (1) |α 0 is an operator in E ∞ . However, as it is sandwiched between ϕ (0) ν | and |ϕ (0) µ which are in D, we can equivalently rewrite (172) as where is the projector onto the subspace D. By definition, Written in this form, Eq. (173) tells us that the biorthogonal set |ϕ (0) µ , ϕ (0) µ | must be chosen to make the N l0 × N l0 matrix α 0 |P DV (1) P D |α 0 diagonal in the subspace D. From now on, we assume that the set of vectors |ϕ (0) µ , ϕ (0) µ | has been chosen in this way (see Appendix B for more details).
Finally, setting ν = µ, we get the value of the firstorder correction to the resonance: The next steps will determine the components of |ψ (1) µ , in D I and C . We recall that for n ≥ 1, where where, here and hereafter, Multiplying (169) from the left by ψ (0) νı |, with ı = 1, 2, 3, and recalling (110), we obtain where Substituting (181) into (180), we obtain where x (1) µ is given by (176). For given values of µ and ν, this is a set of three linear equations in the variables ψ We write it concisely as where we have defined µ0 − x (1) µ δ νµ α ı |D l0 |α 0 , (184b) (ı,  = 1, 2, 3), and the 3 × 3 matrix D with elements A straightforward calculation gives for TE waves, and for TM waves, where z E and z M are given by (119) and (137), respectively. These matrices are invertible, because where f E l (x) and f M l (x) are given by (27) and (28), respectively, and we have denoted by x respectively. Using these equations, we can eventually write the solution of (183), as for TE waves and for TM waves.

Projecting along C
Multiplying (169) from the left by l, m, i|, with l = l 0 and i = 1, 2, 3, 4, we obtain From (110) and l = l 0 , it follows that the last term on the right-hand side is identically zero. Substituting (178c) into (193) we find, after a little calculation, which can be recast into where now we have defined X and Y as Finally, a straightforward calculation gives, for l = l 0 , l, m, 1|ψ (1)

C. Summary of the first-order perturbation theory
In this section we have determined the first-order corrections x (1) µ along the degenerate subspace D 0 . For this, we need to solve the second-order equation (106c). This will be done in the next section.

VII. SECOND-ORDER EQUATIONS I. NON-DEGENERATE CASE
In this section we are going to solve (106c), when the degeneracy is lifted to first order, that is when and µ, ν = 1, · · · , N l0 . To begin with, we use (108) to rewrite (106c) aŝ Next, we proceed as in first-order theory, projecting this equation on the subspaces D 0 to find x (2) µ and ψ (0) ν0 |ψ (1) µ . We remind that, at second order we are interested only in the resonance corrections, so we do not need to determine the full vector |ψ (2) µ .

Projecting along D0
Multiplying (200) from the left by ψ (0) ν0 | and recalling (110), we obtain Now, for clarity, we calculate the five addends on the right-hand side of (201) separately.
Fourth addend [from (201d)]: Fifth addend [from (201e)]: where (161) has been used. Now that we have all the terms, we can use (201) evaluated for ν = µ to obtain x (2) µ . Next we set ν = µ to get ψ (0) ν0 |ψ (1) µ . In the first case we find, In the second case we obtain with ν = µ. It is instructive to rewrite (207) in a compact form using the first-order equation (169) in the bra form Multiplying this equation from the right by |ψ (1) µ , we obtain (210) Substituting this result in (201) we find, after a simple manipulation, This expression is only formal in the sense that it contains unknown coefficients. However, it makes clear what the complicated equation (207) actually means. This completes the calculation for the second-order corrections, when the degeneracy is lifted to first order. If this is not the case, we need a different procedure, which will be developed in the next section.
VIII. SECOND-ORDER EQUATIONS II. DEGENERATE CASE

A. Some preparatory remarks
Now we consider the case when the degeneracy is only partially removed to first order. Without loss of generality, we assume that the first N first-order corrections x (1) µ are equal to each other: where 1 < N ≤ N l0 , and (171) is still valid. Consequently, the initial degenerate subspace D 0 breaks in two parts, denoted by D 0N and D 0M , with N + M = N l0 , and defined by As before, because of the remaining N -fold degeneracy, it is convenient to define a new orthonormal basis in D 0N , denoted by |ψ (0) A0 and defined by where A = 1, · · · , N . The coefficients ϕ (0) Aµ are to be determined. As it should be customary now, we think of this basis as a part of the biorthogonal set |ψ with We build up the perturbation theory as usual, Note that in (217b) the first-order correction x (1) has no label, according to (212). Using the fundamental equation (99),M we can obtain the familiar chain of equationŝ etc., where (108a) has been used. Next, we need to adapt to the present case the normalization condition (167). As shown in Appendix C, the new condition is and n ≥ 1. Therefore, similarly to (177) and (178), we can write now where The zeroth-order equation (219a) is trivially satisfied. The first-order equation (219b), also becomes an identity when projected on D 0N and D 0M . However, it furnishes the coefficients ψ (0) νı |ψ νı |, with ν = 1, · · · , N l0 and ı = 1, 2, 3, we obtain where (108b) has been used. The first term of this sum is because diagonal operatorsD cannot connect D I with C , and (161) has been used. The second term does not require calculations, being simply where the coefficients ϕ (0) Aµ are still to be determined and (214) has been used. Finally, the third and last term is where (161) has been again used. Substituting (224)-(226) into (223), we obtain for ν = 1, . . . , N l0 and ı = 1, 2, 3. For each value of ν, Eq.
(227) can be written as a matrix equation of the form where the 3 × 3 matrix D is defined by (185), that is, l0 |α  , and now Specifically, D is given by (186) for TE waves, and by (187) for TM waves. Therefore, we know that it is in-vertible and we can formally write where D −1 is the matrix inverse of D. Using the standard notation D −1 ({ı, }), (ı,  = 1, 2, 3), to denote the principal submatrix of D −1 that lies between row ı and row , and between column ı and column  [33], we can write for TE waves, and for TM waves, where z E and z M are given by (119) and (137), respectively. Thus, the expression between curly brackets in (230) is completely determined. In the remainder we will indicate it compactly with to rewrite (230) as l |j . Therefore, we do not need to make additional calculations and we can write directly l, m, i|ψ where we have defined In this expression where where x (0) = x σ l0n , with σ = E, M .
C. Second order equations

Projecting along D0N
Let us set ν ≤ N . Multiplying (219c) from the left by ψ (0) ν0 | and using (108b), we obtain The left-hand side of this equation vanishes because of (160b). The four addends on the right-hand side of (243) are calculated as follows.
First addend [from (243a)]: where (170) has been used. Note that in the first line of (244), the coefficients ψ are unknown. However, we will see soon that such term is canceled by an analogous one in the second addend (243b).
Second addend [from (243b)]: As anticipated, the term in this expression, cancels with the same term in (244).
Third addend [from (243c)]: Fourth addend [from (243d)]: Summing all these addends, after straightforward manipulation we eventually obtain N ν=1 where we have defined the N × N matrix M (2) , by the elements Equation (248) is an eigenvalue equation that gives us both the second-order resonance corrections x (2) A , as the eigenvalues of M (2) , and the basis vectors |ψ (0) A0 , as the associated eigenvectors. This completes our calculations.

IX. OBLATE SPHEROID
In this section we apply our theory to nearly spherical dielectric resonators, which are rotationally invariant around the z axis, with h(θ, φ) = h(θ). This permits us to illustrate the use of degenerate perturbation theory in the case in which the degeneration is only partially removed to first order. The unperturbed system is, as always in this work, a dielectric sphere of radius a and refractive index n 1 , surrounded by vacuum or air with n 2 = 1. For practical reasons (the numerical results are more accurate), we will choose n 1 = 2.
As a specific example, we consider as nearly spherical resonator, an oblate spheroid with semiaxes (a(1 + δ), a(1 + δ), a), where 0 < δ 1 quantifies the magnitude of the deformation. The equation of the spheroid in spherical coordinates is r = a where the Taylor expansion truncated at first order, is a good approximation for 0 < δ 1. Thus, in the remainder we will set as deformation function and for the equation defining the approximate oblate spheroid, for both the perturbative and the numerical calculations. To perform the latter, we used the COM-SOL Multiphysics software (Wave Optics Module) [34]. For illustration purposes, we choose as unperturbed resonances (30) and (31) x 7.248 − i 4.325 × 10 −3 , with l 0 = 10 and n 0 = 1, for TE and TM waves, respectively. We take the magnitude of the deformation to be equal to δ = 0.01 and δ = 0.05 for the TE and TM waves, respectively. The choice of n 0 = 1 (first radial mode), is suggested by the fact that higher-order radial numbers (n > 1) mark lossy waves not localized near the surface of the resonator, which are of low practical interest [1].  circles are obtained by direct numerical simulations, and the blue closed circles by solving the eigenvalue equation (248). Deforming the sphere into a spheroid partially lifts the degeneracy, thus yielding l 0 + 1 = 11 distinct resonances, each characterized by a different value of |m| = 0, 1, · · · , l 0 . The remaining twofold degeneracy is due the rotational invariance of the spheroid with respect to the z axis, which implies that the physics is the same for clockwise (m > 0) and counterclockwise (m < 0) waves. Note that waves with |m| < l 0 have a polar angle θ extension, growing with l 0 − |m|. This implies that they are more sensitive to surface deformations.
where δx num has been estimated as the absolute error between the theoretical (exact, |m|-independent) values obtained by solving (29) with l = 10, n = 1, and the (|m|-dependent) numerical results for a perfectly spherical cavity. To give a quantitative estimate of the error, we have also plotted the average relative error (orange dashed lines). Figures 5 and 6 are the same as figs. 3 and 4, respectively, but with δ = 0.05. Overall, all plots exemplify the goodness of the second-order perturbation theory we have developed, even for non-equatorial modes with |m| < l 0 . To produce the plots above, we greatly benefited from the fact that the infinite sums that appears in (207)-(208) and (249), actually contain only a finite number of terms according to the rule where L max is determined by the expansion of the deformation function f (θ, φ) in terms of the spherical harmonics: For example, from (251) it follows that so that L max = 2. The rule (254) is empirically determined. However, it could be rigorously proven by writing the product of spherical harmonics in terms of the Wigner 3jsymbols [35]. Such products appear in the quantities X l m , f (θ, φ)Ψ lm and X l m Y lm (θ, φ)e (θ, φ) , used in (A7b) and (A8b). A practical example of the use of the 3j-symbols in this kind of calculation, can be found in Sec. VI C of Ref. [25], Eqs. (117)-(122).

X. SUMMARY
We have developed a boundary conditions perturbation theory to determine the electromagnetic resonances of nearly spherical dielectric resonators. The threedimensional nature of the resonator and the vector character of the electromagnetic field, dictated the use of vector spherical harmonics for handling the problem, as opposed to the more familiar scalar spherical harmonics, the latter being typically employed in problems with spherical or nearly-spherical symmetry. By imposing standard electromagnetic boundary conditions at the surface of the resonator separating two different dielectric media, we obtained an exact algebraic homogeneous system of linear equations. The mathematical correspondence between linear operators and matrices, allowed us to reformulate the problem in the language of quantum mechanics, and to use the well-known Rayleigh-Schrödinger perturbation theory, to build up a perturbation series for the resonances of the electromagnetic field, up to and including second-order terms. However, as dielectric resonators are de facto open systems, we had to use the mathematical machinery of non-Hermitian operators and biorthogonal bases. We considered both simple and degenerate unperturbed spectra, including the case when degeneracy is not fully removed to first order. For the latter instance, exemplified by the spectrum of an oblate spheroid resonator, we have compared the predictions of our theory with numerical calculations, finding excellent agreement.
The main results are represented by Eqs. (176), (207), (248), and (249). These formulas can be used to calculate the spectrum of the electromagnetic resonances of arbitrarily deformed nearly-spherical dielectric resonators of any size, provided the conditions (50) for the applicability of the perturbation theory are satisfied. Notably, as second-order terms are included, this theory can also be used for the calculation of the spectra of spherical resonators with random surface roughness. This is the case, for example, of helium droplets with thermally excited capillary waves [24]. functions (16) isolating their common denominator b αl (k α a), that is Note that in each of these expressions the dependence on ε enters in different ways in the numerator and the denominator, because with W = Ψ, Φ, Y .