Stéphane Balac

Patrice Féron

Whispering gallery modes volume computation in optical micro-spheres

This report is devoted to the computation of the volume of whispering gallery modes in optical micro-spheres.

We first derive the mathematical expression of TE and TM modes in a homogeneous dielectric optical micro-sphere from the general set of Maxwell's equations. Then we present a method to numerically compute the volume of whispering gallery TE or TM modes which do not require any assumptions on the mode numbers as it is usually the case. Last we present the Matlab Toolbox WGMode written from the present study designed to explore resonance conditions, to visualize whispering gallery TE or TM modes in a micro-sphere and to compute the volume of any given mode.

Introduction

Whispering gallery modes (WGM) are specific resonances of a wave field inside a given cavity with smooth edges. They correspond to waves circling around the cavity, supported by continuous total internal reflection of the cavity surface, that meet the following resonance condition: after one round-trip they return to the same point with the same phase and hence interfere constructively with themselves, forming standing waves, see [START_REF] Righini | Whispering gallery mode microresonators: Fundamentals and applications[END_REF] for details. WGM's have attracted much attention due to potential applications in photonics, quantum electrodynamics, telecommunication, etc. Such applications in technological and scientific fields are e.g. the realization of micro-lasers, narrow filters, optical switching, ultrafine sensing, displacement measurements, high resolution spectroscopy, Raman sources and studies of nonlinear optical effects, see [START_REF] Matsko | Review of applications of whispering-gallery mode resonators in photonics and nonlinear optics[END_REF] for a review.

The quantum and nonlinear properties of an optical microcavitiy are well described by its quality-factor Q and its mode volume V. The spontaneous emission rate enhancement of a two-level system embedded in a cavity is given by the ratio between the quality-factor and the mode volume [START_REF] Purcell | Proceedings of the American Physical Society[END_REF][START_REF] Gérard | Enhanced spontaneous emission by quantum boxes in a monolithic optical microcavity[END_REF] through the Purcell factor FP :

FP = 3 4π 2 Q V λ N0 3 (1.1)
where λ is the resonant wavelength of the cavity and N0 its refractive index. The values of Q and V also well characterize nonlinear microcavities [START_REF] Bravo-Abad | Enhanced nonlinear optics in photonic-crystal microcavities[END_REF][START_REF] Wu | Quantum theory for microcavity enhancement of second harmonic generation[END_REF]. For example, the power bistability threshold P th of a lossless cavity including a Kerr material (with a nonlinear refractive index N2) is given in a first approximation by:

P th = 2πN 2 0 3 √ 3N2λ V Q 2 . (1.2)
In this report, we give a numerical method to calculate the mode volume V for a dielectric spherical WGM microresonator, defined as the integral over the whole space of the energy density normalized by its maximum value inside the micro-sphere, i.e.

V = 1 wmax R 3 w(x1, x2, x3) dx1dx2dx3 (1.3)
where w denotes the energy density given as a function of the position vector x = (x1, x2, x3) by

w(x) = 1 2 ε(x) 2 E(x) E * (x) + 1 2µ0
B(x) B * (x) (1.4) and wmax denotes the maximum value of the energy density inside the micro-sphere and the quantities E * and B * stand respectively for the adjoint (conjugate transpose) of the electric field E and magnetic induction B.

The document is organized as follows. In section 2 we introduce the general framework for this problem as set by Maxwell's equations. The problem turns out to solving a vectorial wave equation in spherical coordinates. The vectorial wave equation can be replaced by Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr a simultaneous system of three scalar equations, but the solution of this system for any component in coordinates systems other than the rectangular coordinates system is most cases impractical because the three components of the unknown vector field are connected in each equation. Thus in section 3 we use Hansen method to obtain the solutions of the vector wave equation in spherical coordinates from the solutions of the corresponding scalar wave equation. We deduce the expression of TE and TM modes in a spherical microresonator in terms of vector spherical harmonics functions and the corresponding modal equations that determine the resonance conditions. The expressions of TE and TM modes are used in section 4 to compute analytically the TE and TM mode volumes for a dielectric spherical WGM microresonator. We have also written a Matlab Toolbox named WGMode dedicated to the study of Whispering gallery modes in spherical microresonator and based on the results shown in this document. The scripts of the WGMode Matlab Toolbox allow to explore resonance conditions for TE or TM modes, to visualize whispering gallery TE or TM modes in a micro-sphere and to compute the volume of any given mode. The features of the WGMode Matlab Toolbox are presented in section 5.

General framework

Maxwell's equations

The behavior of the electromagnetic field at time t and position x in an optical device such as a micro-resonator can be described from the general set of Maxwell's equations:

           curl E(x, t) + ∂ ∂t B(x, t) = 0 div B(x, t) = 0 div D(x, t) = ρ(x, t) curl H(x, t) - ∂ ∂t D(x, t) = j(x, t) (2.1)
where B is the magnetic flux (in units of Tesla, T), H is the magnetic field strength (in units of A/m), D is the electric displacement field (in units of C/m 2 ), E is the electric field (in units of V/m), ρ is the electric charge density (in units of C/m 3 ) and j is the current density (the sum of source currents density js and induced currents density jc in units A/m 2 ). These equations have to be considered in the sense of distributions in R 3 , see e.g. [START_REF] Pinchard | Electromagnétisme Classique et Théorie des Distributions[END_REF][START_REF] Struillou | Analyse de Fourier -Théorie et Applications pour l'Iingénieur et le Physicien[END_REF].

The magnetic flux B and the magnetic field strength H are connected through the relation B(x, t) = µ0 (H(x, t) + M(x, t))

(2.2)
where M is the magnetization vector (in units of A/m) and µ0 denotes the magnetic permeability in vacuum (µ0 = 4π 10 -7 H/m). Similarly, the electric displacement field D and the electric field E are connected through the relation

D(x, t) = ε0E(x, t) + P(x, t) (2.3)
where P is the electric polarisation vector (in units of C/m 2 ) and ε0 denotes the permittivity in vacuum (ε0 = 1 36π 10 -9 F/m). In addition to Maxwell's equations there are constitutive relationships to describe media's properties and effects.

Constitutive relationships for linear materials

When the media under consideration can be assumed to be isotropic and linear with respect to their behavior to the electromagnetic field, the constitutive relationships are in the form

   M = χm H P = ε0χeE jc = σ E (2.4)
Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr where χm and χe are dimensionless quantities known respectively as the magnetic susceptibility and the electric susceptibility and σ is the electric conductivity. The last relation is Ohm's law in a microscopic form. In this study we will assume that all the media under consideration are free of charge (i.e. we assume that ρ = 0) and non conductive media (i.e. we assume that σ = 0). Taking into account (2.4), relations (2.2) and (2.3) can be expressed as

B = µ H D = ε E (2.5)
where the dielectric permittivity ε = ε0(1 + χe) and the magnetic permeability µ = µ0χm are characteristic constants of linear isotropic media.

The three parameters (µ, σ, ε) describe the electromagnetic properties of a material. It is common to write the magnetic permeability µ and the electric permittivity ε as µ = µ0 µr and ε = ε0 εr where µr and εr denote respectively the relative permeability and relative permittivity (also called the dielectric constant). They are dimensionless quantities. For most of the media involved in optical devices, it can be assumed that µr = 1. The speed of light, that is to say the electromagnetic wave velocity in vacuum, is given by

c = 1 √ µ0ε0
and it is usual in optics to set εr = N 2 where N is referred as the optical index of the dielectric medium. For an electromagnetic wave with wavelength λ and pulsation ω = 2πc/λ, we denote by k0 = 2π λ = ω c its wavenumber in vacuum and by k = N k0 its wavenumber in a medium characterized by an optical index N . If we assume that the involved media are linear, isotropic and non dispersive then

• the magnetic induction B is proportional to the magnetic field H and can be expressed as B(x, t) = µ(x) H(x, t) where µ is a real valued scalar function of the space variable;

• the electric displacement vector D can be expressed as D(x, t) = ε(x) E(x, t) where ε is a real valued scalar function of the space variable.

Using the above constitutive relations, it is possible to write Maxwell's equations (2.1) in terms of electric field E and magnetic induction B alone as:

             curl E(x, t) + ∂ ∂t B(x, t) = 0 div B(x, t) = 0 div (ε(x) E(x, t)) = 0 curl 1 µ(x) B(x, t) -ε(x) ∂ ∂t E(x, t) = 0 (2.6)

Continuity conditions for Maxwell equation

The electromagnetic field in an isotropic and linear homogeneous medium are "regular" and partial differential operators in Maxwell equations (2.1) can be considered in the usual sense of differential calculus, see e.g. [START_REF] Jackson | Classical Electrodynamics Third Edition[END_REF][START_REF] Stratton | Electromagnetic Theory[END_REF]. In such a case, the electromagnetic fields on both sides of an interface between two media with different electromagnetic characteristics are interrelated through the following Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr conditions:

E ∧ n = 0 (2.7a) B • n = 0 (2.7b) ε E • n = ρs (2.7c) 1 µ B ∧ n = -jΣ (2.7d)
where jΣ is the surface electrical current density at the interface and ρs is the surface charge density, n denotes the unit vector normal to the interface between the two different involved media and the brackets [ ] stand to denote the jump of the quantity inside the brackets across the boundary.

Optical micro-spheres

Optical micro-spheres, see Fig. 1, are high quality factor optical resonators used for active and passive photonic applications such as microlaser cavities and filters. The radius of the micro-spheres considered typically ranges from around ten up to several hundred of micrometers. In the present study, it is assumed that the micro-sphere is made of a homogeneous isotropic linear medium. Namely, throughout the document, we assume that the micro-sphere is made of a homogeneous dielectric material where µ = µ0 and ε = N 2 ε0 with the optical index N being a positive real number. Typical values for N is around 1.45 for silicate glass. Moreover we assume that the micro-sphere is uncharged (ρ = 0), non conductive (σ = 0) and not the place of electrical currents js = 0. The wavelength range of interest is 800 to 1600 nano-meters. where ω is the angular frequency (in units of radians per second), E0 and B0 are the field amplitude and φ0 is the phase angle (in units of radians).

Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

It is very convenient in practice to represent sinusoidal fields in complex notation thanks to Euler's formula. A phasor representation of the electric field E(x, t) (the same holds for the magnetic induction B) is the complex number E(x) with a magnitude E0(x) and a phase φ0 such that E(x) = E0(x) exp(iφ0(x)).

If we multiply E(x) by exp(iωt) and apply Euler's formula,

E(x) exp(iωt) = E0(x) exp(i(ωt + φ0(x)) = E0(x) cos(i(ωt + φ0(x)) + E0(x) sin(i(ωt + φ0(x)).
The real part of E(x) exp(iωt) is the electric field E(x, t)

E(x, t) = Re E(x) exp(iωt) . (2.8)
In the sequel, we will use the following phasor representations

E(x, t) = Re E(x) exp(iωt) where E(x) = E0(x) exp(iφ0(x)) and B(x, t) = Re B(x) exp(iωt) where B(x) = B0(x) exp(iφ0(x))
2.3.2. Maxwell's equations for sinusoidal electromagnetic waves From equations (2.6) we deduce that the complex fields E and B satisfy the following set of time independent equations :

curl E(x) + iω B(x) = 0 (2.9a) div B(x) = 0 (2.9b) div (ε(x) E(x)) = 0 (2.9c) curl 1 µ(x) B(x) -i ωε(x) E(x) = 0 (2.9d)
In the case of a homogeneous dielectric domain characterized by a constant electric permittivity ε and a magnetic permeability µ0, taking the curl of the curl equations in (2.6) we obtain:

curl curl E = - ∂ ∂t curl B = -µ0ε ∂ 2 E ∂t 2 and curl curl B = µ0ε ∂ ∂t curl E = -µ0ε ∂ 2 B ∂t 2 .
Then by using the vector identity, see e.g. [START_REF] Stratton | Electromagnetic Theory[END_REF],

curl (curl V) = ∇ (div V) -∆V
where V is any vector function of the space variables, and taking into account that the electromagnetic field is divergence free in a homogeneous dielectric medium, see equations (2.6), we find that the electromagnetic field satisfies the wave equations: If we introduce the wave number k such that k 2 = ω 2 εµ0 then the complex vector fields E and B satisfies a system of equations in the form

∂ 2 ∂t 2 E(x, t) -µε ∆E(x, t) = 0 and ∂ 2 ∂t 2 B(x, t) -µε ∆B(x, t) =
∆C(x) + k 2 C(x) = 0 (2.12) div C(x) = 0 (2.13)

Whispering Gallery Modes in an isotropic dielectric micro-sphere

An isotropic dielectric micro-sphere can retain the light that has been injected into it from a coupling waveguide (most of the time a micro-taper with strong evanescent field is required because far field electromagnetic energy could not be captured by a dielectric microsphere). Whispering Gallery Modes (WGMs) are particular form of the electromagnetic field inside an axisymmetic dielectric waveguide such as micro-sphere when some resonance conditions are satisfied. Actually, in an isotropic optical micro-sphere, the "dielectric wall" (corresponding to the difference of dielectric properties inside and outside the micro-sphere) prevents the trapped light from being scattered easily. When resonance condition is satisfied, the electromagnetic field of WGMs can be reinforced by coupling in light and forming spatial coherent modes. A dielectric micro-sphere can accommodate many different discrete resonance frequencies with respect to its dimension and materials. From an experimental point of view, the modes properties of the micro-sphere can be evaluated by measuring and analyzing the coupling characteristics and resulting interference patterns from a micro-taper output. The WGMs are excited either by evanescent field or by the modes coupling between the WGMs.

We introduce the spherical coordinate system as shown in Fig. 2 where

• the radius or radial distance r is the Euclidean distance between the origin O set at the microsphere center and the current point P ;

• the inclination (or polar angle) θ is the angle between the zenith direction and the position vector OP measured between 0 and π radian ;

• the azimuth (or azimuthal angle) ϕ is the signed angle measured between 0 and 2π radian from the azimuth reference direction to the orthogonal projection of the position vector OP on the reference plane.

We denote by er, e θ , eϕ the spherical unit vectors and by r = r er = OP the radial vector where r = OP 2 is the radial distance.

The analytic solutions of WGMs can be obtained by solving vector Helmhotz equations (2.12) for the electromagnetic field (E, B) in the domain occupied by the dielectric sphere and in the outer domain, together with appropriate boundary conditions at the interface deduced from (2.7a)-(2.7d).

The vectorial wave equation (2.12) can be replaced by a simultaneous system of 3 scalar equations, but the solution of this system for any component in coordinates systems other than the rectangular coordinates system is most cases impractical because the 3 components of the unknown vector field C are connected in each equation. For instance, the vectorial Laplacian in spherical coordinates reads

∆C =         1 r ∂ 2 (rCr) ∂r 2 + 1 r 2 ∂ 2 Cr ∂ϕ 2 + 1 r 2 sin 2 (ϕ) ∂ 2 Cr ∂θ 2 + cotan(ϕ) r 2 ∂Cr ∂ϕ - 2 r 2 ∂Cϕ ∂ϕ - 2 r 2 sin(ϕ) ∂C θ ∂θ - 2Cr r 2 - 2 cotan(ϕ) r 2 Cϕ 1 r ∂ 2 (rCϕ) ∂r 2 + 1 r 2 ∂ 2 Cϕ ∂ϕ 2 + 1 r 2 sin 2 (ϕ) ∂ 2 Cϕ ∂θ 2 + cotan(ϕ) r 2 ∂Cϕ ∂ϕ - 2 r 2 cotan(ϕ) sin(ϕ) ∂C θ ∂θ + 2 r 2 ∂Cr ∂ϕ - 1 r 2 sin 2 (ϕ) Cϕ 1 r ∂ 2 (rC θ ) ∂r 2 + 1 r 2 ∂ 2 C θ ∂ϕ 2 + 1 r 2 sin 2 (ϕ) ∂ 2 C θ ∂θ 2 + cotan(ϕ) r 2 ∂C θ ∂ϕ + 2 r 2 sin(ϕ) ∂Cr ∂θ + 2 cotan(ϕ) r 2 sin(ϕ) ∂Cϕ ∂θ - 1 r 2 sin 2 (ϕ) C θ        
Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr However, when dealing with the spherical coordinates system the method of Hansen can be used for the determination of 3 independent vector solutions of (2.12). The method of Hansen is for instance detailed in [START_REF] Stratton | Electromagnetic Theory[END_REF]. It permits to deduce solutions of the vector wave equation in spherical coordinates directly from the solutions of the corresponding scalar wave equation.

In the next section we focus on the resolution of the scalar wave equation (i.e. the component-wise vector wave equation (2.12)). In section 3.2 we will consider the way the vectorial solution to the vector wave equation (2.12) can be deduced from the solution to the scalar wave equation by the method of Hansen.

The scalar wave equation in spherical coordinates

The scalar wave equation reads

∆Ψ(x) + k 2 Ψ(x) = 0. ( 3.1)
In spherical coordinates, equation (3.1) is separable. We are looking for solutions in the form

Ψ(r, θ, ϕ) = f1(r) × f2(θ) × f3(ϕ) (3.2)
where f1, f2 and f3 denote 3 complex valued functions of the real variable. The angles ϕ and θ are the azimuthal angle and the polar angle respectively, see Fig. 2. Relative to the spherical coordinates system, the Laplace operator has the form [START_REF] Stratton | Electromagnetic Theory[END_REF]

∆Ψ = 1 r 2 ∂ ∂r r 2 ∂ ∂r Ψ + 1 r 2 1 sin θ ∂ ∂θ sin(θ) ∂ ∂θ Ψ + 1 sin 2 θ ∂ 2 ∂ϕ 2 Ψ . (3.3)
In spherical coordinates, the partial differential equation (3.1) lends itself to being separated into a system of 3 ordinary differential equations: upon substituting ansatz (3.2) in the wave equation (3.1) one finds that the 3 unknown functions f1, f2 and f3 satisfy

r 2 f 1 (r) + 2r f 1 (r) + (k 2 r 2 -p 2 )f1(r) = 0 (3.4a) 1 sin(θ) d dθ (sin(θ) f 2 (θ)) + (p 2 - q 2 sin 2 (θ) ) f2(θ) = 0 (3.4b) f 3 (ϕ) + q 2 f3(ϕ) = 0 (3.4c)
Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

where the real parameters p and q are the separation constants.

Angular dependency

Since the optical properties of the micro-sphere are independent of the azimuthal angle ϕ, it is necessary for f3 to be a periodic function with period 2π. Then, the solution to equation (3.4c) reads

f3(ϕ) = C1 e imϕ + C2 e -imϕ (3.5)
where the separation constant q coincides with an integer m ∈ Z and C1 and C2 denote two complex constants. The change of variable η = cos(θ) in equation (3.4b) leads to the equation

(1 -η 2 ) g 2 (η) -2η g 2 (η) + (p 2 - m 2 1 -η 2 ) g2(η) = 0 (3.6)
where the new unknown function g2 is defined by the relations f2(θ) = g2(cos(θ)). The solutions to the linear second order ordinary differential equation (ODE) (3.6) are the socalled hypergeometric functions, see [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF]

g2(η) = P m (η) := (-1) m (1 -η 2 ) m 2 d m dη m P (η)
where P denotes the Legendre polynomial of degree defined for all x ∈] -1, 1[ by

P (x) = 1 2 ! d dx (x 2 -1) .
The associated Legrendre functions are sometimes defined without the multiplicative constant (-1) m which is known as the Cordon-Shortley phase factor. Following the standard convention, we set for negative values of the order

P -m = (-1) m ( -m)! ( + m)! P m .
For a fixed integer m, the associated Legrendre functions satisfy the following orthogonality conditions for all k, ∈ N such that 0 m <

P m P m dx = 2( + m)! (2 + 1)( -m)! .
Finally we find the functions product f2(θ) f3(ϕ) in the form of

Y m (θ, ϕ) = C m P m (cos(θ)) e imϕ - m , ∈ N (3.7)
where Y m is known as the Spherical Surface Harmonics of degree and order m, see [START_REF] Atkinson | Spherical Harmonics and Approximations on the Unit Sphere: An Introduction[END_REF]. The normalization constant C m is taken to be

C m = (2 + 1) 4π ( -m)! ( + m)! . (3.8)
Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

With this normalization convention, the complex conjugate of the Spherical Surface Harmonics of degree and order m is

Y m (θ, ϕ) = (-1) m Y -m (θ, ϕ).
The Spherical Surface Harmonics form a complete set of orthonormal functions and thus they form an orthonormal basis of the Hilbert space of square-integrable functions on the unit sphere [START_REF] Atkinson | Spherical Harmonics and Approximations on the Unit Sphere: An Introduction[END_REF][START_REF] Struillou | Analyse de Fourier -Théorie et Applications pour l'Iingénieur et le Physicien[END_REF]. Namely, on the unit sphere any square-integrable function can be expanded as a linear combination of Spherical Surface Harmonics as:

g(θ, ϕ) = +∞ =0 m=- g m Y m (θ, ϕ) (3.9)
where the equality holds in L 2 ([0, π] × [0, 2π], C) and

g m = 2π 0 π 0 g(θ, ϕ) Y m (θ, ϕ) sin θ dθ dϕ.

Radial dependency

There remains to achieve the identification of the radial function f1. It satisfies the following Spherical Bessel equation deduced from (3.4a)

r 2 f 1 (r) + 2rf 1 (r) + (k 2 r 2 -( + 1))f1(r) = 0. (3.10)
By the change of variable x = kr, equation (3.10) for the new unknown function g1 defined by g1(x) = f1(x/k) reads

x 2 g 1 (x) + 2x g 1 (x) + (x 2 -( + 1)) g1(x) = 0. ( 3.11)
The two linearly independent solutions to the linear second order ordinary differential equation (3.11) are known as the spherical Bessel functions of the first and second kinds [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF][START_REF] Gradstein | Table of Integrals, Series, and Products. Table of Integrals, Series, and Products Series[END_REF][START_REF] Temme | Special Functions: An Introduction to the Classical Functions of Mathematical Physics[END_REF] denoted j and y . They are related to the ordinary Bessel functions J and Y for x > 0 by the relations

j (x) = π 2x J + 1 2 (x), (3.12
)

y (x) = π 2x Y + 1 2 (x) = (-1) +1 π 2x J --1 2 (x). (3.13)
When dealing with the problem of scattering of electromagnetic waves by a sphere, it is convenient to introduce the Riccati-Bessel functions of the first and second kinds denoted ψ and χ . They are related to the Spherical Bessel functions and to Bessel functions for x > 0 by the relations

ψ (x) = x j (x) = πx 2 J + 1 2 (x), (3.14)
χ (x) = -x y (x) = - πx 2 Y + 1 2 (x). (3.15)
We deduce that the radial function f1 can be expressed as a linear combination of the Spherical Bessel functions as

f1(r) = α j (kr) + β y (kr)
where α and β denote two complex constant numbers. In terms of the Riccati-Bessel functions, the unknown function f1 can be expressed as

f1(r) = α ψ (kr) kr -β χ (kr) kr .
Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

We conclude that the general solution of the scalar wave equation (3.1) in spherical coordinates, expressed either in terms of the Spherical Bessel functions or in terms of the Ricatti-Bessel functions, reads

Ψ(r, θ, ϕ) = +∞ =0 m=- (α j (k r) + β y (k r)) Y m (θ, ϕ) = +∞ =0 m=- α ψ (k r) kr -β χ (k r) kr Y m (θ, ϕ).
(3.16)

Hansen method for solving the vector wave equation in spherical coordinates

We have shown in the previous section that the characteristic solution of the scalar wave equation (3.1) in spherical coordinates is

Ψ m (r, θ, ϕ) = z (k r) Y m (θ, ϕ)
where z denotes one of the two linearly independent solutions (j or y ) to the spherical Bessel equation (3.10). We are now concerned by solving the vector wave equation

∇(div C) -curl curl C + k 2 C = 0 (3.17)
in a domain which is either the inner or the outer of a ball. Note that for solenoidal (i.e. divergence free) fields, the vectorial wave equation (3.17) can be simplified into

∆C + k 2 C = 0. (3.18)
In spherical coordinates, the 3 independent vector solutions of (3.17) can be constructed form the characteristic solution Ψ m of the scalar wave equation (3.1) by Hansen method as described in [START_REF] Stratton | Electromagnetic Theory[END_REF], chap. 7. One solution, denoted L m in the sequel, can be found by taking the gradient of Ψ m . It reads

L m (r, θ, ϕ) = ∇Ψ m (r, θ, ϕ) =        ∂ ∂r Ψ m (r, θ, ϕ) 1 r ∂ ∂θ Ψ m (r, θ, ϕ) 1 r sin(θ) ∂ ∂ϕ Ψ m (r, θ, ϕ)        =        ∂ ∂r z (k r) Y m (θ, ϕ) 1 r z (k r) ∂ ∂θ Y m (θ, ϕ) 1 r sin(θ) z (k r) ∂ ∂ϕ Y m (θ, ϕ)        . (3.19)
A second solution is found in the form of the vector

M m (r, θ, ϕ) = ∇Ψ m (r, θ, ϕ) ∧ r er = L m (r, θ, ϕ) ∧ r =      0 1 sin(θ) ∂ ∂ϕ Ψ m (r, θ, ϕ) - ∂ ∂θ Ψ m (r, θ, ϕ)      =      0 1 sin(θ) z (k r) ∂ ∂ϕ Y m (θ, ϕ) -z (k r) ∂ ∂θ Y m (θ, ϕ)      . (3.20)
Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

The third independent solution is

N m (r, θ, ϕ) = 1 k curl M m (r, θ, ϕ) = 1 k         - 1 r sin(θ) ∂ ∂θ sin(θ) ∂ ∂θ Ψ m (r, θ, ϕ) + 1 sin(θ) ∂ 2 ∂ϕ 2 Ψ m (r, θ, ϕ) - 1 r ∂ ∂r -r ∂ ∂θ Ψ m (r, θ, ϕ) 1 r ∂ ∂r r sin(θ) ∂ ∂ϕ Ψ m (r, θ, ϕ)        
.

The first component of N m can be rewritten as

(N m )1(r, θ, ϕ) = r k -∆Ψ m (r, θ, ϕ) + 1 r 2 ∂ ∂r r 2 ∂ ∂r Ψ m (r, θ, ϕ) .
From the scalar wave equation (3.1) we have ∆Ψ m = -k 2 Ψ m so that

(N m )1(r, θ, ϕ) = r k k 2 Ψ m (r, θ, ϕ) + 1 r 2 ∂ ∂r r 2 ∂ ∂r Ψ m (r, θ, ϕ) = 1 k ∂ 2 ∂r 2 rΨ m (r, θ, ϕ) + kr Ψ m (r, θ, ϕ). Moreover, since Ψ m (r, θ, ϕ) = z (r)Y m (θ, ϕ) we deduce that N m (r, θ, ϕ) =        1 k d 2 dr 2 rz (k r) + k 2 r Y m (θ, ϕ) 1 kr ∂ ∂r r z (k r) ∂ ∂θ Y m (θ, ϕ) 1 kr sin(θ) ∂ ∂r r z (k r) ∂ ∂ϕ Y m (θ, ϕ)        .
Finally, since z denotes one of the two linearly independent solutions to the spherical Bessel equation (3.10), we may rewrite N m as

N m (r, θ, ϕ) =        ( + 1) kr z (k r) Y m (θ, ϕ) 1 kr ∂ ∂r r z (k r) ∂ ∂θ Y m (θ, ϕ) 1 kr sin(θ) ∂ ∂r r z (k r) ∂ ∂ϕ Y m (θ, ϕ)        . ( 3.21)
The vector functions L m , M m and N m have several notable properties that follows directly from their definition. The vector functions L m is irrotational (its curl is zero) whereas M m and N m are solenoidal (their divergence is zero), i.e.

curl L m = 0, div M m = 0, div N m = 0. (3.22)
Moreover the 3 vector functions L m , M m and N m are related to each other through the relations

M m = L m ∧ r, N m = 1 k curl M m (3.23)
deduced from their definition and by

M m = 1 k curl N m (3.24)
deduced from the vector wave equation (3.17) for M m .

Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

In order to simplify the expression of the three independent solutions to the vector wave equation (3.17), it is convenient to introduce the Vector Spherical Harmonics defined from the Scalar Spherical Harmonics Y m as

Z m = Y m er, (3.25)
Y m = r ∇Y m , (3.26)
X m = ∇Y m ∧ r. (3.27)
Note that several conventions are used in the literature to define the Vector Spherical Harmonics and they are actually defined using a different convention in [START_REF] Treussart | Étude expérimentale de l'effet laser dans des micro-sphères de silice dopées avec des ions néodyme[END_REF]. Here, the radial distance r and radial vector r are included in (3.25)-(3.27) so as to guarantee that the dimension of the Vector Spherical Harmonics are the same as the ordinary Spherical Harmonics and that the Vector Spherical Harmonics do not depend on the radial spherical coordinate. We find that the three independent solutions of the vector wave equation (3.17) are expressed in terms of the Vector Spherical Harmonics as

L m (r, θ, ϕ) = ∂ ∂r z (k r) Y m (θ, ϕ) er + z (k r) ∇Y m (θ, ϕ) = kz (k r) Z m (θ, ϕ) + 1 r z (k r) Y m (θ, ϕ), (3.28)
M m (r, θ, ϕ) = z (k r) ∇Y m (θ, ϕ) ∧ r = z (k r) X m (θ, ϕ), (3.29)
N m (r, θ, ϕ) = ( + 1) kr z (k r) Y m (θ, ϕ) er + ∂ ∂r r z (k r) ∇Y m (θ, ϕ) = ( + 1) kr z (k r) Z m (θ, ϕ) + 1 kr ∂ ∂r r z (k r) Y m (θ, ϕ) (3.30)
where z denotes either the spherical Bessel function of first or second kind. The solution to the vector wave equation (3.17) can be represented as

C(r, θ, ϕ) = ∞ =0 m=- α m L m (r, θ, ϕ) + β m M m (r, θ, ϕ) + γ m N m (r, θ, ϕ) (3.31)
where α m , β m and γ m denotes complex valued constants.

TE and TM whispering gallery modes in a micro-sphere

When the desired solution to the vector wave equation (3.17) is solenoidal, which is the case in the present study for the magnetic induction B and electric field E inside and outside the sphere, see equations (2.9b) and (2.9d), it can be expanded in terms of M m and N m only [START_REF] Stratton | Electromagnetic Theory[END_REF]. We can therefore distinguish two types of particular electromagnetic fields to our problem:

(i) Transverse Electric (TE) modes where the electric field E is collinear to M m , i.e.

E(r, θ, ϕ) = A TE i/o M m (r, θ, ϕ) (3.32)
where A TE i/o denotes a complex constant number with a different value depending on the domain (index i for inside or index o for outside the sphere).

(ii) Transverse Magnetic (TM) modes where the electric field E is collinear to N m , i.e.

E(r, θ, ϕ) = A TM i/o N m (r, θ, ϕ) (3.33)
where A TM i/o denotes a complex constant number with a different value depending on the domain (index i for inside or index o for outside the sphere).

Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr From equation (2.9a), the magnetic induction is given by B = i ω curl E so that from relations (3.23) and (3.24) we deduce that: (i) for TE modes the magnetic induction B is collinear to ik ω N m , i.e.

B(r, θ, ϕ) = A TE i/o ik ω N m (r, θ, ϕ); (3.34) (ii) for TM modes the magnetic induction B is collinear to ik ω M m , i.e. B(r, θ, ϕ) = A TM i/o ik ω M m (r, θ, ϕ). (3.35)
Since M m is orthogonal to the unit radial vector er, in TE modes the electric field is parallel to the surface of the microsphere whereas in TM modes, the magnetic induction is parallel to the surface of the microsphere. We conclude that the electromagnetic field for TE modes has the following form:

E(r, θ, ϕ) = A TE i/o z (kr) X m (θ, ϕ) = A TE i/o r (kr) kr X m (θ, ϕ), (3.36) B(r, θ, ϕ) = A TE i/o ik ω ( + 1) kr z (k r) Z m (θ, ϕ) + 1 kr ∂ ∂r r z (k r) Y m (θ, ϕ) = A TE i/o ik ω ( + 1) r (k r) k 2 r 2 Z m (θ, ϕ) + r (k r) kr Y m (θ, ϕ) , (3.37)
and for TM modes it has the form:

E(r, θ, ϕ) = A TM i/o ( + 1) kr z (k r) Z m (θ, ϕ) + 1 kr ∂ ∂r r z (k r) Y m (θ, ϕ) = A TM i/o ( + 1) r (k r) k 2 r 2 Z m (θ, ϕ) + r (k r) kr Y m (θ, ϕ) , (3.38) B(r, θ, ϕ) = A TM i/o ik ω z (kr) X m (θ, ϕ) = A TM i/o ik ω r (kr) kr X m (θ, ϕ), (3.39)
where z denotes the spherical Bessel function of first kind j or second kind y and r denotes the Ricatti-Bessel function of first kind ψ or second kind χ and and m denote the mode numbers.

The Ricatti-Bessel functions of second kind like the Bessel functions of the second kind have a singularity at the origin (i.e. at the sphere center). They are therefore not suited to describe the electromagnetic field inside the microsphere. Inside the sphere, TM and TE modes are expressed in terms of the Ricatti-Bessel functions of first kind alone, i.e. r = ψ .

Outside the sphere, TM and TE modes coincide with waves that propagates outward from the sphere. With the +iωt convention used here for sinusoidal varying fields, see (2.8), such waves are expressed as the following linear combination of the Ricatti-Bessel functions of first kind and second kind corresponding to the so-called second Ricatti-Bessel function of the third kind

ζ = ψ + iχ (3.40)
Finally, for TE modes the electromagnetic field is given in spherical coordinates by

E(r, θ, ϕ) =        A TE i ψ (kr) kr X m (θ, ϕ) if r R -A TE o ζ (k0r) k0r X m (θ, ϕ) if r > R (3.41)
Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr and

B(r, θ, ϕ) =          A TE i ik ω ( + 1) ψ (k r) k 2 r 2 Z m (θ, ϕ) + ψ (k r) kr Y m (θ, ϕ) if r R -A TE o ik0 ω ( + 1) ζ (k0 r) k 2 0 r 2 Z m (θ, ϕ) + ζ (k0 r) k0r Y m (θ, ϕ) if r > R (3.42)
where X m , Y m and Z m denote the vector spherical harmonics and ψ and ζ denote respectively the Ricatti-Bessel functions of first and third kinds. For TM modes the electromagnetic field reads

E(r, θ, ϕ) =          A TM i ( + 1) ψ (k r) k 2 r 2 Z m (θ, ϕ) + ψ (k r) kr Y m (θ, ϕ) if r R -A TM o ( + 1) ζ (k0 r) k 2 0 r 2 Z m (θ, ϕ) + ζ (k0 r) k0r Y m (θ, ϕ) if r > R (3.43) and B(r, θ, ϕ) =        A TM i ik ω ψ (kr) kr X m (θ, ϕ) if r R -A TM o ik0 ω ζ (k0r) k0r X m (θ, ϕ) if r > R (3.44)

The modal equation

The electromagnetic field inside and outside the microsphere are connected through the interface conditions at the sphere boundary deduced from boundary conditions (2.7a)-(2.7d). They read

∀(θ, ϕ) ∈ [0, π] × [0, 2π] E(R, θ, ϕ) ∧ er = 0, (3.45a) B(R, θ, ϕ) • er = 0, (3.45b) ε E(R, θ, ϕ) • er = 0, (3.45c) B(R, θ, ϕ) ∧ er = 0. (3.45d)
In the next two sections we will study the conditions at which equations (3.45a)-(3.45d) are satisfied by TE modes and TM modes respectively.

3.4.1.

Modal equation for TE modes According to (3.32) and (3.34), the electromagnetic field (E, B) for TE modes is given by

E(R, θ, ϕ) = A TE i/o r (kR) kR X m (θ, ϕ), (3.46)
B(R, θ, ϕ) = A TE i/o ik ω ( + 1) r (k R) k 2 R 2 Z m (θ, ϕ) + r (k R) kR Y m (θ, ϕ) . (3.47)
Since X m is orthogonal to er, condition (3.45c) is always satisfied. Condition (3.45a) implies that we must have

A TE i/o r (kR) kR = 0 (3.48)
that is to say

A TE i ψ (kR) kR = -A TE o ζ (k0R) k0R . (3.49) Condition (3.45b) implies that A TE i/o ( + 1) r (k R) kR 2 Z m (θ, ϕ) • er + A TE i/o r (k R) R Y m (θ, ϕ) • er = 0
Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr and thanks to (3.48), the condition reduces to

A TE i/o r (kR) = 0 (3.50)
that is to say to

A TE i ψ (kR) = -A TE o ζ (k0R). (3.51) Condition (3.45d
) is identically satisfied when relations (3.48) and (3.50) are taken into account. Thus, boundary conditions for TE modes give rise to the following system of equations to be fulfilled

   A TE i ψ (kR) kR + A TE o ζ (k0R) k0R = 0 A TE i ψ (kR) + A TE o ζ (k0R) = 0 . (3.52)
Whenever the determinant of the linear system is nonzero, its unique solution is the zero solution. There exists nontrivial solutions if and only if the determinant is zero, that is to say for k such that k k0

ψ (kR) ψ (kR) = ζ (k0R) ζ (k0R) . (3.53)
In such a case, there exists an infinite set of solutions to the linear system (3.52) in the form

(A TE i , A TE o ) where A TE o = -A TE i k0 k ψ (kR) ζ (k0R) . (3.54)
The value of A TE i for a given experiment can be computed from the condition that the total power in the device is imposed or known. Equation (3.53) is referred as the modal equation for TE whispering gallery modes.

Modal equation for TM modes

E(R, θ, ϕ) = A TM i/o ( + 1) r (k R) k 2 R 2 Z m (θ, ϕ) + r (k R) kR Y m (θ, ϕ) , (3.55)
B(R, θ, ϕ) = A TM i/o i ω r (kR) R X m (θ, ϕ). (3.56)
Since X m is orthogonal to er, condition (3.45b) is always satisfied. Condition (3.45d) implies that we must have

A TM i/o r (kR) = 0 (3.57)
that is to say

A TM i ψ (kR) = -A TM o ζ (k0R). (3.58) Condition (3.45a) implies A TM i/o ( + 1) r (k R) k 2 R 2 Z m (θ, ϕ) ∧ er + A TM i/o r (k R) kR Y m (θ, ϕ) ∧ er = 0.
Since Z m is collinear to er, the first term in the brackets is always zero and the condition reduces to

A TM i/o r (kR) k = 0 (3.59) that is to say k0 A TM i ψ (kR) = -k A TM o ζ (k0R). ( 3
A TM i ψ (kR) + A TM o ζ (k0R) = 0 A TM i k0 ψ (kR) + A TM o k ζ (k0R) = 0 . (3.61)
Whenever the determinant of the linear system is nonzero, its unique solution is the zero solution. There exists nontrivial solutions if and only if the determinant is zero, that is to say when k is such that k0 k

ψ (kR) ψ (kR) = ζ (k0R) ζ (k0R) . (3.62)
In such a case, there exists an infinite set of solutions to the linear system (3.61) in the form

(A TM i , A TM o )
where Note that this underscores a mode degeneracy since there exists 2 + 1 modes with the same values of and n, and therefore with the same resonance frequency, but with a different value for m and therefore with a different expression for the electromagnetic field. It can be shown that the mode number is equal to the number of wavelengths taken to travel around the sphere and that the radial mode number n is equal to the number of intensity maxima of the mode in the radial direction er The index m is called the azimuthal mode number. It can take 2 + 1 values from -to and it is related to the sinusoidal variation of the mode with the azimuthal angle ϕ. Moreover, -|m| + 1 is the number of intensity maxima in the polar direction e θ . Thus modes with index n = 1 correspond to the best confined modes in the radial direction and modes for which m = are the best confined in the polar direction. The mode satisfying these conditions and corresponding to the highest value of for which the modal equation (3.65) has a solution is termed the fundamental mode.

A TM o = -A TM i ψ (kR) ζ (k0R) . ( 3

Simplification of the modal equations

We have established in the previous section that the modal equations for TE and TM modes were in the form In practice, for most applications, the micro-sphere radius R is large relative to the wavelength λ of the optical wave (and therefore the mode index is large, see remark 1). As Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr a consequence, the radiative leakage of the energy is small and a good approximation consists in assuming the radiative part outside the sphere as negligible. This approximation leads to consider the following form for the modal equation [START_REF] Righini | Whispering gallery mode microresonators: Fundamentals and applications[END_REF][START_REF] Treussart | Étude expérimentale de l'effet laser dans des micro-sphères de silice dopées avec des ions néodyme[END_REF]

P ψ (kR) ψ (kR) = ζ (k0R) ζ (
P ψ (kR) ψ (kR) = χ (k0R) χ (k0R) (3.65)
where χ denotes Riccati-Bessel function of the second kind with order . This approximation relies on asymptotic expansions of Bessel's functions for large order, see e.g. [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF] formula 10.19.2 that implies that ζ (z) ≈ iχ (z) for large . One should note that this assumption implies that Im(k) = 0 and k ∈ R.

Numerical resolution of the modal equations

For numerical purposes it is convenient to express the modal equation (3.65) in terms of Bessel function of first and third kind. For all ν ∈ R and x ∈ C we have [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF][START_REF] Gradstein | Table of Integrals, Series, and Products. Table of Integrals, Series, and Products Series[END_REF]

ψν (x) = πx 2 J ν+ 1 2 (x) (3.66) χν (x) = πx 2 Y ν+ 1 2 (x) (3.67)
It follows that

ψ ν (x) = π 8x J ν+ 1 2 (x) + πx 2 J ν+ 1 2 (x), (3.68)
χ ν (x) = π 8x Y ν+ 1 2 (x) + πx 2 Y ν+ 1 2 (x). ( 3
J ν (x) = Jν-1(x) - ν x Jν (x) and Y ν (x) = Yν-1(x) - ν x Yν (x)
so that we can express the modal equation (3.65) as

Y -1 2 (k0R) Y + 1 2 (k0R) -P J -1 2 (kR) J + 1 2 (kR) = 1 k0R - P kR . (3.71)
For TE modes we have P = N so that the second hand side of equation (3.71) is zero, whereas for TM modes we have P = 1/N and therefore the second hand side of equation (3.71) is

k 0 R (1 -1 N 2 )
. For a given mode number , the wavelengths λ for which resonance occurs are obtained by solving equation (3.71). That is to say they are obtained by looking for the zeros of the modal function

F : λ ∈ ]0, +∞[ -→ Y -1 2 (k0R) Y + 1 2 (k0R) -P J -1 2 (kR) J + 1 2 (kR) - 1 k0R - P kR . (3.72)
It is known, see [START_REF] Treussart | Étude expérimentale de l'effet laser dans des micro-sphères de silice dopées avec des ions néodyme[END_REF], that for a given λ, the mode numbers are such that

2π(R + δP ) λ < + 1 2 < N 2π(R + δP ) λ (3.73)
Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

where δP ≈ λ 2πN P √ N 2 -1 .
The bounds for given by (3.73) imply that for a fixed value of the wavelengths λ for which resonance occurs are such that 2πR

+ 1 2 -α < λ < 2πRN + 1 2 -αN (3.74)
where we have set α = P/(N √ N 2 -1). In this study, typical values for the physical parameters are λ = 810 nm for the wavelength, N = 1.45 for the optical index of the micro-sphere and R = 50 µm for the radius of the micro-sphere. The mode number varies between 194 and 281, see [START_REF] Treussart | Étude expérimentale de l'effet laser dans des micro-sphères de silice dopées avec des ions néodyme[END_REF]. We have k0 = 7.75 10 6 and k = 1.12 10 (kR) involve quantities with very comparable values which is a condition for high accuracy results.

We can speed up the computations for solving equation (3.71) by using the following recurrence relation for Bessel functions of first kind [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF][START_REF] Gradstein | Table of Integrals, Series, and Products. Table of Integrals, Series, and Products Series[END_REF]

Jν-1(x) + Jν+1(x) = 2ν x Jν (x). (3.75)
We have the same recurrence relation for Bessel functions of second kind. If we define the sequence of functions (uν )ν such that uν = Jν+1/Jν , then for all x ∈ R where x do not coincide with a zero of the Bessel functions of first kind, the values of uν (x) can be computed from the value of u0(x) through the recurrence relation

uν (x) = 2ν x - 1 uν-1(x) . ( 3

Numerical results

A first numerical comparison with existing results

In F. Treussart Phd thesis [START_REF] Treussart | Étude expérimentale de l'effet laser dans des micro-sphères de silice dopées avec des ions néodyme[END_REF] is given a micro-sphere of raduis 10 µm and optical index 1.45 is considered on p. 34 and the TE wave functions for the index = 100 are drawn. We give in Table 1 the values of the resonance frequencies for such a micro-sphere. The higher frequency (n = 1) is found to be 844.4558 nm. The corresponding value of the size parameter x = 2πR/λ is 74.405141 nm. The value proposed in [START_REF] Treussart | Étude expérimentale de l'effet laser dans des micro-sphères de silice dopées avec des ions néodyme[END_REF] is 74.4064. The 5th higher frequency (n = 5) is found to be 697.3817 nm. The corresponding value of the size parameter x = 2πR/λ is 90.096790 nm. The value proposed in [START_REF] Treussart | Étude expérimentale de l'effet laser dans des micro-sphères de silice dopées avec des ions néodyme[END_REF] is 90.0955. The values in F. Treussart Phd thesis [START_REF] Treussart | Étude expérimentale de l'effet laser dans des micro-sphères de silice dopées avec des ions néodyme[END_REF] are therefore in a good agreement with the one obtained by solving (3.65) by the method outlined in section 3.6.

We have depicted in Figure 3 the variation of the electric field as a function of the normalized radial distance r/R. These figures are in a good agreement with the figures given in F. Treussart Phd thesis [START_REF] Treussart | Étude expérimentale de l'effet laser dans des micro-sphères de silice dopées avec des ions néodyme[END_REF]

An additional numerical experiment

We consider a micro-sphere of radius R = 25 µm and optical index N = 1.453. We proceed in two steps. First, assume that we are interested with resonance conditions for a wavelength λ of the optical wave around 810 nm. We solve the modal equation We can now start the second computational step. We fix the mode number to be = 271 and we look for the wavelength around 810 nm for which the resonance condition is realized. We solve the modal equation (3.71) for = 271 in order to determine the values of λ for which the resonance occurs. Equation (3.71) is numerically solved with Matlab using the fsolve command. The bounds for the mode numbers given by (3.74) indicate that the area of interest for λ in units of nm is the interval [580.5908, 844.9405].

F λ : ∈ ]0, +∞[ -→ Y -1 2 (k0R) Y + 1 2 (k0R) -N J -1 2 (kR) J + 1 2 (kR) . ( 3
In Fig. 5 We have depicted in Fig. 6 the radial function dependency of the TE modes for several resonance wavelengths as given by Table 3. On the top line is depicted the radial function dependency of the TE modes for λ = 808.6189 nm. We have also depicted the radial function dependency of the TE modes for the 2nd, 3rd and 10th resonance wavelength to compare the behavior of the field inside the sphere when the wavelength decreases.

Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr We have depicted in Fig. 7 the euclidean norm of the electric field inside the sphere in the plane of azimuthal angle ϕ = 0 for TE modes with degree and order ( , m) = (271, 271),( , m) = (271, 270) and ( , m) = (271, 10). We can see that the intensity of the mode decreases very quickly with m.

Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

Computation of the volume of a whispering gallery mode

The volume of a whispering gallery mode in a micro-sphere is defined, see [START_REF] Kippenberg | Nonlinear Optics in Ultra-high-Q Whispering-Gallery Optical Microcavities[END_REF][START_REF] Treussart | Étude expérimentale de l'effet laser dans des micro-sphères de silice dopées avec des ions néodyme[END_REF], as the integral over the whole space of the energy density normalized by its maximum value inside Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr the micro-sphere, i.e.

V = 1 wmax R 3 w(x1, x2, x3) dx1dx2dx3 (4.1)
where w denotes the energy density given as a function of the position vector x = (x1, x2, x3) by

w(x) = 1 2 ε(x) 2 E(x) E * (x) + 1 2µ0 B(x) B * (x) (4.2)
and wmax denotes the maximum value of the energy density inside the micro-sphere. In relation (4.2), the quantity E * (resp. B * ) stands for the adjoint (conjugate transpose) of E (resp. B) so that

E(x) E * (x) = E(x) • E(x) = E(x) 2 2 = |E1(x)| 2 + |E2(x)| 2 + |E3(x)| 2 .
When we assume that the energy losses by diffraction and diffusion can be neglected, the energy conservation law implies that the contribution of the electric field and the magnetic field to the energy density are equal, see [START_REF] Jackson | Classical Electrodynamics Third Edition[END_REF] and Appendix B. Namely, we have

R 3 ε(x) 2 E(x) 2 dx = R 3 1 2µ0 B(x) 2 dx. (4.3)
As a consequence, the volume of a WGM is most often computed in the literature [START_REF] Oraevsky | Review: Whispering-gallery waves[END_REF][START_REF] Oxborrow | Traceable 2-D Finite-Element simulation of the whispering-gallery modes of axisymmetric electromagnetic resonators[END_REF][START_REF] Kippenberg | Nonlinear Optics in Ultra-high-Q Whispering-Gallery Optical Microcavities[END_REF][START_REF] Treussart | Étude expérimentale de l'effet laser dans des micro-sphères de silice dopées avec des ions néodyme[END_REF][START_REF] Ward | WGM microresonators: sensing, lasing and fundamental optics with microspheres[END_REF] from the following formula:

V = 1 ε0N 2 E 2 max R 3 ε(x) E(x) 2 dx (4.4)
where Emax denotes the maximum value of the Euclidean norm of the electric field. Formula (4.4) is well suited for the computation of the mode volume of a TE mode in a spherical optical micro-resonator because of the simple expression of the electric field for TE mode, see (3.36). It should be noted however that the normalization constant is not the same in (4.1) and in (4.4).

For TM modes, computation of the mode volume from (4.4) is a little more tricky because it is the magnetic induction that has the simplest expression, see (3.39). Thanks to (4.3), it is however possible to express the integral in (4.4) in terms of B. Unfortunately, the normalization constant in (4.4) can't be expressed in terms of B and changing it for max x∈R 3 B(x) 2 would not be consistent when comparing volumes of TE and TM modes. Thus, we will compute the volume of TM mode from the following formula

V = c 2 N 2 E 2 max R 3 B(x) 2 dx (4.5)
where c is the speed of light in free space.

Volume of a TE mode

For TE modes we found that the electromagnetic field is given in spherical coordinates by

E(r, θ, ϕ) =        A TE i ψ (kr) kr X m (θ, ϕ) if r < R -A TE o χ (k0r) k0r X m (θ, ϕ) if r > R (4.6)
and

B(r, θ, ϕ) =          A TE i ik ω ( + 1) ψ (k r) k 2 r 2 Z m (θ, ϕ) + ψ (k r) kr Y m (θ, ϕ) if r < R -A TE o ik ω ( + 1) χ (k0 r) k 2 0 r 2 Z m (θ, ϕ) + χ (k0 r) k0r Y m (θ, ϕ) if r > R (4.7)
Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

where X m , Y m and Z m denote the vector spherical harmonics as defined by (3.25)-(3.27) and ψ and χ denote respectively the Ricatti-Bessel functions of first and second types.

One should note that due to the assumption introduced in Section 3.5 p. 17, k is real and k > 0.

4.1.1.

Computation of the volume integral Since the electric field decays very quickly outside the micro-sphere, the integral (B.6) in spherical coordinates can be approximated by

V = 1 ε0N 2 L 0 π 0 2π 0 ε(r) E(r, ϕ, θ) 2 2 E 2 max r 2 sin(θ) dϕdθdr (4.8)
where L is a positive number, large enough so that the electromagnetic field can be neglected at a radial distance greater than L. Using Chasles's theorem, we split the integration domain in two parts such that the first integral is over the domain inside the micro-sphere and the second integral is over the exterior domain. Accordingly we set V = Vi +Vo. The first integral over the interior domain reads

Vi = 1 k 2 |A TE i | 2 E 2 max R 0 π 0 2π 0 ψ 2 (kr) X m (ϕ, θ) 2 2 sin(θ) dϕdθdr = 1 k 2 |A TE i | 2 E 2 max R 0 ψ 2 (kr) dr π 0 2π 0 X m (ϕ, θ) 2 2 sin(θ) dϕdθ . (4.9)
The second integral term in (4.9) is known to be

π 0 2π 0 X m (ϕ, θ) • X * m (ϕ, θ) sin(θ) dϕdθ = ( + 1).
Therefore, we have

Vi = ( + 1) k 2 |A TE i | 2 E 2 max R 0 ψ 2 (kr) dr = ( + 1) k |A TE i | 2 E 2 max π 2 R 0 r J 2 + 1 2 (kr) dr. ( 4.10)
The integral in (4.10) can be evaluated exactly thanks to the following formula, see [START_REF] Gradstein | Table of Integrals, Series, and Products. Table of Integrals, Series, and Products Series[END_REF] formula 5.54 p. 629,

x Jp(αx)

2 dx = x 2 2 Jp(αx) 2 -Jp-1(αx)Jp+1(αx) .
We conclude that Vi = ( + 1)πR 2 4k

|A TE i | 2 E 2 max J + 1 2 (kR) 2 -J -1 2 (kR)J + 3 2 (kR) . (4.11)
The integral over the exterior domain reads

Vo = 1 N 2 |A TE o | 2 E 2 max L R π 0 2π 0 χ 2 (k0r) k 2 0 X m (ϕ, θ) 2 2 sin(θ) dϕdθdr = 1 N 2 k 2 0 |A TE o | 2 E 2 max L R χ 2 (k0r) dr π 0 2π 0 X m (ϕ, θ) 2 2 sin(θ) dϕdθ = ( + 1) N 2 k 2 0 |A TE o | 2 E 2 max L R χ 2 (k0r) dr = ( + 1) N 2 k0 |A TE o | 2 E 2 max π 2 L R r Y 2 + 1 2 (k0r) dr (4.12)
Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

The integral in (4.12) can be evaluated exactly thanks to the following formula, see [START_REF] Gradstein | Table of Integrals, Series, and Products. Table of Integrals, Series, and Products Series[END_REF] formula 5.54 p. 629,

x Yp(αx)

2 dx = x 2 2 Yp(αx) 2 -Yp-1(αx)Yp+1(αx) .
We conclude that

Vo = ( + 1)π 4N 2 k0 |A TE o | 2 E 2 max L 2 Y + 1 2 (k0L) 2 -L 2 Y -1 2 (k0L)Y + 3 2 (k0L) -R 2 Y + 1 2 (k0R) 2 + R 2 Y -1 2 (k0R)Y + 3 2 (k0R) . (4.13)
For TE mode, according to relation (3.54) we have

A TE o = -A TE i k0 k ψ (kR) χ (k0R) = A TE i k0 k J + 1 2 (kR) Y + 1 2 (k0R)
.

Since k = N k0, we finally conclude that

Vo = ( + 1)π 4N 3 k0 |A TE i | 2 E 2 max J + 1 2 (kR) 2 Y + 1 2 (k0R) 2 L 2 Y + 1 2 (k0L) 2 -L 2 Y -1 2 (k0L)Y + 3 2 (k0L) -R 2 Y + 1 2 (k0R) 2 + R 2 Y -1 2 (k0R)Y + 3 2 (k0R) (4.14) = ( + 1)π 4N 3 k0 |A TE i | 2 E 2 max J + 1 2 (kR) 2 L 2 Y + 1 2 (k0L) 2 Y + 1 2 (k0R) 2 -L 2 Y -1 2 (k0L)Y + 3 2 (k0L)
Y + 1 2 (k0R) 2 -R 2 + R 2 Y -1 2 (k0R)Y + 3 2 (k0R)
Y + 1 2 (k0R) 2 (4.15)
4.1.2. Computation of the maximum value of the electric field inside the micro-sphere Inside the micro-sphere at a point with spherical coordinates (r, θ, ϕ), the electric field for a TE mode reads

E(r, θ, ϕ) = A TE i ψ (kr) kr X m (θ, ϕ).
The maximum value of the euclidean norm of the electric field over the micro-sphere is given by

Emax = A TE i sup (r,θ,ϕ)∈[0,R]×[0,π]×[0,2π] ψ (kr) kr X m (θ, ϕ) 2 = A TE i sup r∈[0,R] ψ (kr) kr × sup (θ,ϕ)∈[0,π]×[0,2π] X m (θ, ϕ) 2 (4.16)
We first compute the maximum over [0, R] of |j (kr)| = |ψ (kr)/kr|. It follows from the known behavior of the spherical Bessel function of the first kind, see e.g. [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF][START_REF] Gradstein | Table of Integrals, Series, and Products. Table of Integrals, Series, and Products Series[END_REF], that the global maximum of |j | coincides with the first local maximum of j . Thus, we are looking for the first zero of the derivative of the spherical Bessel function of the first kind and order denoted j . For all x ∈ R, we have [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF][START_REF] Gradstein | Table of Integrals, Series, and Products. Table of Integrals, Series, and Products Series[END_REF]

j (x) = j -1 (x) -( + 1) j +1 (x) 2 + 1 and j (x) = 0 ⇐⇒ j -1 (x) -( + 1)j +1 (x) = 0 ⇐⇒ ( + 1) J + 3 2 (x) -J -1 2 (x) = 0. (4.17)
Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

Thus the problem turns out to solving the nonlinear equation (4.17). More precisely, we have to compute the first positive root of equation (4.17). Numerical methods for solving a nonlinear equation require either the knowledge of an interval where the root is isolated, or a first guess for this solution. An approximation based on an asymptotic expansion for the first zero (denoted a ,1 ) of j for large is given in [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF] (see formula 10.1.59 on p. 441):

a ,1 = ( + 1 2 ) + 0.8086165 ( + 1 2 ) 1/3 -0.236680 ( + 1 2 ) -1/3 -0.20736 ( + 1 2 ) -1 + 0.0233 ( + 1 2 ) -5/3 . (4.18)
Actually, a ,1 can be considered as a good initial guess for the first zero of j even for small . For instance for = 5 we have a ,1 = 6.7606 whereas an accurate value for the first zero of j computed with the symbolic computation software Maple is found to be 6.7564. Thus, from a computational point of view we use the matlab command fsolve to solve (4.17) with an initial guess for the solution given by (4.18).

We then have to compute the maximum value of the euclidean norm of the vectorial spherical harmonic X m inside the micro-sphere. We have X m = ∇Y m ∧ r and

Y m (θ, ϕ) = C ,m P m (cos(θ)) e imϕ where C ,m = (2 + 1)( -m)! 4π( + m)! .
It follows that

∇Y m (θ, ϕ) ∧ r = 1 sin(θ) ∂Y m ∂ϕ (θ, ϕ) e θ - ∂Y m ∂θ (θ, ϕ) eϕ = C ,m e imϕ im sin(θ) P m (cos(θ)) e θ + sin(θ) (P m ) (cos(θ)) eϕ and X m (θ, ϕ) 2 2 = C 2 ,m m 2 sin 2 (θ) P m (cos(θ)) 2 + sin 2 (θ) (P m ) (cos(θ)) 2 .
We are looking for the maximum of X m (θ, ϕ) 2 2 , which is actually independent of ϕ, for all θ ∈ [0, π]. The derivative of the associated Legendre function P m is given by, see [START_REF] Gradstein | Table of Integrals, Series, and Products. Table of Integrals, Series, and Products Series[END_REF] (formula 8.733 p. 965),

∀x ∈ ] -1, 1[ (P m ) (x) = 1 1 -x 2 P m (x) -( -m + 1) P m +1 (x) .
Therefore we have

X m (θ, ϕ) 2 2 = C 2 ,m sin 2 (θ) (m 2 + 1) P m (cos(θ)) 2 + ( -m + 1) 2 P m +1 (cos(θ)) 2 -2( -m + 1) P m (cos(θ)) P m +1 (cos(θ)) .
Unfortunately, the previous expression is not well suited for numerical computation purposes. Indeed for larger than 150 the values of the associated Legendre function P m is not anymore representable in the floating point arithmetic in Matlab. As a consequence, we introduce the Schmidt semi-normalized associated Legendre functions S m instead of the associated Legendre function P m . They are related to each others by the relation

∀x ∈] -1, 1[ S m (x) =      P (x) if m = 0 (-1) m 2( -m)! ( + m)! P m (x) if m > 0 . (4.19)
Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

We set

D ,m =      1 if m = 0 (-1) m ( + m)! 2( -m)! if m > 0
so that for all x ∈] -1, 1[ we have P m (x) = D ,m S m (x). Moreover one can readily check that

D +1,m = D ,m + 1 + m + 1 -m .
In terms of the Schmidt semi-normalized associated Legendre functions, we have

X m (θ, ϕ) 2 2 = C 2 ,m D 2 ,m sin 2 (θ) (m 2 + 1) S m (cos(θ)) 2 + (( + 1) 2 -m 2 ) S m +1 (cos(θ)) 2 -2 ( + 1) 2 -m 2 S m (cos(θ)) S m +1 (cos(θ)) (4.20)
where

C 2 ,m D 2 ,m =      2 + 1 4π if m = 0 2 + 1 8π if m > 0 .
Since there is no obvious way for determining the maximum of X m (θ, 0) 2 2 as given by (4.20), we proceed by using a brute force method that consists in comparing the values of X m (θ, 0) 2 2 over a sufficiently accurate subdivision of the interval [0, π].

Volume of a TM mode

For TM modes the electromagnetic field reads

E(r, θ, ϕ) =          A TM i ( + 1) ψ (k r) k 2 r 2 Z m (θ, ϕ) + ψ (k r) kr Y m (θ, ϕ) if r < R -A TM o ( + 1) χ (k0 r) k 2 0 r 2 Z m (θ, ϕ) + χ (k0 r) k0r Y m (θ, ϕ) if r > R (4.21) and B(r, θ, ϕ) =        A TM i ik ω ψ (kr) kr X m (θ, ϕ) if r < R -A TM o ik0 ω χ (k0r) k0r X m (θ, ϕ) if r > R (4.22)
where X m , Y m and Z m denote the vector spherical harmonics as defined by (3.25)-(3.27) and ψ and χ denote respectively the Ricatti-Bessel functions of first and second types.

Computation of the volume integral

The integral (4.5) in spherical coordinates reads

V = L 0 π 0 2π 0 B(x) 2 2 E 2 max r 2 sin(θ) dϕdθdr (4.23)
where L is a positive number, large enough so that the electromagnetic field can be neglected at a radial distance greater than L. Using Chasles's theorem, we split the integration domain in two parts V = Vi + Vo such that the first integral holds over the domain inside the microsphere and the second integral holds over the exterior domain. The first integral reads Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

Vi = 1 k 2 |A TM i | 2 E 2 max R 0 π 0 2π 0 ψ 2 (kr) X m (ϕ, θ) 2 2 sin(θ) dϕdθdr = 1 k 2 |A TM i | 2 E 2 max R 0 ψ 2 (kr) dr π 0 2π 0 X m (ϕ, θ)
The second integral term in (4.24) is known to be [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF][START_REF] Gradstein | Table of Integrals, Series, and Products. Table of Integrals, Series, and Products Series[END_REF] π 0 2π 0 X m (ϕ, θ) X * m (ϕ, θ) sin(θ) dϕdθ = ( + 1).

Therefore we have

Vi = ( + 1) k 2 |A TM i | 2 E 2 max R 0 ψ 2 (kr) dr = k ( + 1) k 2 |A TM i | 2 E 2 max π 2 R 0 r J 2 + 1 2 (kr) dr. (4.25)
The integral in (4.25) can be evaluated exactly thanks to the following formula, see [START_REF] Gradstein | Table of Integrals, Series, and Products. Table of Integrals, Series, and Products Series[END_REF] formula 5.54 p. 629,

x Jp(αx)

2 dx = x 2 2 Jp(αx) 2 -Jp-1(αx)Jp+1(αx) .
We get Vi = ( + 1)πR 2 4k

|A TM i | 2 E 2 max J + 1 2 (kR) 2 -J -1 2 (kR)J + 3 2 (kR) .
(4.26)

The integral over the domain exterior to the sphere reads

Vo = c 2 N 2 R 3 B(x) 2 2 E 2 max dx = 1 k 2 |A TM o | 2 E 2 max L R π 0 2π 0 χ 2 (k0r) X m (ϕ, θ) 2 2 sin(θ) dϕdθdr = 1 k 2 |A TM o | 2 E 2 max L R χ 2 (k0r) dr π 0 2π 0 X m (ϕ, θ) 2 2 sin(θ) dϕdθ = ( + 1) k 2 |A TM o | 2 E 2 max L R χ 2 (k0r) dr = ( + 1) k 2 |A TM o | 2 E 2 max k0π 2 L R r Y 2 + 1 2 (k0r) dr. (4.27)
The integral in (4.27) can be evaluated exactly thanks to the following formula, see [START_REF] Gradstein | Table of Integrals, Series, and Products. Table of Integrals, Series, and Products Series[END_REF] formula 5.54 p. 629,

x Yp(αx)

2 dx = x 2 2 Yp(αx) 2 -Yp-1(αx)Yp+1(αx) .
We conclude that

Vo = ( + 1)πk0 4k 2 |A TM o | 2 E 2 max L 2 Y + 1 2 (k0L) 2 -L 2 Y -1 2 (k0L)Y + 3 2 (k0L) -R 2 Y + 1 2 (k0R) 2 + R 2 Y -1 2 (k0R)Y + 3 2 (k0R) . (4.28)
For TM mode, according to relation (3.54) we have

A TM o = -A TM i ψ (kR) χ (k0R) = A TM i k k0 J + 1 2 (kR) Y + 1 2 (k0R) .
Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr Therefore Vo = ( + 1)π 4k max the maximum value of the euclidean norm of the electric field. Since the electromagnetic field for a WGM is at its peak inside the cavity, we deduce from (3.38) that

|A TM i | 2 E 2 max J + 1 2 (kR) 2 Y + 1 2 (k0R) 2 L 2 Y + 1 2 (k0L) 2 -L 2 Y -1 2 (k0L)Y + 3 2 (k0L) -R 2 Y + 1 2 (k0R) 2 + R 2 Y -1 2 (k0R)Y + 3 2 (k0R) (4.29) = ( + 1)π 4k |A TM i | 2 E 2 max J + 1 2 (kR) 2 L 2 Y + 1 2 (k0L) 2 Y + 1 2 (k0R) 2 -L 2 Y -1 2 (k0L)Y + 3 2 (k0L) Y + 1 2 (k0R) 2 -R 2 + R 2 Y -1 2 (k0R)Y + 3 2 (k0R) Y + 1 2 (k0R) 2 . ( 4
E 2 max (r, θ) = C 2 ,m k 2 r 2 2 ( + 1) 2 |ψ (kr)| 2 k 2 r 2 P m (cos θ) 2 + |ψ (kr)| 2 sin(θ)(P m ) (cos θ) 2 + m sin(θ) P m (cos θ) 2 . (4.31)
Analytical computation of the maximum value of E 2 max over [0, R] × [0, π] is out of reach. Therefore, in the Matlab toolbox we proceed by using numerical optimization methods and to this purpose, we express E 2 max in terms of Schmidt semi-normalized associated Legendre functions S m and Bessel's functions:

E 2 max (r, θ) = C 2 ,m D 2 ,m k 2 r 2 π 2 2 ( + 1) 2 |J + 1 2 (kr)| 2 kr S m (cos θ) 2 + (kr) 1 2 J -1 2 (kr) -(kr) -1 2 J + 1 2 (kr) 2 (m 2 + 1) S m (cos θ) 2 (4.32) -2 ( + 1) 2 -m 2 S m (cos θ) S m +1 (cos θ) + (( + 1) 2 -m 2 ) S m +1 (cos θ) 2 .
In order to compute the global maximum of E 2 max over [0, R] × [0, π] we use Matlab fminunc routine designed to find minimum of unconstrained multivariable functions. Namely, a Trust Region Method algorithm is chosen which requires to provide the gradient of E 2 max . The initial guess for the optimization algorithm is obtained by comparing the values of E 2 max over a coarse grid of the domain [0, R] × [0, π].

Numerical results

We have written a Matlab program to compute the volume of any whispering gallery (WG) TE or TM mode on the basis of the formulas obtained in the last section. The program is available from the author. The Matlab script entitled volmod is documented in the next section. To illustrate the numerical computation of the volume of a WG mode, we have considered again a sphere of radius R = 25 µm and optical index N = 1.453. According to the results presented in section 3.7, we have considered the TE mode with indexes ( , m) = (271, 271) for a wavelength λ = 808.62 nm. The exterior integration domain was bounded at a distance L = 1.2 R. The mode volume was found to be 319.216 µm 3 . This value is very closed to the one obtained by F. Treussart in his Phd Thesis by means of asymptotic expansions formulas. Using formula (1.58) of Treussart Phd [START_REF] Treussart | Étude expérimentale de l'effet laser dans des micro-sphères de silice dopées avec des ions néodyme[END_REF], we find an Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr approximated mode volume of 315.4740 µm 3 . The detail of the calculations leading to F. Treussart approximate formula is given in Appendix.

The main advantage of the present work is that no assumption on the values of the mode indexes and m is done. For instance, we have obtained that the WG TE-mode volume for ( , m) = (271, 270) is 433.353 µm 3 and for ( , m) = (271, 250) it is 602.255 µm 3 . In Fig. 8, we have depicted the variation of the volume (in units of µm 3 ) of the WG mode = 271 as a function of index m ∈ [0, ]. We have to mention that the computation of the volume of a WG mode is very fast. For instance, computations leading to Fig. 8 which necessitated the computation of the volume of 272 modes took 75 s on a AMD-A8 Personal desktop Computer with 8 Go RAM. This means that the computation of the volume for one mode took less that 0.3 s. Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

WGMode : a Matlab Toolbox dedicated to the study of whispering gallery modes in optical micro-spheres

The WGMode Matlab Toolbox is constituted of various Matlab scripts to study whispering gallery modes in optical micro-spheres. It is governed by the CeCILL-C license under French law and abiding by the rules of distribution of free software. It can be used, modified and/or redistributed under the terms of the CeCILL-C license as circulated by CEA, CNRS and INRIA at the following URL http://www.cecill.info. The scripts of the WGMode Matlab Toolbox allow to explore resonance conditions for TE or TM modes, to visualize whispering gallery TE or TM modes in a micro-sphere and to compute the volume of any given mode. Namely, the WGMode Matlab Toolbox contains the following Matlab scripts :

• ELLRES computes, for a given wavelength, the values of the mode index for which a resonance occurs

• WVLRES computes, for a given couple of mode index ( , m), the values of the wavelength for which a resonance occurs

• VOLMOD computes the volume of a whispering gallery TE or TM modes in a microsphere

• PLTMOD plot the whispering gallery modes

The toolbox also provides a collection of special functions:

• In order to illustrate the use of the WGMode Matlab Toolbox, we consider the case of a micro-sphere of radius 100 µm and optical index 1.453. We are interested in resonance conditions around a wavelength of 1480 nm.

First, we use the script entitled ellres to determine the mode index for which resonance occurs at 1480 nm. The mode indexes are not integer numbers. We round up the propound values to the nearest integer. For instance, the largest index will be 602. Then we determine the exact wavelength around 1480 nm for which resonance occurs for a mode index of 602. This can be achieved with the script entitled wvlres. The wvlres Matlab script provides a draw of the modal function which enables to localize the position of its zeros corresponding to the resonance wavelengths and a list of these values computed by mean of the fsolve Matlab command. The script wvlres also proposes to draw the radial behavior of any of the computed modes at resonance. An expert mode can be activated by setting the variable expert mode to 1. This enables the user to access extra plotting functionality such as the possibility to represent the radial behavior of a mode at a frequency other than a resonance frequency for comparison purposes.

>> wvlres

The script entitled volmod achieves the computation of the volume of a given mode. For instance, when we want to know the volume of the TE mode defined by the indexes = 602 and m = 600 we can proceed as follows: The plotting area can be delimited to the area of interest as illustrated below.

>> pltmod

Conclusion

In this report we have proceeded to a mathematical study of whispering gallery (WG) modes in an optical micro-sphere starting from the general set of Maxwell equations. We have obtained the general expression of TE and TM modes using the method of Hansen for solving the vectorial wave equation in spherical coordinates. From the knowledge of the mathematical expression of TE and TM WG modes we have obtained a general analytical expression for the volume of a WG mode in terms of the spherical Bessel functions and spherical surface harmonics. We have written a Matlab Toolbox that implements the WG mode volume formulas and allows to compute with high accuracy the volume of any WG mode whatever are the indexes of the mode. This Toolbox is of interest since up to now only estimations of the volume of modes for mode index m ≈ based on asymptotic formulas were available.

The second integral term in (A.6) is known to be π 0 2π 0 X m (θ, ϕ) • X * m (θ, ϕ) sin(θ) dϕdθ = ( + 1).

Therefore by using Chasles's theorem we obtain where a1 denotes the first zero of ψ . An approximate value of a1 for large is given by (see [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF] formula 9.5.4 p. 371):

V = ( + 1) k 2 |A T E i | 2 E 2 max R 0 ψ 2 (kr) dr + |A T E o | 2 |A T E i | 2
a1 ≈ ( + 1/2) + 1.8557571 ( + 1/2) 1/3 + 1.033150 ( + 1/2) -1/3 -0.00397 ( + 1/2) -1 -0.0908 ( + 1/2) -5/3 + 0.043 ( + 1/2) -7/3 . (A.8)

To illustrate the quality of the approximation, we have drawn the graph of the function f for a micro-sphere of radius R = 67.5 µm and optical index N = 1.46 for wavelength around 1550 nm. In Fig. A1 we have depicted the graphs of function f over [0, kL] for L = 1.2 R and ψ 2 over [0, a1] where a1 denotes the first zero of ψ , for = 387 and λ = 1548.93 nm (it corresponds to the higher wavelength value for this mode index , that is to say to n = 1). We can observe that the approximation is very good.

However the quality of the approximation decreases very quickly when the mode index n increases. In Fig. A2 we have depicted the graphs of function f over [0, kL] for L = 1.2 R and ψ 2 over [0, a1] for = 377 and λ = 1548.20 nm (it corresponds to the second higher wavelength value for this mode index , that is to say to n = 2) and for = 368 and λ = 1550.28 nm (it corresponds to the 5th higher wavelength value for this mode index , that is to say to n = 5).

Thus, the approximation (A.7) works well for the higher wavelength resonance value (mode index n = 1) but the error increases quickly when n becomes higher due to the oscillations of the Bessel functions ψ after its first zero and this behavior is not taken into account.

Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr The main advantage of the approximation formula (A.7) is that the right-hand side integral can be computed exactly as follows. First, using formula 5.54.2 p. 629 of [START_REF] Gradstein | Table of Integrals, Series, and Products. Table of Integrals, Series, and Products Series[END_REF] we have Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

Figure 1 .

 1 Figure 1. Picture of an optical micro-sphere made at CNRS FOTON Laboratory.

  0. It follows that the complex vector fields E and B satisfy the vectorial Helmholtz equation in each domain of constant electric permittivity ∆E(x) + ω 2 εµ0E(x) = 0 and ∆B(x) + ω 2 εµ0B(x) = 0 (2.10) Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr together with the conditions div E(x) = 0 and div B(x) = 0. (2.11)

Figure 2 .

 2 Figure 2. Definition of the spherical coordinate system used.

  According to (3.33) and (3.35), the electromagnetic field (E, B) for TM modes is given by

. 60 )

 60 Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr Condition (3.45c) is identically satisfied when relations (3.57) and (3.59) are taken into account. Thus, boundary conditions for TM modes give rise to the following system of equations with unknowns A TM i and A TM o to be fulfilled

  .63) The value of A TM i for a given experiment can be computed from the condition that the total power in the device is imposed or known. Equation (3.62) is referred as the modal equation for TM whispering gallery modes. Remark 1 A TE mode (resp. TM mode) of the sphere is obtained by solving the modal equation (3.53) (resp. (3.62)). As a consequence, a mode is described in terms of three integers : the integer involved in the modal equation, one integer n used to label the solutions of the modal equation for a fixed value of and the integer m involved in the expression of the electromagnetic field of the mode, see (3.46)-(3.47) and (3.55)-(3.56).

  k0R)(3.64) with P = k/k0 = N for TE modes and P = k0/k = 1/N for TM modes where N is the optical index of the micro-sphere, k0 = 2π/λ where λ is the wavelength, k = k0N and R is the micro-sphere radius. In equation (3.65), ψ and ζ denote Riccati-Bessel functions of the first kind and of the third kind respectively.Note that we can either solve equation (3.64) for determining the mode numbers for a given wavelength λ or conversely, solve equation (3.64) for determining the wavelengths λ for a given mode number .

+ 1 2

 1 .69) Taking into account relations (3.66) to (3.69), the modal equation (3.65) can be reformulated asY functions satisfy for all ν ∈ R the recurrence relation[START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF][START_REF] Gradstein | Table of Integrals, Series, and Products. Table of Integrals, Series, and Products Series[END_REF]

2 (

 2 kR) through the recurrence relation (3.76).

Figure 3 .

 3 Figure 3. Radial function dependency of the TE modes for the wavelength 844.4558 nm (n = 1) on the left and 697.3817 nm (n = 5) on the right.

( 3 .

 3 71) for λ = 810 nm in order to determine the values of for which the resonance occurs. Equation (3.71) is numerically solved with Matlab using the fsolve command. The bounds for the modes number given by (3.73) indicate that the area of interest in units of nm is the interval [194.37, 282.65]. In Fig. 4 is depicted the modal function (3.72) for TE modes (where P = N ) on the interval [0, 300]. It reads

Figure 6 .

 6 Figure 6. Radial function dependency of the TE modes for the 1th, 2nd, 3rd, 10th resonance wavelength (from up left to right down).

Figure 7 .

 7 Figure 7. Norm of the electric field inside the sphere in the plane of azimuthal angle ϕ = 0 for TE modes with degree and order ( , m) = (271, 271),( , m) = (271, 270) and ( , m) = (271, 10).

Figure 8 .

 8 Figure 8. Volume in µ-m of the WG modes with index = 271 as a function of index m ∈ [0, ].

  SLEGEND associated Legendre function with Schmidt semi-normalization • SBESSELJ spherical Bessel function of the first kind • SBESSELY spherical Bessel function of the second kind • SBESSELH spherical Bessel function of the third kind (Hankel function) • DSLEGEND derivative of the associated Legendre function with Schmidt seminormalization • DSBESSELJ derivative of the spherical Bessel function of the first kind • DSBESSELY derivative of the spherical Bessel function of the second kind • DSBESSELH derivative of the spherical Bessel function of the third kind (Hankel function) • HANSENM Hansen solution M m of the spherical vectorial wave equation • HANSENN Hansen solution N m of the spherical vectorial wave equation Last, the WGMode Matlab Toolbox contains low level internal functions: • MODEQEL Modale equation for TE or TM modes for a fixed wavelength and a variable mode index • MODEQVW Modale equation for TE or TM modes for a fixed mode index and a variable wavelength • PLOTD plot a discontinuous function removing the draw line at function discontinuity jumps • VSHNOR computes the square of the euclidean norm of the vector spherical harmonics

  Micro-sphere radius [micro-m] = 100 Optical index of the micro-sphere = 1.453 Mode index ell = 602 TE or TM mode ? (TE/TM) : TE

  Draw the radial behavior of the mode for some of the computed wavelength? (y/n) : y Wavelength number of the mode: 1 Draw for another wavelength? (y/n) : y Wavelength number of the mode: 20 Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr Draw for another wavelength? (y/n) for ell = 602 and λ [nm] = 1479.for ell = 602 and λ [nm] = 1220.583

5 |E|

 5 >> volmodMicro-sphere radius [micro-m] = 100 Optical index of the micro-sphere = 1.453 Wavelength [nm] = 1479.8985 TE or TM mode ? (TE/TM) : TE Mode number ell = 602 Mode number m = 600 Mode volume (in micron.m^3)= 12589.4749 >> Last, the script entitled pltmod draws the norm of both the electric field and magnetic induction for a mode defined by the indexes and m in a plane given in spherical coordinates by either a user defined azimuth ϕ or a user defined inclination angle θ, see Figure 2 on p.8. >> pltmod Micro-sphere radius [micro-m] R = 25 Optical index of the micro-sphere = 1.453 TE or TM mode : TE Wavelength [nm] = 804.57 Mode number ell = 271 Mode number m = 271 Plot required in an azimuthal plane phi=cste (A) or in a polar plane theta=cste (P) = A Azimuthal angle [radian] : phi = 0 Minimum radial distance (% R) = 0 Maximal radial distance (% R) = 1 Minimum polar angle [radian] = 0 Maximal polar angle [radian] = pi >> Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr z = r cos(θ) x 10 -5 x = r sin(θ) x 10 -

5 |E|

 5 Micro-sphere radius [micro-m] R = 25 Optical index of the micro-sphere = 1.453 TE or TM mode : TE Wavelength [nm] = 804.57 Mode number ell = 271 Mode number m = 271 Plot required in an azimuthal plane phi=cste (A) or in a polar plane theta=cste (P) = A Azimuthal angle [radian] : phi = 0 Minimum radial distance (% R) = 0.8 Maximal radial distance (% R) = 1 Minimum polar angle [radian] = pi/2-pi/8 Maximal polar angle [radian] = pi/2+pi/8 >> z = r cos(θ) x 10 -5 x = r sin(θ) x 10 -radius [micro-m] R = 25 Optical index of the micro-sphere = 1.453 TE or TM mode : TE Wavelength [nm] = 804.57 Mode number ell = 271 Mode number m = 271 Plot required in an azimuthal plane phi=cste (A) or in a polar plane theta=cste (P) = P Polar angle [radian] : theta = pi/2 Minimum radial distance (% R) = 0 Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr Maximal radial distance (% R)

ψ 2

 2 (s) if s ∈ [0, kR] k0 k 2 ψ (kR) 2 χ (k0R) 2 χ 2 (s/N ) if s ∈ [kR, kL]since according to (3.54) we have

Figure A1 .

 A1 Figure A1. Left: graphs of function f over [0, kL] for L = 1.2 R (red line) and ψ 2 over [0, a 1 ] (black line) where a 1 denotes the first zero of ψ . Right: zoom on the area of interest.

Figure A2 .

 A2 Figure A2. Graphs of function f over [0, kL] for L = 1.2 R (red line) and ψ 2 over [0, a 1 ] (black line) where a 1 denotes the first zero of ψ . Left: for = 377, λ = 1548.20 nm and n = 2. Right: for = 368, λ = 1550.28 nm and n = 5.

2 ψ 2 ( 1 0ψ 2

 2212 (s) 2 -ψ -1 (s) ψ +1 (s) . (A.9) Since ψ (a1) = J + 1 a1) = 0, we deduce that a (s) ds = -a1 2 ψ -1 (a1) ψ +1 (a1).

  7 . It follows that k0R = 387.850 and kR = 562.383.

	Moreover,
	Y194.5(387.850) ≈ -0.021806760787795,	and J194.5(387.850) ≈ 0.037705042371296
	Therefore the evaluation of the ratios Y -1 2	(k0R)/Y + 1 2	(k0R) and J -1 2	(kR)/J + 1 2

  p. 34.

	x 10 -3	Radial function dependency for ell = 100 and λ [nm] = 844.4558			Radial function dependency for ell = 100 and λ [nm] = 697.3817
								0.015
								0.01
								0.005
	8							0
	6
								-0.005
	4
	2							-0.01
	0	0.2	0.4	0.6	0.8	1	1.2	1.4	0	0.2	0.4	0.6	0.8	1	1.2	1.4
				r/R								r/R
		Technical Report -CNRS UMR 6082 FOTON -http://foton.cnrs.fr

Table 1 .

 1 Altogether we have computed 34 possible values for the mode number solution to the modal equation (3.71) for λ = 810 nm. They are given in table 2. There are only 15 of them that satisfy the bound conditions for given by 194.3741 < < 282.6521, see(3.73). None of them corresponds to an integer value. However, we obtained that for the wavelength λ = 810 nm, the maximal value for can be round up to 271. Resonance frequencies in units nm for = 100. Modal function F λ as defined in (3.77) for the TE modes on the interval [0, 300] on the left and on interval[194.37, 282.65] on the right.

	n	λ	n	λ
	1 844.4558 2 794.0584
	3 755.9579 4 724.4811
	5 697.3817 6 673.5653
	7 652.2256 8 632.441
	Technical Report -CNRS UMR 6082 FOTON -http://foton.cnrs.fr

.77)

We have also provided a zoom on the area of interest

[194.37, 282.65

].

Table 2 .

 2 The 34 possible values for the mode numbers .

  is depicted the modal function F as given by (3.72) for = 271 in the interval [580.5908, 844.9405]. Altogether we have computed 21 possible values for the wavelength λ solution to the modal equation (3.71) for = 271. They are given in table3. We found that, around 810 nm, for a mode number = 271, the resonance occurs for a wavelength of 808.6189 nm. Figure 5. Modal function F as defined in (3.72) for the TE modes on the interval [580.5908, 844.9405].

				Modale function for TE modes
		2
		1.5
		1
		0.5
		0
		-0.5
		-1
		-1.5
		550 -2	600	650	700 λ [nm]	750	800	850
	n	λ	n	λ	n	λ		n	λ
	1 808.6189 2 782.8448 3 762.5956 4 745.3336
	5 730.0277 6 716.1493 7 703.4341 8 691.5913
	9 680.5414 10 670.1525 11 660.2922 12 650.9607
	13 642.0786 14 633.5929 15 625.451 16 617.6526
	17 610.1715 18 602.9812 19 596.0553 20 589.3408
	21 582.8113
	Technical Report -CNRS UMR 6082 FOTON -http://foton.cnrs.fr

Table 3 .

 3 The 21 possible values for the wavelength λ (in units of nm) for which the resonance occurs for the mode number = 271.

  .30) 4.2.2. Computation of the maximum value of the electric field inside the micro-sphere It remains to compute E 2

  Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

	>> ellres
	Micro-sphere radius [micro-m] = 100
	Optical index of the micro-sphere = 1.453
	Wavelength [nm] = 1480
	TE or TM mode ? (TE/TM) : TE
	Modes number ell :	17	494.3939
	1	425.8768	18	499.3159
	2	429.7629	19	504.3468
	3	433.6548	20	509.4962
	4	437.5904	21	514.7759
	5	441.5826	22	520.1993
	6	445.6311	23	525.783
	7	449.7355	24	531.5469
	8	453.8967	25	537.5159
	9	458.1165	26	543.7219
	10	462.3977	27	550.2064
	11	466.7431	28	557.0256
	12	471.1563	29	564.2598
	13	475.6411	30	572.0302
	14	480.2019	31	580.5386
	15	484.8438	32	590.1758
	16	489.5723	33	601.9616

Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

Appendix A. F. Treussart approximation formula for volume of modes with large mode number

In this appendix we detail the way to obtain F. Treussart approximation formula for mode volume (formula 1.58 p. 51 of his Phd thesis [START_REF] Treussart | Étude expérimentale de l'effet laser dans des micro-sphères de silice dopées avec des ions néodyme[END_REF]) and we analyze and quantify the various approximations leading to this formula.

According to [START_REF] Treussart | Étude expérimentale de l'effet laser dans des micro-sphères de silice dopées avec des ions néodyme[END_REF], the volume of a whispering gallery mode in a micro-sphere is defined as the integral over the whole space of the energy density normalized by its maximum value inside the micro-sphere, i.e.

where w denotes the energy density given as a function of the position vector x = (x1, x2, x3) by

and wmax denotes the maximum value of the energy density inside the micro-sphere. In relation (A.2), the quantity E * (resp. B * ) stands for the adjoint (conjugate transpose) of E (resp. B).

When we assume that the energy losses by diffraction and diffusion can be neglected, the energy conservation law implies that the contribution of the electric field and the magnetic field to the energy density in (A.1) are equal [START_REF] Jackson | Classical Electrodynamics Third Edition[END_REF]. As a consequence,

Changing the normalization constant wmax to max x∈S(0,R)

Treussart obtains the following expression for the mode volume:

where Emax denotes the maximum value of the euclidean norm of the electric field.

In his Phd thesis [START_REF] Treussart | Étude expérimentale de l'effet laser dans des micro-sphères de silice dopées avec des ions néodyme[END_REF], F. Treussart only deals with the case of TE modes. For TE modes the electric field is given in spherical coordinates by

where X m denotes the vector spherical harmonics and ψ and χ denote respectively the Ricatti-Bessel functions of first and second types. The integral (A.4) in spherical coordinates then reads

where L is a positive number large enough so that the electromagnetic field can be neglected at a radial distance greater than L.

Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

We can simplify the result as follows. We have (see [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF] formula 10.1.21 p. 439)

and

Since ψ (a1) = 0 by product of the 2 above identities we deduce that

and we conclude that

The result can also be expressed in terms of the spherical Bessel function j as follows

We now have to compute E 2 max where Emax denotes the maximum value of the electric field inside the micro-sphere. We have

We first compute the maximum over [0, R] of |j (kr)| = |ψ (kr)/kr|. It follows from the known behavior of the spherical Bessel function, see e.g. [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF][START_REF] Gradstein | Table of Integrals, Series, and Products. Table of Integrals, Series, and Products Series[END_REF], that the global maximum of |j | coincides with the first local maximum of j . Thus, we are looking for the first zero a 1 of the derivative of the spherical Bessel function of order . An approximation based on an asymptotic expansion for the first zero (denoted a ,1 ) of j for large is given in [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF] (see formula (10.1.59) p. 441):

Actually, a ,1 can be considered as an approximation for the first zero of j even for small . For instance for = 5 we have a ,1 = 6.7606 whereas an accurate value for the first zero of j computed with the symbolic computation software Maple is found to be 6.7564. Thus,

Let us now consider the computation of sup

From (3.27) and (3.7), we have

Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

where P m denotes the Associated Legendre function of degree and order m and the constant C m is defined in (3.8), namely

Let us consider the special case when m = . We have

so that

We deduce that

Finally, we conclude that for m = the volume of the TE mode can be approach by

For large values of we can use Stirling's formula to approach the factorial terms. We obtain

F. Treussart formula is deduced by using the relation k = 2πN/λ and for large values of the approximation ( + 1)

Appendix B. Mode volume formula in lossless media

The volume of a whispering gallery mode in a micro-sphere is defined, see [START_REF] Kippenberg | Nonlinear Optics in Ultra-high-Q Whispering-Gallery Optical Microcavities[END_REF][START_REF] Treussart | Étude expérimentale de l'effet laser dans des micro-sphères de silice dopées avec des ions néodyme[END_REF], as the integral over the whole space of the energy density normalized by its maximum value inside the micro-sphere, i.e.

Technical Report -CNRS UMR 6082 FOTONhttp://foton.cnrs.fr

where w denotes the energy density given as a function of the position vector x = (x1, x2, x3) by

and wmax denotes the maximum value of the energy density inside the micro-sphere. In relation (B.2), the quantity E * (resp. B * ) stands for the adjoint (conjugate transpose) of E (resp. B) so that

In this appendix, we show that when we assume that the energy losses by diffraction and diffusion can be neglected, the energy conservation law implies that the contribution of the electric field and the magnetic field to the energy density in (B.1) are equal, i.e.

To prove this result, we use the following expression deduced from harmonic Maxwell's equation (2.9d)

It follows from (B.3) that

where N is the micro-sphere optical index, B(0, R) denotes the ball with radius R and center the origin and B (0, R) its complement in R 3 . Then using Green's formula for the curl operator in each domain B(0, R) and B (0, R), we get

where S(0, R) denotes the sphere with radius R and center the origin, n its unit outward normal vector and the brackets [ ] mean that we have to take into account the jump across the surface of the quantity inside the brackets. From harmonic Maxwell's equations (2.9a) and (2.9d), we deduce that

Therefore,