Charlotte Herzog

Ombeline De

La Rochefoucauld

Fabrice Harms

Philippe Zeitoun

Xavier Granier

Unifying the refocusing algorithms and parameterizations for traditional and focused plenoptic cameras

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

A plenoptic camera allows to acquire and separate spatial and directional information of the light coming from a scene [1,2]. It allows applications such as refocusing at different depths from the one where the image has been acquired. In the literature, different refocusing algorithms are presented for several optical plenoptic configurations [3 -6]. We have previously shown the continuity between these optical designs, and the similarities and differences between the associated algorithms [START_REF] Herzog | Proc. SPIE 10677, Unconventional Optical Imaging[END_REF]. Here we propose a unique parameterization of the light rays in a plenoptic setup, allowing the development of a unique refocusing algorithm valid for any plenoptic configurations, based on the integration of étendues' intersections in the image-space. With this method, we aim at refocusing images at any distances from the camera, without the discontinuity due to the change of optical configurations.

Unifying the refocusing algorithms and parameterizations for traditional and focused plenoptic cameras

Charlotte Herzog 1,2,3 , Xavier Granier 2,3,4 , Fabrice Harms 1 , Philippe Zeitoun 5 & Ombeline de La Rochefoucauld 1

2 .LOA 1 .

 21 Parametrization of the Light-Field Conclusions • Unification of the two historical algorithms. • Algorithm based on the intersection of étendues instead of using pixels as points. • First implementation of algorithm → validation of the idea. • Future work: improving execution time. VOXEL project EU's Horizon 2020 nº 665207 Validity domains of historical refocusing algorithms Bijection between object-space and image-space for each microlens k For each pixel 𝑿 in the image to reconstruct & For each microlens k: • Transform of the étendue defined by pixel X and main lens aperture U : For each pixel i on the sensor: • étendues of intersection 𝑄 𝑅 & 𝐾 𝐼 defined by microlens k and pixel i • add the weighted contribution of pixel i to reconstructed pixel X Quality of reconstructed images depends on both the optical design and the algorithm used.

1

  Imagine Optic, Orsay/Bordeaux, France; 2 Laboratoire Photonique, Numérique et Nanosciences, Talence, France 3 Inria Bordeaux Sud-Ouest; 4 UMS Archéovision 3D SHS, Pessac, France; 5 Laboratoire d'Optique Appliquée, UMR7639, Palaiseau, France

Intersection of étendues in the image space

Integration is approximated on a reduced set of samples


Stéphane Balac
email: stephane.balac@univ-rennes1.fr

WGMode : A Matlab toolbox for whispering gallery modes volume computation in spherical optical micro-resonators

Keywords: optical micro-resonators, whispering gallery modes, mode volume

We present a Matlab toolbox aimed at studying Whispering Gallery Modes (WGM) in spherical optical micro-resonators. Whispering gallery modes are specific resonances of an optical wave trapped inside a dielectric cavity with smooth edges due to continuous total internal reflection of light. The novelty of the work stands in the implementation of analytical formulas to compute the volume of a WGM in a spherical cavity while only approximate asymptotic expressions or purely numerical computations based on the Finite Element Method were previously available. The toolbox also provides Matlab programs to compute and visualize WGM.

from Maxwell's equations using Hansen's method [2]. The expression of the mode field involves Vector Spherical Harmonics and requires the finding of roots of a non-linear equation with Bessel's functions. From the analytical expression of the mode field, the mode volume integral can be computed analytically. These formulas are implemented in the WGMode toolbox to compute the volume of a WGM.

Additional comments including Restrictions and Unusual features:

The WGMode toolbox has been tested with Matlab version 8.1 (R2013a), version 9 (R2016a), version 9.4 (R2018a). The WGMode toolbox also works under GNU Octave software (tested with version 4.0.0). However, under GNU Octave the scripts pltmod and volmod written to display the mode field (electric field and magnetic induction parts) and to compute the mode volume provide the results after an excessive time due to the computation of the Vector Spherical Harmonics that is not efficient under Octave.

Introduction

Optical Whispering gallery modes (WGM) are specific resonances of electromagnetic waves inside a dielectric cavity (the resonator) with smooth surface. They correspond to light-waves circling around the cavity, almost perfectly guided round by optical total internal reflection, that meet a 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). The unique combination of strong temporal and spatial confinement of light in an optical micro-resonator makes these systems of interest for various applications in optics such as lasing, frequency comb generation, sensing, etc. For a review on micro-resonators' fundamentals, fabrication processes and applications, we refer to [1,2,[START_REF] Benson | Frontiers in Planar Lightwave Circuit Technology[END_REF][START_REF] Ward | WGM microresonators: sensing, lasing and fundamental optics with microspheres[END_REF].

Among the various shapes for the dielectric cavity (disks, rings, toroids, photonic crystal cavities, cylindrical optical fibers, etc.), although less easily integrable, spherical micro-resonators have the advantage of being the simplest three-dimensional WGM resonator and they can exhibit very high quality factors. Spherical micro-resonators are typically ten to hundred micrometers in diameter, and they are often fabricated in silica glass by simply melting the tip of a conventional single-mode fiber [2,[START_REF] Benson | Frontiers in Planar Lightwave Circuit Technology[END_REF][START_REF] Ward | WGM microresonators: sensing, lasing and fundamental optics with microspheres[END_REF][START_REF] Righini | Whispering gallery mode microresonators: Fundamentals and applications[END_REF].

Properties of WGM in an optical microcavity are well described by the quality factor Q of the mode that accounts for energy losses and is related to the temporal confinement of light and by its mode volume V that accounts for spatial confinement of light. Depending on the shape of the dielectric cavity, numerical computation of these two quantities can be complicated. Spherical cavity is the only rare exception.

A part of optical losses in an optical micro-resonator is caused by the fact that, unlike a flat surface, the total internal reflection from a curved surface does not exist and leads to a radiation of the wave from the dielectric cavity [START_REF] Righini | Whispering gallery mode microresonators: Fundamentals and applications[END_REF][START_REF] Talebi | Analytical modeling of quality factor for shell type microsphere resonators[END_REF][START_REF] Borselli | Beyond the Rayleigh scattering limit in high-Q silicon microdisks: Theory and experiment[END_REF]. This phenomenon is taken into account through the radiative quality factor which is defined as the ratio of the real part of the wave number of a given WGM and its imaginary part. Other energy losses in a micro-resonator include surface scattering losses, surface absorption losses, material absorption losses and are taken into account in the definition of the intrinsic quality factor. These additional contributions to the quality factor are difficult to handle mathematically and to compute analytically [START_REF] Talebi | Analytical modeling of quality factor for shell type microsphere resonators[END_REF][START_REF] Borselli | Beyond the Rayleigh scattering limit in high-Q silicon microdisks: Theory and experiment[END_REF][START_REF] Oraevsky | Review: Whispering-gallery waves[END_REF]. For this reason, in the WGMode toolbox only the radiative quality factor is computed.

The WGMode toolbox was primarily written to compute the mode volume V of a WGM in a spherical micro-resonator. The mode volume is defined as the ratio of the total energy stored inside the resonator for that mode and the maximum energy density of the mode [9, Sec. 1.6], [10, Sec. 2.7]:

V = 1 M R 3 ε(r) 2 E(r) 2 + 1 2µ(r) B(r) 2 dr (1)
where E and B are respectively the electric field and magnetic induction of the mode, µ and ε are respectively the magnetic permeability and electric permittivity, M = max 2 is the maximum of the energy density of the optical mode and denotes the euclidean norm in C 3 . Note that the definition (1) of the mode volume is open to question since integrating over the entire space leads to a diverging integral since modes of an optical micro-resonator grow exponentially at infinity (see below the expression (3)-( 4) of a TE mode). This issue is mostly ignored in the literature and overcome by integrating over a finite integration domain. Such an approach can result in relatively small errors when dealing with modes of high quality factor. We refer to [START_REF] Muljarov | Exact mode volume and Purcell factor of open optical systems[END_REF] for a detailed discussion of the issue and nevertheless we adopt formula (1) as expression of the mode volume.

r∈R 3 ε(r) 2 E(r) 2 + 1 2µ(r) B(r)
Mode volume provides relevant information on WGM features. It is related to the free spectral range (FSR): lower mode volume means higher FSR. For lasing applications of microresonators it is also related to the lasing characteristics since stimulated as well as spontaneous emission rates are inversely proportional to the mode volume whereas the lasing threshold is directly proportional to the mode volume [START_REF] Wu | High-q, low-mode-volume microsphere-integrated fabry-perot cavity for optofluidic lasing applications[END_REF]. For sensing applications of micro-resonators, the mode volume is a parameter connected with the coupling of the resonator and related to its sensitivity to external parameters [START_REF] Barucci | Optical microbubble resonators with high refractive index inner coating for bio-sensing applications: An analytical approach[END_REF].

The novelty of the WGMode toolbox is to provide tools to compute the volume of a WGM in a spherical optical micro-resonator based on an analytical approach to evaluate the volume integral in (1). This ensures both accurate results and fast computation times. Alternative approaches investigated in the literature to compute the volume of a WGM include approximate formulas based on asymptotic expansions of the mode field [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] Kippenberg | Nonlinear optics in ultra-high-Q whispering-gallery optical microcavities[END_REF][START_REF] Oraevsky | Review: Whispering-gallery waves[END_REF] and purely numerical computations based on the Finite Element Method (FEM) [START_REF] Srinivasan | Cavity Q, mode volume, and lasing threshold in small diameter AlGaAs microdisks with embedded quantum dots[END_REF][START_REF] Oxborrow | Traceable 2-D Finite-Element simulation of the whisperinggallery modes of axisymmetric electromagnetic resonators[END_REF]. The advantage of mode volume formulas based on asymptotic expansions of the WGM field is their simplicity. However, they are based on strong assumptions that limit their domain of validity to a restricted range of possible values for the mode indexes and geometric size of the cavity. Opposite is the use of the FEM to compute the WGM field and mode volume. The FEM does not introduce geometric limitations or assumptions on the mode indexes of the WGM. However, computation can take a long time compared to the use of analytical formulas, but more importantly it is not at all obvious to determine the mode indexes of the WGM from a FEM computation which can be considered as a severe disadvantage for the study of spherical optical micro-resonators.

The paper is organized as follows. In section 2, we recall the expression of WGM field in a spherical optical micro-resonator as it can be found e.g. in [START_REF] Oraevsky | Review: Whispering-gallery waves[END_REF][START_REF] Righini | Whispering gallery mode microresonators: Fundamentals and applications[END_REF]. In section 3, we show how the mode volume integral (1) can be computed analytically starting from the expression of the WGM field given in Section 2. Section 4 is devoted to an overview of the WGMode toolbox, including a comprehensive test run. Finally, we conclude in Section 5 by presenting benchmark results comparing the three approaches mentioned above for computing the mode volume of a WGM in a spherical micro-resonator.

Expression of WGM in a spherical micro-resonator

For a spherical dielectric cavity, an analytical expression of the electromagnetic field

E(r, t) = Re E(r) exp(iωt) , B(r, t) = Re B(r) exp(iωt) (2)
corresponding to a WGM (ω denotes the resonance frequency) is obtained from Maxwell's harmonic equations using the so-called Hansen's method [START_REF] Stratton | Electromagnetic Theory[END_REF]. We denote by N the optical index of the dielectric cavity and by R its radius. We set k 0 = ω/c and k = k 0 N . Using Hansen's method, it is found that there exist two kinds of electromagnetic fields with different polarization termed TE (transverse electric) and TM (transverse magnetic) defined in spherical coordinates (r, θ, ϕ), see Fig. 1, as follows [START_REF] Oraevsky | Review: Whispering-gallery waves[END_REF][START_REF] Righini | Whispering gallery mode microresonators: Fundamentals and applications[END_REF]. For TE modes, the mode field can be expressed as

E(r, θ, ϕ) =        ψ (kr) kr X m (θ, ϕ) if r ≤ R k 0 k ψ (kR) ζ (k 0 R) ζ (k 0 r) k 0 r X m (θ, ϕ) if r > R (3)
and

B(r, θ, ϕ) =                        ik ω ( + 1) ψ (kr) k 2 r 2 Z m (θ, ϕ) + ψ (kr) kr Y m (θ, ϕ) if r ≤ R ik 2 0 k ω ψ (kR) ζ (k 0 R) ( + 1) ζ (k 0 r) k 2 0 r 2 Z m (θ, ϕ) + ζ (k 0 r) k 0 r Y m (θ, ϕ) if r > R (4)
where , m are integers such that -≤ m ≤ , X m , Y m and Z m denote the Vector Spherical Harmonics defined from the Scalar Spherical Harmonics Y m as

Z m = Y m e r , Y m = r ∇Y m , X m = ∇Y m ∧ r (5)
and ψ and ζ denote respectively the Riccati-Bessel functions of first and third kinds. They are related to the Bessel function J + 1 2 and Hankel function H

respectively by [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF][START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF]

ψ (x) = πx 2 J + 1 2 (x), ζ (x) = πx 2 H (2) + 1 2 (x). (6)
The Spherical Surface Harmonics of degree and order m is defined as, see [START_REF] Atkinson | Spherical Harmonics and Approximations on the Unit Sphere: An Introduction[END_REF],

Y m (θ, ϕ) = C m P m (cos(θ)) e imϕ -≤ m ≤ , ∈ N (7)
where the normalization constant C m is taken to be

C m = (2 + 1) 4π ( -m)! ( + m)! . (8)
Note that the mode field (E, B) given by ( 3)-( 4) is defined up to a normalization constant. Note also that Hankel function

H (2) + 1 2
, and its Riccati-Bessel counterpart ζ used in (3)-( 4), express outward-propagating spherical-wave solutions of the spherical wave equation for the sign convention +iωt used here for the frequency in (2).

For TM modes, we have and

E(r, θ, ϕ) =                ( + 1) ψ (kr) k 2 r 2 Z m (θ, ϕ) + ψ (kr) kr Y m (θ, ϕ) if r < R ψ (kR) ζ (k 0 R) ( + 1) ζ (k 0 r) k 2 0 r 2 Z m (θ, ϕ) + ζ (k 0 r) k 0 r Y m (θ, ϕ) if r > R (9)
B(r, θ, ϕ) =        ik ω ψ (kr) kr X m (θ, ϕ) if r ≤ R ik 0 ω ψ (kR) ζ (k 0 R) ζ (k 0 r) k 0 r X m (θ, ϕ) if r > R (10)
Here again, the mode field (E, B) given by ( 9)-( 10) is defined up to a normalization constant.

Resonance frequencies ω are obtained by solving the so-called modal equation [START_REF] Oraevsky | Review: Whispering-gallery waves[END_REF][START_REF] Righini | Whispering gallery mode microresonators: Fundamentals and applications[END_REF] P ψ (kR)

ψ (kR) = ζ (k 0 R) ζ (k 0 R) (11)
where P = N for TE modes and P = 1/N for TM modes. As a consequence, a mode is described in terms of three integers: the integer involved in the modal equation [START_REF] Muljarov | Exact mode volume and Purcell factor of open optical systems[END_REF], one integer n used to label the solutions of the modal equation [START_REF] Muljarov | Exact mode volume and Purcell factor of open optical systems[END_REF] for a fixed value of and the integer m involved in the expression of the electromagnetic field of the mode, see (3)-( 4) and ( 9)- [START_REF] Kippenberg | Nonlinear optics in ultra-high-Q whispering-gallery optical microcavities[END_REF]. Note that this underscores a mode degeneracy since there exist 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 polar 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 e r . 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 ϕ. Negative values of m signify counter-propagating modes. 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 ( 11) has a solution, is termed the fundamental mode. WGM refers in the literature to modes with polar index close to the highest possible value for which the modal equation ( 11) has a solution, together with a small difference -|m| and n close to unity [START_REF] Oraevsky | Review: Whispering-gallery waves[END_REF].

Frequencies ω obtained as solutions to the modal equation ( 11) are complex valued. It is usual in the study of WGM in optical micro-resonators to introduce the radiative quality factor as the ratio of the real part of the resonance frequency ω and its imaginary part:

Q = Re(ω) Im(ω) . ( 12
)
The smaller losses by radiation are, the smaller Im(ω) is and therefore the higher is the radiative quality factor. We wish to emphasize that the quality factor of an optical micro-resonator usually takes into account other losses phenomena, that can be predominant [START_REF] Talebi | Analytical modeling of quality factor for shell type microsphere resonators[END_REF][START_REF] Borselli | Beyond the Rayleigh scattering limit in high-Q silicon microdisks: Theory and experiment[END_REF][START_REF] Oraevsky | Review: Whispering-gallery waves[END_REF]. Therefore, it is likely to happen that the radiative quality factor is not close to the true, experimentally measured, quality factor.

In practice, spherical optical micro-resonators often correspond to resonators with a large radius compared to the resonance wavelength λ = ω/c of interest in the optical domain. For such resonators, the solutions to the modal equation [START_REF] Muljarov | Exact mode volume and Purcell factor of open optical systems[END_REF] corresponding to WGM are related to large values of the polar mode index . For example, if we consider a 50 µm radius silica sphere at a wavelength λ ≈ 1.5 µm, we obtain the highest values of greater than 300, see Section 4. As a consequence, for such a mode the radiative leakage of the energy is small and a good approximation consists in assuming the radiative part outside the sphere as negligible. From a mathematical point of view, this approximation consists in considering the following asymptotic expansions of Bessel's functions for large order , see e.g. [17, formula 10.19.2] : ζ (z) ≈ iχ (z) where χ denotes Riccati-Bessel function of the second kind with order related to Bessel function

Y + 1 2 by χ (x) = - πx 2 Y + 1 2 (x). ( 13
)
Analytical formulas for modes volume in spherical optical micro-resonators presented in the next section and implemented in the WGMode toolbox were obtained under this assumption, mathematically restrictive but very realistic from a practical viewpoint. This approximation leads to consider the following form for the modal equation

P ψ (kR) ψ (kR) = χ (k 0 R) χ (k 0 R) (14)
that admits real solution k 0 .

Computation of the volume of a whispering gallery mode

The volume of a whispering gallery mode is mathematically defined by formula (1). Here we make the assumption that the dielectric cavity is nonmagnetic, i.e. its magnetic permeability µ is the one of vacuum (µ = µ 0 ). 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 is equal [START_REF] Jackson | Classical Electrodynamics[END_REF]:

R 3 ε(r) 2 E(r) 2 dr = R 3 1 2µ 0 B(r) 2 dr. ( 15
)
As a consequence, the volume of a WGM is most often computed in the literature [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] Kippenberg | Nonlinear optics in ultra-high-Q whispering-gallery optical microcavities[END_REF][START_REF] Ward | WGM microresonators: sensing, lasing and fundamental optics with microspheres[END_REF][START_REF] Oraevsky | Review: Whispering-gallery waves[END_REF][START_REF] Oxborrow | Traceable 2-D Finite-Element simulation of the whisperinggallery modes of axisymmetric electromagnetic resonators[END_REF] from the following formula:

V = 1 ε 0 N 2 E 2 max R 3 ε(r) E(r) 2 dr ( 16
)
where E max denotes the maximum value of the Euclidean norm of the electric field. Formula ( 16) 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 [START_REF] Benson | Frontiers in Planar Lightwave Circuit Technology[END_REF]. It should be noted however that the normalization constant is not the same in (1) and in [START_REF] Stratton | Electromagnetic Theory[END_REF].

For TM modes the mathematical expression of the electric field given by ( 9) is much complicated than for TE modes. That is why computation of the mode volume from ( 16) is not recommended. It is much preferable to use relation [START_REF] Oxborrow | Traceable 2-D Finite-Element simulation of the whisperinggallery modes of axisymmetric electromagnetic resonators[END_REF] to express the volume of a TM mode as

V = c 2 N 2 E 2 max R 3 B(r) 2 dr ( 17
)
where c is the speed of light in free space. Unfortunately, the normalization constant in ( 16) can't be expressed in terms of B and changing it for max r∈R 3 B(r) 2 would not be consistent when comparing volumes of TE and TM modes.

Volume of a TE mode

Because of the evanescent behavior of the electromagnetic field outside the micro-sphere for WGM, we bound the integration domain in ( 16) at a radial distance L > R appropriately selected. For TE modes, the mode volume is therefore given in spherical coordinates as

V = 1 ε 0 N 2 L 0 π 0 2π 0 ε(r) E(r, ϕ, θ) 2 E 2 max r 2 sin(θ) dϕdθdr. ( 18
)
For TE modes, given the approximation of the mode field outside the sphere discussed in p. 7, the electric field is given in spherical coordinates by

E(r, θ, ϕ) =        ψ (kr) kr X m (θ, ϕ) for 0 ≤ r ≤ R ψ (kR) χ (k 0 R) χ (k 0 r) k 0 r X m (θ, ϕ) for R < r ≤ L (19)
Using Chasles's theorem, we split the integral in [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF] into two parts such that the first integral corresponds to the part V i of the mode volume over the cavity and the second integral to the part V o of the mode volume over the exterior domain:

V = V i + V o . We have V i = 1 k 2 E 2 max R 0 π 0 2π 0 ψ 2 (kr) X m (ϕ, θ) 2 sin(θ) dϕdθdr = 1 k 2 E 2 max R 0 ψ 2 (kr) dr π 0 2π 0 X m (ϕ, θ) 2 sin(θ) dϕdθ . (20)
As a consequence of the orthogonality properties of Vector Spherical Harmonics, we have

π 0 2π 0 X m (ϕ, θ) 2 sin(θ) dϕdθ = ( + 1). ( 21
)
Therefore,

V i = ( + 1) k 2 1 E 2 max R 0 ψ 2 (kr) dr = ( + 1) k 1 E 2 max π 2 R 0 r J 2 + 1 2 (kr) dr.
This last integral can be calculated thanks to [21, formula 5.54 p. 629]

x J p (αx) 2 dx = x 2 2 J p (αx) 2 -J p-1 (αx)J p+1 (αx) .
We obtain

V i = ( + 1)πR 2 4k 1 E 2 max J + 1 2 (kR) 2 -J -1 2 (kR)J + 3 2 (kR) . ( 22
)
The integral over the exterior domain reads

V o = 1 k 2 0 N 2 E 2 max ψ (kR) 2 χ (k 0 R) 2 L R π 0 2π 0 χ 2 (k 0 r) X m (ϕ, θ) 2 sin(θ) dϕdθdr = 1 N 2 k 2 0 E 2 max ψ (kR) 2 χ (k 0 R) 2 L R χ 2 (k 0 r) dr π 0 2π 0 X m (ϕ, θ) 2 sin(θ) dϕdθ = ( + 1) N 2 k 0 E 2 max ψ (kR) 2 χ (k 0 R) 2 π 2 L R r Y 2 + 1 2 (k 0 r) dr (23)
The last integral can be evaluated exactly thanks to [21, formula 5.54 p. 629]

x Y p (αx) 2 dx = x 2 2 Y p (αx) 2 -Y p-1 (αx)Y p+1 (αx) .
We obtain

V o = ( + 1)π 4N 2 k 0 E 2 max ψ (kR) 2 χ (k 0 R) 2 L 2 Y + 1 2 (k 0 L) 2 -L 2 Y -1 2 (k 0 L)Y + 3 2 (k 0 L) -R 2 Y + 1 2 (k 0 R) 2 + R 2 Y -1 2 (k 0 R)Y + 3 2 (k 0 R) = ( + 1)π 4N 3 k 0 E 2 max J + 1 2 (kR) 2 Y + 1 2 (k 0 R) 2 L 2 Y + 1 2 (k 0 L) 2 -L 2 Y -1 2 (k 0 L)Y + 3 2 (k 0 L) -R 2 Y + 1 2 (k 0 R) 2 + R 2 Y -1 2 (k 0 R)Y + 3 2 (k 0 R) . (24)
It remains to compute E max , the maximum value of the euclidean norm of the electric field over the dielectric cavity, given by

E max = sup r∈[0,R] ψ (kr) kr × sup (θ,ϕ)∈[0,π]×[0,2π] X m (θ, ϕ) . ( 25
)
We first consider the computation of the maximum of |ψ (kr)/kr| = |j (kr)|, where j denotes the spherical Bessel function of first kind, for r ∈ [0, R]. It is known, see e.g. [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF][START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[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 j where [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF][START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF]

∀x ∈ R j (x) = j -1 (x) -( + 1) j +1 (x) 2 + 1 .
We have j (x) = 0 if and only if

( + 1) J + 3 2 (x) -J -1 2 (x) = 0. ( 26
)
Thus the problem turns out to compute the first positive root of equation ( 26). 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 by [18, formula 10.1.59, 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 . (27)
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 obtained by symbolic computation is 6.7564. Thus, from a computational point of view we use the Matlab command fsolve to solve (26) with an initial guess given by (27). Next, we have to compute the maximum value of X m (θ, ϕ) for all (θ, ϕ) ∈

[0, π] × [0, 2π]. Since X m = ∇Y m ∧ r, we have X m (θ, ϕ) 2 = C 2 ,m m 2 sin 2 (θ) P m (cos(θ)) 2 + sin 2 (θ) (P m ) (cos(θ)) 2
which is a quantity independent of ϕ. From the following relation for the derivative of the associated Legendre function P m [21, formula 8.733 p. 965]

∀x ∈ ] -1, 1[ (P m ) (x) = 1 1 -x 2 P m (x) -( -m + 1) P m +1 (x)
we deduce that

X m (θ, ϕ) 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(θ)) . ( 28
)
Analytical computation of the maximum value of X m (θ, 0) 2 for θ ∈ [0, π] seems out of reach. Therefore, in the WGMode toolbox we proceed by using numerical optimization methods. Unfortunately, expression (28) is not well suited for numerical computation purposes. Indeed for larger than 150 the values of the associated Legendre function P m are 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 functions P m . They are related to each other by the relation [START_REF] Winch | Geomagnetism and Schmidt quasi-normalization[END_REF] ∀x ∈] -

1, 1[ S m (x) =      P (x) if m = 0 (-1) m 2( -m)! ( + m)! P m (x) if m > 0 . (29)
We set D ,0 = 1 and D ,m = (-1) m ( + m)!/2( -m)! for m > 0 so that for all x ∈] -1, 1[ we have P m (x) = D ,m S m (x). In terms of the Schmidt semi-normalized associated Legendre functions, we have

X m (θ, ϕ) 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(θ)) (30)
where

C 2 ,m D 2 ,m =      2 + 1 4π if m = 0 2 + 1 8π if m > 0 .
In the WGMode toolbox, the global maximum of X m (θ, 0) 2 over [0, π] is computed numerically using Matlab fminunc routine. A Trust Region Method algorithm is chosen which requires to provide the gradient of X m (θ, 0) 2 . The initial guess for the optimization algorithm is obtained by comparing the values of X m (θ, 0) 2 over a coarse subdivision of the interval [0, π].

Volume of a TM mode

For TM modes, the magnetic induction is given by

B(r, θ, ϕ) =        i ω ψ (kr) r X m (θ, ϕ) if r < R i ω ψ (kR) χ (k 0 R) χ (k 0 r) r X m (θ, ϕ) if r > R (31)
and the integral [START_REF] Olver | NIST Handbook of Mathematical Functions[END_REF] in spherical coordinates reads

V = c 2 N 2 E 2 max L 0 π 0 2π 0 B(r) 2 r 2 sin(θ) dϕdθdr ( 32
)
where L is a positive number chosen such that the magnetic induction can be neglected at a radial distance greater than L. (This assumption is made possible from the evanescent behavior of the WGM field outside the micro-resonator.) Using Chasles's theorem, we split the integral into two parts such that the first integral holds over the cavity and the second integral holds over the exterior domain. The first integral reads

V i = 1 k 2 E 2 max R 0 π 0 2π 0 ψ 2 (kr) X m (ϕ, θ) 2 sin(θ) dϕdθdr (33)
and we recognize the same integral as in [START_REF] Jackson | Classical Electrodynamics[END_REF]. Therefore, we directly deduce that

V i = ( + 1)πR 2 4k E 2 max J + 1 2 (kR) 2 -J -1 2 (kR)J + 3 2 (kR) . ( 34
)
The integral over the exterior domain reads

V o = 1 k 2 E 2 max L R π 0 2π 0 B(r, θ, ϕ) 2 sin(θ) dϕdθdr = 1 k 2 E 2 max ψ (kR) χ (k 0 R) L R π 0 2π 0 χ 2 (k 0 r) X m (ϕ, θ) 2 sin(θ) dϕdθdr. (35)
Once again, the integral has already been computed, see [START_REF] Delves | A numerical method for locating the zeros of an analytic function[END_REF]. Therefore, we immediately obtain

V o = ( + 1)π 4kE 2 max J + 1 2 (kR) 2 Y + 1 2 (k 0 R) 2 L 2 Y + 1 2 (k 0 L) 2 -L 2 Y -1 2 (k 0 L)Y + 3 2 (k 0 L) -R 2 Y + 1 2 (k 0 R) 2 + R 2 Y -1 2 (k 0 R)Y + 3 2 (k 0 R) . (36)
It remains to compute E 2 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 ( 9) that

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 . ( 37
)
Analytical computation of the maximum value of E 2 max over [0, R] × [0, π] is out of reach. Therefore, in the WGMode 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 (38) -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, π].

Overview of the Matlab Toolbox WGMode

Numerical issues addressed in the WGMode toolbox

Computation of WGM in a spherical optical micro-resonator must address the following issue. Resonance wavelengths (corresponding to the visible spectrum) are obtained by solving the modal equation ( 11) that relies on the polar mode index ∈ N. If one sets any value for the integer , there is a strong possibility that the solutions to the modal equation ( 11) (if they exist) do not correspond to wavelengths in the visible spectrum and even fewer to WGM. In practice, one is interested by WGM around a given wavelength, say λ ≈ 1.5 µm. However, if one sets the value of λ and solves the modal equation [START_REF] Muljarov | Exact mode volume and Purcell factor of open optical systems[END_REF] for the polar mode index there will be little chance that the value of will be an integer. To overcome this issue, we have to proceed in two steps.

1. In a first step, we set the value of the wavelength λ corresponding to a rough approximation of the wavelength of the WGM sought after and we solve the modal equation for the polar mode index . We obtain a set of solutions that are not integer.

2. In a second step, we choose one of the solutions (say the highest for the fundamental mode) and we round it to the nearest integer. We solve once again the modal equation with this integer value of to determine the resonance wavelengths. Because of the continuity of the solutions to the modal equation with respect to the equation's parameters, the first solution (corresponding to a radial mode index n = 1) will be close to the target value.

Another important issue lies in the solving of the modal equation [START_REF] Muljarov | Exact mode volume and Purcell factor of open optical systems[END_REF]. Indeed, because Riccati-Bessel function of the third kind (just as Hankel's functions) is complex valued, solutions to the modal equation ( 11) are complex numbers (the imaginary part of which takes into account losses and is required to compute the Q-factor). From a numerical point of view, it is much more difficult to solve a non-linear equation in the complex plane than in the real line. Therefore, we proceed in two steps. We first solve the approximated modal equation ( 14) on the real line. Then, each real solution of the approximated modal equation ( 14) is used as an initial guess for the iterative method implemented in the Matlab fsolve routine to solve the modal equation [START_REF] Muljarov | Exact mode volume and Purcell factor of open optical systems[END_REF].

For numerical purposes it is more convenient to express the modal equation [START_REF] Srinivasan | Cavity Q, mode volume, and lasing threshold in small diameter AlGaAs microdisks with embedded quantum dots[END_REF] in terms of Bessel functions of first and second kind. We obtain

Y -1 2 (k 0 R) Y + 1 2 (k 0 R) -P J -1 2 (kR) J + 1 2 (kR) = 1 k 0 R - P kR . ( 39
)
For a given mode number , the wavelengths λ for which resonance occurs are obtained by solving equation (39), that is to say by looking for the zeros of the modal function

F : λ ∈ ]0, +∞[ -→ Y -1 2 (k 0 R) Y + 1 2 (k 0 R) -P J -1 2 (kR) J + 1 2 (kR) - 1 k 0 R - P kR . (40)
In the WGMode toolbox, the zeros of this non-linear function on the real line are computed thanks to Matlab's fzero function. It is very useful from a practical point of view to use the following bounds for the solutions corresponding to WGM provided in [9, p. 25]. For a given λ, the mode numbers corresponding to solutions to the modal equation ( 14) are such that

2π(R + δ P ) λ < + 1 2 < N 2π(R + δ P ) λ ( 41
)
where

δ P ≈ λ 2πN P √ N 2 -1 .
There is no solution to [START_REF] Srinivasan | Cavity Q, mode volume, and lasing threshold in small diameter AlGaAs microdisks with embedded quantum dots[END_REF] greater than the upper bound value in (41). Conversely, there are some solutions lower than the lower bound value, but they do not correspond to WGM. Relation (41) implies that for a fixed value of the wavelengths λ for which resonance occurs are such that 2πR

+ 1 2 -α < λ < 2πRN + 1 2 -αN (42)
where α = P/(N √ N 2 -1). The WGMode toolbox uses relations (41) and (42) in the root finding process but one can easily ignore these bounds and provide its own bounds by setting the expert mode variable to the value 1 in the main program. In order to be sure to compute all the solutions to the approximated modal equation ( 14), the interval defined by the above lower and upper bounds is discretized with a very small step-size and sign changes in the values taken by the non-linear function (40) are detected. Then, Matlab's fzero routine is used to accurately compute each solution in the interval it has been located.

To compute the quality factor, one has to solve the exact modal equation ( 11) in the complex plane. This is achieved by using Matlab's fsolve routine with an initial guess for the iterative algorithm used in fsolve given by the (real) obtained by solving the approximated modal equation ( 14) modified by adding a very small imaginary part (Matlab documentation recommends to set a complex initial point when using the non-linear equation solver fsolve in the complex plane). This approach is implemented in the script qfact.m of the WGMode toolbox. Alternatively, since solving the modal equation [START_REF] Muljarov | Exact mode volume and Purcell factor of open optical systems[END_REF] amounts to locate the zeros of an analytical function, this can be achieved by computing by quadrature contour integrals in the complex plane according to the Delves-Lyness Method [START_REF] Delves | A numerical method for locating the zeros of an analytic function[END_REF][START_REF] Ioakimidis | Quadrature methods for the determination of zeros of transcendental functions -a review[END_REF]. This approach is implemented in the script qfactn.m of the WGMode toolbox.

We would like to draw attention to the fact that the solutions to the modal equation ( 11) are complex numbers with a very small imaginary part compared to its real part and that solving equation ( 11) represents a stiff problem from a numerical point of view. The user should pay attention to the information returned by Matlab fsolve routine used to solve equation [START_REF] Muljarov | Exact mode volume and Purcell factor of open optical systems[END_REF] to decide on the reliability of the radiative quality factor value. Moreover, even if computed accurately, this value can differ significantly from the experimentally measured quality factor because of intrinsic loss sources not taken into account in the radiative quality factor but that can be predominant compared to radiative losses.

Particular attention has been paid to the numerical computation of E 2 max the maximum value of the euclidean norm of the electric field by optimization methods. Since optimization methods such as the ones implemented in Matlab fminunc routine can't make a distinction between local maxima and the global maximum and because in general convergence is only ensured when the initial guess is sufficiently close to the global maximum, we first have to compute a reliable initial guess. To this purpose, we compare the values of E 2 max over a coarse grid of the optimization domain (the interval [0, π] for TE modes or the set [0, R] × [0, π] for TM modes). Then, the global maximum of E 2 max is computed by using the fminunc routine for unconstrained optimization problem where a Trust Region algorithm is chosen. The Trust Region method requires that the gradient of E 2 max is provided. Since the fminunc routine could sometimes converge to a local maximum, we compare the maximum provided by fminunc to the value of the objective function with the initial guess. If the two values differ too much, we consider more reliable and safe to use a brute force method to compute the global maximum that consists in comparing the values of E 2 max over a very fine grid of the optimization domain. Finally, to overcome the possible absence of the optimization toolbox in the user Matlab installation, we test if the toolbox is available and if this is not so, we use the above mentioned brute force method to compute the global maximum.

WGMode toolbox contents

The WGMode toolbox is constituted of various Matlab scripts aimed at studying WGM in spherical optical micro-resonators. They allow to explore resonance conditions for TE or TM modes, to visualize WGM in a micro-sphere and to compute the volume and quality factor of any given mode. Namely, the WGMode toolbox contains the following Matlab scripts:

-ELLRES computes, for a given wavelength, the values of the mode index for which a resonance occurs (note that the solutions are real numbers); -WVLRES computes, for a given polar mode index the values of the wavelength for which a resonance occurs;

-VOLMOD computes the volume of a whispering gallery mode in a microsphere;

-QFACT computes for a given whispering gallery mode its quality factor by solving the modal equation by fsolve;

-QFACTN computes for a given whispering gallery mode its quality factor by solving the modal equation by the Delves-Lyness Method;

-PLTMOD plots a whispering gallery mode.

The general organization of the WGMode toolbox is depicted in Fig. 2. The toolbox also provides a collection of special functions:

- Note that an expert mode can be activated by setting the variable expert mode to 1 in the ellres program. It provides a plot 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 , the zero of which are the sought out solutions.

As it can be seen above, the polar mode indexes obtained as solutions to the modal equation [START_REF] Muljarov | Exact mode volume and Purcell factor of open optical systems[END_REF] are not integers. Thus, we have to round up the values to the nearest integer. But by doing so, the couple (λ, ) does not anymore satisfy the modal equation. For each rounded value for , we have to compute the resonance wavelength by solving once again the modal equation. This time, is fixed and we look for λ. Computation is achieved with the Matlab script wvlres. As an example, we consider the first value for (302.418), rounded to 302. An expert mode can be activated by setting the variable expert mode to 1 in the wvlres program. It allows to compare the resonance wavelengths to the one obtained from the Eikonale approach, see e.g. [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], which basically relies on asymptotic expansion for Bessel's functions of large order. It also provides a plot of the modal function F defined in (40).

It is possible thanks to the pltmod script in the WGMode toolbox to draw the WGM field. Three different plots of the euclidean norm of the electric field E(r, θ, ϕ) and magnetic induction B(r, θ, ϕ) are proposed. They correspond to a 1D-plot in the radial direction, to a 2D-plot in an azimuthal plane ϕ = cste and to a 2D-plot in a polar plane θ = cste. We have depicted in Fig. 3 the radial variation of the electric field and magnetic induction for (θ, ϕ) = ( π 2 , 0) corresponding to the TE mode with polar index = 302 and radial index n = 5.

It has a resonance wavelength of 1365.36812 nm. Note that the unit is arbitrary since the mode field is defined up to a multiplicative constant. One can clearly see the five expected maxima in the radial direction for a mode with radial index n = 5. One can also see that both the electric field and magnetic induction are continuous across the sphere border as expected for TE modes. This would not be the case for a TM mode. The pltmod script also allows to draw the WGM field in an azimuthal plane ϕ = cste as shown in Fig. 4 where we have represented the euclidean norm of the electric field and magnetic induction in the plane ϕ = 0 for the WGM with indexes = 302, m = 300, n = 1 (λ = 1502.024 nm) in the case of a spherical cavity with radius R = 50 µm and optical index N = 1.5. Note that the norm of the electric field (resp. magnetic induction) is given in arbitrarily units since the WGM field is defined up to a multiplicative constant. Since a zoom in the area of interest in the Matlab graphical window induces a loss in the image rendering, an expert mode is provided in the pltmod script (obtained by turning the variable expert mode to the value 1 in the script). When this option is chosen it is possible to indicate the bounds of the area to plot (bounds for the polar angle θ and radial distance r). Fig. 4 was obtained using this feature with the following bounds: r

∈ [0.8 R, 1.1 R] and θ ∈ [0.8 π 2 , 1.2 π 2 ]
. One can see in Fig. 4 that for the chosen data, we have -|m| + 1 = 3 which is the number of intensity maxima in the polar direction e θ as stated on p. 6. Similarly, it is possible with the pltmod script to draw the WGM field in an azimuthal plane θ = cste.

With all the information at hand, we are now in position to compute the mode volume. For the dielectric micro-sphere with radius R = 50 µm and optical --------------------------------------------------------------------- ---------------------------------------------------------------------Mode volume (in micron.m^3)= 2259.6416 CPU time for mode volume computations = 1.58 sec.

Always for this particular WGM, using the qfact script, the quality factor is found to be 1.578 10 17 . Note that this quality factor is not in accordance with the quality factor usually observed experimentally because radiative losses are much lower than intrinsic losses for such a micro-resonator. The tolerance values used to end iterations in the fsolve routine are set by default to 10 -9 . These tolerance values can be modified using the expert mode by setting the expert mode variable to the value 1 in the qfact script. Because solving equation [START_REF] Muljarov | Exact mode volume and Purcell factor of open optical systems[END_REF] represents a stiff problem from a numerical point of view, the user should pay attention to the information returned by Matlab fsolve routine to decide on the reliability of the quality factor value. -------------------------------------------------------------------- --------------------------------------------------------------------- One should note that the ratio of the imaginary part and the real part of the resonance is lower than the floating point accuracy and therefore the value returned for the quality factor should be considered with caution. When comparing to the result returned by the Delves-Lyness method, see below, we obtain a very good agreement for the real part of the resonance wavelength but significantly different results for the imaginary part.

Numerical validation

We have compared our computational method for mode volume V in spherical optical micro-resonators to the values obtained when using the Finite Elements Method (FEM) and to the values provided by the following formula [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] based on asymptotic approximation of the WGM field

V ≈ 2π 3 2 λ 2πN 3 z 3 1 2 √ j (z 1 ) j (z 1 ) 2 (43)
where j denotes the spherical Bessel function of the first kind with order , z 1 its first zero and z 1 the first zero of j . This formula is valid for the mode volume of WGM with indexes such that = |m| and n = 1. We refer to [START_REF] Oxborrow | Traceable 2-D Finite-Element simulation of the whisperinggallery modes of axisymmetric electromagnetic resonators[END_REF] for presentation of the way the mode volume of WGM can be computed by the Finite Element method. To perform the computations, we have used the opensource Finite Element software FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF]. We have taken into account the axial symmetry of the micro-sphere to reduce the problem to a vectorial formulation for the magnetic field set in a cylindrical two-dimensional section as explained in [START_REF] Oxborrow | Traceable 2-D Finite-Element simulation of the whisperinggallery modes of axisymmetric electromagnetic resonators[END_REF]. Computations were achieved with Matlab 9.0 on an Intel Core i5 desktop computer (4 Cores, 4 Threads @2.50 GHz) with 8 Go memory. We first consider the case of the micro-sphere of radius 50 µm and optical index 1.5 discussed in section 4. Using the WGMode toolbox, we found that the TE mode with indexes = m = 302, n = 1 has a resonance wavelength λ = 1502.0244 nm and a volume V = 2259.64 µm 3 . Computation of the mode volume took 1.5 sec. Using the approximated formula (43), we obtain a mode volume of 2258.76 µm 3 . The values provided by the FEM are a resonance wavelength of 1502.0195 nm and a mode volume of 2259.79 µm 3 . The FEM computation with Lagrange P 3 Finite Elements and 30474 mesh elements took 324.89 sec. Note that the FEM can compute several modes at once (10 in our simulation) with a low over-cost. The relative difference on the mode volume computed by the FEM and the WGMode toolbox is 6.6 10 -5 .

When considering the TM mode with the same indexes = m = 302, n = 1, we found a resonance wavelength λ = 1498.4979 nm and volume V = 2334.89 µm 3 with the WGMode toolbox. Computation of the mode volume took 4.4 sec. Compared to the computation of the TE mode volume, this higher computation time is explained by the computation of the maximum value of the electric field that is more complicated for TM modes, see Section 3.2. For TM modes the only comparison is with the FEM. The FEM provided a resonance wavelength λ = 1498.4923 nm and a mode volume V = 2335.59 µm 3 . The relative difference between the mode volume given by the two methods is 3.0 10 -4 .

As a second benchmark test, we consider a micro-sphere with radius 25 µm and optical index 1.453 and the TE mode with indexes = 271, m = 271. The mode volume obtained by formula (43) is 319.03 µm 3 . Using the WGMode toolbox we found a resonance wavelength of 808.6104 nm and a mode volume of 319.21 µm 3 . The values provided by the FEM are a resonance wavelength of 808.6073 nm and a mode volume of 319.30 µm 3 . The FEM computation with Lagrange P 3 Finite Elements and 25998 mesh elements took 131.7 sec. (10 modes are computed at the same time) whereas using the WGMode toolbox, the mode volume was computed in 1.6 sec. The relative difference between the mode volume given by the two methods is 2.8 10 -4 .

As a last benchmark test, we consider a micro-sphere with radius 25 µm and optical index 1.453 and the TM mode with indexes = 271, m = 270 and n = 2. Formula (43) can not be used for such modal index values. Using the WGMode toolbox we found a resonance wavelength of 780.9564 nm and for the mode volume 542.83 µm 3 . The computation took 4 sec. The values provided by the FEM are a resonance wavelength of 780.9533 nm and a mode volume of 539.88 µm 3 . The computation with Lagrange P 3 Finite Elements and 25998 mesh elements took 262.0 sec. (10 modes are computed at the same time). The relative difference between the mode volume given by the two methods is 5.3 10 -3 .

Conclusion

We have presented a Matlab toolbox aimed at studying whispering gallery modes (WGM) in spherical optical micro-resonators. Its originality comes from analytical formulas used to compute the volume of WGM which provide fast and accurate results. Such a software development has been made possible because of the existence of analytical expressions for the electromagnetic field of a WGM in a spherical dielectric cavity and the existence of exact quadrature formulas for various integrands involving Bessel's functions in the literature. It would be difficult to extend the toolbox to other geometries of the dielectric cavity. However, a future development of the WGMode toolbox could be envisaged to cover optical micro-bubble resonators [START_REF] Barucci | Optical microbubble resonators with high refractive index inner coating for bio-sensing applications: An analytical approach[END_REF].

Figure 1 :

 1 Figure 1: Definition of the spherical coordinate system used.

  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 semi-normalization; -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;

Figure 2 :

 2 Figure 2: General organization of the WGMode toolbox. R denotes the sphere radius, N its optical index. λ is the resonance wavelength of the mode and , m its indexes. The scripts ELLRES and WVLRES are used to compute the resonant mode features around a wavelength target value. The scripts QFACT, VOLMOD can be used to compute key attributes of the mode. PLTMOD allows to display the mode.

  >> wvlresMicro-sphere radius [micro-m] = 50 Optical index of the micro-sphere = 1.5 TE or TM mode ? (TE/TM) : TE Mode index ell = 302 Wavelength interval of interest found [nm] = [1041.6229 , 1564.7545] Radial mode number and resonance wavelength [nm]:

Figure 3 :

 3 Figure 3: Radial behavior of the TE mode with polar index = 302 and radial index n = 5 of a micro-sphere with radius R = 50 µm and optical index N = 1.5.

  -Local minimum possible. fminunc stopped because the size of the current step is less than the default value of the step size tolerance. <stopping criteria details> -

  >> qfact Micro-sphere radius [micro-m] = 50 Optical index of the micro-sphere = 1.5 TE or TM mode ? (TE/TM) : TE Mode index ell = 302 Approximated real value of the resonance wavelength [nm] = 1502 -

  -Non-linear solver fsolve informations Equation solved. fsolve completed because the vector of function values is near zero as measured by the selected value of the function tolerance, and the problem appears regular as measured by the gradient. <stopping criteria details>

  -Complex resonance wavelength for the mode in [nm] = 1502.024424398681 -9.518987898596565e-15i Quality factor for the mode = 1.578e+17

  radius [micro-m] = 50 Optical index of the micro-sphere = 1.5 TE or TM mode ? (TE/TM) : TE Mode index ell = 302 Approximated real value of the resonance wavelength [nm] = 1502 One root (wavelength resonance) of the modal equation located in the disk centered at 1502 nm with radius 1 nm Complex resonance wavelength for the mode in [nm] = 1502.024424398681 -2.593625978591534e-17i Quality factor for the mode = 5.791e+19

Acknowledgments

The author would like to thanks Patrice Féron, from Foton Laboratory, ENSSAT-University of Rennes 1, for the numerous discussions on optical microresonators and for his constant support.

-HANSENN Hansen solution N m of the spherical vectorial wave equation.

Lastly, the WGMode toolbox contains low level internal functions:

-MODEQEL Modal equation [START_REF] Srinivasan | Cavity Q, mode volume, and lasing threshold in small diameter AlGaAs microdisks with embedded quantum dots[END_REF] for TE or TM modes for a fixed wavelength and a variable mode index;

-MODEQWV Modal equation ( 14) for TE or TM modes for a fixed mode index and a variable wavelength;

-CMODEQWV Modal equation [START_REF] Muljarov | Exact mode volume and Purcell factor of open optical systems[END_REF] for TE or TM modes for a fixed mode index and a variable wavelength;

-PLOTD plots a discontinuous function removing the draw line at function discontinuity jumps;

-VSHNOR computes the square of the euclidean norm of the vector spherical harmonics;

-FOPTE/FOPTM objective function (and its gradient) required to compute the maximum of the square of the euclidean norm of the electric field for TE/TM modes;

-MODEFIELD computes the electromagnetic field (electric field and magnetic induction) of a given mode;

-INTGDL integrand in the Delves-Lyness Method used to solve the modal equation in the complex plane.

A comprehensive test run

In order to illustrate the use of the WGMode toolbox, we consider a microsphere of radius 50 µm and optical index 1.5, see [START_REF] Righini | Whispering gallery mode microresonators: Fundamentals and applications[END_REF]. We are interested in resonance conditions around a wavelength of 1500 nm. First, we use the script entitled ellres to determine the solutions to the approximated modal equation [START_REF] Srinivasan | Cavity Q, mode volume, and lasing threshold in small diameter AlGaAs microdisks with embedded quantum dots[END_REF] for wave numbers k and k 0 corresponding to a wavelength of 1500 nm and to a TE mode. index N = 1.5, let us consider the WGM characterized by the mode indexes = 302, m = and n = 1 for which the resonance wavelength has been found to be λ = 1502.024424 nm. Thanks to the volmod script in the WGMode toolbox, we find that the volume of this WGM is 2259.6416 µm 3 as detailed below.


Daniel Gandolfo

Christian Maes

Jean Ruiz

Senya Shlosman

Glassy states: the free Ising model on a tree

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

It is well known that the Ising model on a regular Cayley tree undergoes a second order phase transition at the critical temperature T cr , below which the Gibbs states µ + and µ -, corresponding to + and -boundary conditions, are different extremal states. Unlike the usual Z d -lattice case, on the tree the behavior of the free state µ ∅ , corresponding to empty boundary conditions, is very rich. On Z d we have µ ∅ = 1 2 (µ + + µ -), while on the tree that is trivially true only in the uniqueness regime. Moreover, the free state is extremal for temperatures T below the critical temperature T cr , until a certain spin-glass temperature T SG , below which it stops to be extremal. The question of finding the temperature range of the ergodicity of the state µ ∅ was open for twenty years, and was solved by Bleher, Ruiz and Zagrebnov in their 1995 paper [START_REF] Bleher | On the purity of the limiting Gibbs state for the Ising model on the Bethe lattice[END_REF]. Soon after, a simpler argument was provided by Ioffe, [START_REF] Ioffe | A note on the extremality of the disordered state for the Ising model on the Bethe lattice[END_REF]. For a closely related communication theory problem, see the "Broadcasting on trees" paper by Evans, Kenyon, Peres and Schulman [START_REF] Evans | Broadcasting on trees and the Ising model[END_REF].

In the present paper we study the free state at temperatures below T SG .

A principal static feature of the spin glass phase is the presence of infinitely many pure states; see the discussion in [START_REF] Newman | Short-range spin glasses: results and speculations[END_REF] and the references therein. By a gauge transformation the spin glass on the tree and the ferromagnet are equivalent, except that random boundary arXiv:1709.00543v3 [cond-mat.stat-mech] 6 May 2019 conditions for the ferromagnet correspond to fixed boundary conditions for the spin-glass, as was for example discussed in [START_REF] Chayes | A mean field spin glass with short rang interaction[END_REF]. We show that the same phenomenon also happens in the ferromagnetic Ising model on the tree, i.e., without randomness in the interaction. Namely, the free state of the Ising model below the spin-glass temperature has a decomposition into extremal states which involves uncountably many extremal states. That is why we call this state 'glassy'. Our result answers an older question of Arnout van Enter, [15].

The next section contains the more detailed question with some notations and definitions.

Our main decomposition result is given in Section 3. The remaining sections contain further details and proofs.

II. NOTATIONS AND DEFINITIONS

Let T k = (V, E) be the Cayley tree with branching ratio k ≥ 2. We consider the nearest neighbor Ising model, where spins σ x = ±1 have a Gibbs distribution at temperature T = 1/β with boundary condition η in a finite volume Λ given by

µ(σ) = Z -1 exp    β x,y J xy σ x σ y + β x,y J xy σ x η y    (1)
Both sums run over nearest neighbors pairs, the first being over the pairs x ∈ Λ, y ∈ Λ, and the second one runs over sites x ∈ Λ, y / ∈ Λ. The infinite-volume Gibbs distributions are obtained as the convex hull of the set of all possible limit point when Λ grows to cover all the vertices of the tree. The set of Gibbs distributions for a fixed temperature and interaction is convex, and its extreme points are called pure states: they cannot be decomposed into other states.

In the ferromagnetic case we put J xy = 1, while in the spin-glass model the interaction is random: J xy = ±1 with probability 1/2 independently for any pair x, y . We can also consider the random interaction

J xy =    -1 with probability p +1 with probability 1 -p (2)
which interpolates between the ferromagnetic model p = 0 and the spin-glass model p = 1/2.

It is well known that in the ferromagnetic case the phase transition happens at the critical temperature T cr = 1/ arctanh(1/k). For T > T cr there is a unique infinite-volume Gibbs distribution and below that critical temperature the spontaneous magnetization m * (T ) becomes nonzero, while plus and minus boundary conditions give rise to different states:

σ 0 T ± = ±m * (T ), m * (T ) > 0 iff T < T cr .
We use the convention that the boundary condition by which the infinite-volume Gibbs state is obtained is indicated by a subscript to the expectation • , and superscripts will be used to indicate further parameters like temperature or the choice of model in [START_REF] Bleher | On the purity of the limiting Gibbs state for the Ising model on the Bethe lattice[END_REF].

There is yet another special temperature, T SG = 1/ arctanh 1/ √ k , called the spin-glass temperature. It will often appear in what follows below, and it has various interpretations.

For the ferromagnetic model it is known that the free state, the infinite-volume Gibbs distribution obtained by putting η ≡ 0 in (1), is extreme for T > T SG while it is not for T < T SG ; see [START_REF] Bleher | On the purity of the limiting Gibbs state for the Ising model on the Bethe lattice[END_REF], [START_REF] Ioffe | A note on the extremality of the disordered state for the Ising model on the Bethe lattice[END_REF], [START_REF] Evans | Broadcasting on trees and the Ising model[END_REF], [START_REF] Martinelli | Glauber dynamics on trees: boundary conditions and mixing times[END_REF]. Hence the question that motivates the present paper: what is the decomposition of the free state • T ∅ (at these lower temperatures) into extremal ones? Part of the question is also to understand why that extremality of the free state exactly stops at the spin-glass temperature, which in its origin characterizes a transition to glassy behavior in the spin-glass model. For example consider the spin-glass state obtained via plus-boundary conditions and look at the (random) local magnetization at the origin:

σ 0 T,SG + = m(T, {J xy })
It depends on the independent random variables J xy , which take values ±1 with equal probability. The spin-glass temperature T SG is that temperature below which m(T, {J xy }) is a fluctuating quantity with a non-trivial distribution. For example, the second moment, the Edwards-Anderson parameter, is positive only there:

q EA (T ) = E[m 2 (T, {J xy })] > 0 iff T < T SG ,
where the expectation E[•] is taken over the randomness {J xy }.

III. DECOMPOSITION OF THE FREE STATE

In this section we present the decomposition of the free state into pure states and we explain that a continuum of them enters into it, at least at low temperatures. Presumably, it is the case for all temperatures in the spin-glass region.

A. The construction of the decomposition

For any Gibbs state µ corresponding to temperature T we have that

µ (•) = µ(dσ) [ • T σ ]
where σ is a spin configuration drawn from µ and used as boundary condition at infinity. The Gibbs distributions • T σ are µ-almost surely extreme. Hence, we have here a decomposition of µ into extreme Gibbs distributions, but obviously the states • T σ might well be the same for different σ-s. Nevertheless this decomposition is nontrivial when µ is not extremal. It remains then to see how many and which different extremal states we get.

For the Ising model free state µ T ∅ we will use the Edwards-Sokal representation, [START_REF] Georgii | The random geometry of equilibrium phases[END_REF]:

µ T ∅ (•) = E q(T ) µ ES (•|n) (3)
(In contrast to the more standard notation we prefer here to call q(T ) = 1 -tanh 1/T the probability of removed (or closed) bonds.) On the tree, the random cluster measure is generated by independent bond percolation and n is the resulting random bond configuration over which we take the expectation E q . The open bonds generate a partition of the tree into maximal connected components. The measure µ ES (dσ|n) is supported by the spin configurations σ which are constant on each connected component as specified by n; these constants take values ±1, independently with probability 1/2.

We can rewrite (3) on the tree by ordering the components using the following definitions.

Let D be a subset of bonds of the tree T k , and consider the two Ising spin configurations σ D,± on T k defined as:

σ D,+ 0 = +1, σ D, - 0 = -1 (4)
and

σ D,± x = -σ D,± y for (x, y) ∈ D σ D,± x = +σ D,± y for (x, y) / ∈ D
That is, we fix the value of the spin at the root (say 0) to be +1 or -1, and the nearest neighbor spins alternate iff the corresponding bond belongs to the set D.

By • T σ D,ω we denote the Gibbs state of the ferromagnetic Ising model at inverse temperature T with the boundary condition σ D,ω where ω = ±1 corresponds to the way the spin at the origin is chosen.

Let p ∈ (0, 1). Take the set D to be random: every bond decides to be in D with probability p independently of the other bonds. Denote by E p the expectation with respect to that process.

Proposition III.1. The following decomposition of the free state for the ferromagnetic Ising model on the tree holds for all temperatures T :

• T ∅ = 1 2 E p(T ) • T σ D,+ + • T σ D,-, (5)
where

p(T ) = 1 2 [1 -tanh 1/T ].
Proof. We apply the Edwards-Sokal representation (3) of the Ising model. Start by noting that

• T ∅ = E q(T ) µ ES (•|n) = E q(T ) • T σ µ ES (dσ|n),
Consider now for a given bond collection n the atomic measure µ ± (dσ|n) = 1 2 (δ σ n+ + δ σ n-), which to n assigns two configurations defined by the relations (4), each with probability 1 2 . In other words, fixing the spin of the origin and fixing n determines all the other spin values, where in particular neighboring connected components of open bonds alternate their spin. However the resulting spin configuration would have the same distribution as in the Edwards-Sokal representation with twice as large probability of closed bonds: with p(T ) = q(T )/2,

E q(T ) [µ ES (dσ|n)] = E p(T ) [µ ± (dσ|n)],
Bringing all that together we conclude that on the tree (3) reduces to the formula [START_REF] Carlson | Bethe lattice spin glass: the effects of a ferromagnetic bias and external fields. I. Bifurcation analysis[END_REF].

Note also that the states • T σ D,+ , • T σ D,-can be obtained as thermodynamic limits of the finite-volume Gibbs states with the boundary conditions σ D,+ , σ D,-. These limits exist for

• T ∅ -almost all boundary conditions σ D,+ , σ D,-.
To see that this decomposition ( 5) is non-trivial for T < T SG it suffices to show that when

T < T SG , • T σ D,+ = • T σ D,-for E p(T ) -typical sets D.
That follows from the relations

E p(T ) [ • T σ D,+ ] > 0 > E p(T ) [ • T σ D,-], (6)
which are a special case of the Theorem 1.1 of [START_REF] Evans | Broadcasting on trees and the Ising model[END_REF].

Of course, the states • T σ D,+ may coincide for different sets D -this is the case at high temperature. In the next subsection we will show that at low temperatures there is a continuum of different states • T σ D,+ as we vary over the sets D, see also [START_REF] Gandolfo | A Manifold of Pure Gibbs States of the Ising Model on a Cayley tree[END_REF]. There we have shown that if the set D consists of bonds sufficiently separated from one another, then the configuration σ D,+ is a stable ground state. Of course, for our random configuration D this is not the case; for example we will see in D pairs of bonds sharing a vertex, which will happen with positive density. However, quite often we will see just isolated bonds, well separated, once p(T ) is small enough, see the next subsection for more details.

Remark. At all temperatures T the free-state two-point function σ 0 σ x T ∅ → 0 goes to zero when x goes to infinity, yet for T < T SG the free state is not extreme. In particular, the magnetization in increasingly large volumes has a variance that does not go to zero with the size of the volume. The interesting tail-observable which shows that the free state is not extreme is related to the Edwards-Anderson parameter. Here is the simplest version: take the magnetization at the origin,

M (τ ) = σ 0 T τ
in the infinite-volume Gibbs distribution with boundary condition τ ; that τ is drawn from the free state at temperature T . For T < T SG the random variable M (τ ) has a non-trivial distribution. , for all N large enough. We call the bond b to be r-strongly-frustrated (or just r-frustrated) in the state

• T =0,B N (b) σ D,+
, if the event σ x σ y = -1 happens with probability one in the state

• T =0,B N (b) σ D,+
, as well as the events σ x σ y = 1 for all bonds b = (x , y ) within distance r from the bond b, again for all N large enough.

For example, the above will hold if b ∈ D, while D is a deterministic configuration composed from isolated bonds which are sufficiently far away from each other, see [START_REF] Gandolfo | A Manifold of Pure Gibbs States of the Ising Model on a Cayley tree[END_REF], [START_REF] Gandolfo | A Manifold of Pure Gibbs States of the Ising Model on the Lobachevsky plane[END_REF].

What we want to show now is that if D is random, and b ∈ D, then it is very likely that b is r-frustrated, provided p is small enough (depending on r). Once we show that, our claim about mutual singularity will be proven, since for two independent configurations D , D we will be able to find arbitrarily large disjoint sets of r-frustrated bonds.

So let D be random, and b ∈ D. Our first observation is that the probability of D having other bonds at distance 2r from b is quite small, provided p is small enough. That would be the end of the story if the configuration σ D,+ would be a ground state configuration. Indeed, in that case the state

• T,B N (b) σ D,+
would be a small perturbation of the configuration σ D,+ once T is low, uniformly in N.

However, the configuration σ D,+ is not a ground state configuration a.s., so the state

• T =0,B 2r (b) σ D,+
might have other frustrated bonds in B 2r (b) ; moreover, it even can happen that b itself is not frustrated in this state. We will show now that all this is highly unlikely, once p is small enough. So suppose the set of frustrated bonds of the state

• T =0,B 2r (b) σ D,+ is not the singleton {b} . That can happen iff there is a contour γ, [γ ∪ Int (γ)] ∩ B 2r (b) = ∅, crossing certain number ≥ k + 1 of bonds of T k , such that |γ ∩ D| ≥ 2 .
Consider the set T γ of the bonds of T k which are inside γ, and the set L γ of bonds of T k the contour γ is intersecting. Together they form a finite tree S γ , which has the same branching number k as our infinite tree T k .

The set L γ is the set of all leaves of the tree S γ . Let n γ be the number of nodes of S γ inside γ, and Lγ ⊂ L γ be the intersection L γ ∩ D. So we have

|T γ | = n γ -1, |L γ | = , Lγ ≥ 2 ,
and

|L γ | = 1 + n γ (k -1) .
Evidently, the probability of seeing such a tree

Pr S γ , L γ , Lγ ≤ p /2 , so Pr (S γ , L γ ) ≤ 2 p /2 .
The number of trees S with n inner nodes does not exceed k 2n . Thus the probability that a given bond b 1 is a leaf of such a tree with leaves is bounded from above by 2 k 2 /(k-1) p /2 , which is exponentially small in for p small enough. So we can apply the standard Peierls argument to shows that the probability of the event

{there is a contour γ, such that γ ∩ B 2r (b) = ∅ or B 2r (b) ⊂ Int (γ)}
goes to zero as p → 0, which ends the proof.

To conclude, we point out for clarity that the "Gibbs ground" states lim T →0 • T σ D,+ constructed from the p-random bond configurations D, are typically nontrivial measures, i.e. they have infinite supports, a.s. This is in contrast with the ground states constructed in [START_REF] Gandolfo | A Manifold of Pure Gibbs States of the Ising Model on a Cayley tree[END_REF], where the corresponding Gibbs ground states are supported by a single ground state configuration. However, as is explained above, a vast majority of the frustrated bonds under typical ground state • T =0 σ D,+ are isolated bonds, once p is small. This fact is the source for the decomposition [START_REF] Carlson | Bethe lattice spin glass: the effects of a ferromagnetic bias and external fields. I. Bifurcation analysis[END_REF] to have a continuum of extremal components.

IV. DOUBLE-TEMPERATURE ISING MODEL

We already mentioned in the introduction that the phase diagram of the ferromagnetic Ising model is essentially determined by the critical temperature T cr = 1/ arctanh(1/k), and the spin-glass temperature T SG = 1/ arctanh 1/ √ k . A clarification of the situation can however be obtained by enlarging the objects in ( 5) into a two-temperature setting.

We consider two-temperature states, with T 2 the bulk temperature and T 1 the boundary temperature,

ν(T 1 , T 2 ) := • T 2 σ D(T 1 ),+
which is the infinite-volume Gibbs distribution at temperature T 2 with the boundary condition taken to be the spin configuration [START_REF] Carlson | Critical behavior of the Bethe lattice spin glass[END_REF] where D is drawn from E p(T 1 ) , the Bernoulli bond percolation process with parameter 1 -p(T 1 ). Of course, one may wonder whether the thermodynamic limits • T 2 σ D(T 1 ),+ exist. We are not going to prove it; what is said below holds for any limit point of that family. Note that (5) contains these states ν(T, T ) with

T 1 = T 2 =
T -and that is why it is useful to speak about the temperature T 1 , but of course the relevant parameter is the density p(T 1 ). The following is therefore presented in the (p, T )-plane, which is also the setting of [START_REF] Carlson | Critical behavior of the Bethe lattice spin glass[END_REF].

Consider the curve

T SG (p) = max{ 1 arctanh[ 1 k(1-2p) ]
, 0}.

Note that T SG (p) > 0 when k(1 -2p) > 1. IV.2. See also [START_REF] Carlson | Bethe lattice spin glass: the effects of a ferromagnetic bias and external fields. I. Bifurcation analysis[END_REF][START_REF] Carlson | Bethe lattice spin glass: the effects of a ferromagnetic bias and external fields. II. Magnetized spin-glass phase and the de Almeida-Thouless line[END_REF] for a more qualitative discussion.

E p ( σ x T σ D,± ) = 0 2. When T < T SG (p), E p ( σ x T σ D,+ ) = -E p ( σ x T σ D,-) > 0 Let us now look at the second moment, E p ([ σ x T σ D,+ ] 2 ) = E p ([ σ x T σ D,-] 2 ),
From Fig. 2 it is clear that the model shows an interesting and non-trivial behavior even on the line

T 1 = ∞.
That case is treated in [START_REF] Bleher | The Bethe Lattice spin glass at zero temperature[END_REF].

To understand the nature of the spin-glass temperature, we remark that the composition T SG (p(T )), which we abbreviate as

T SG (T ) := T SG (p(T )) = 1 arctanh 1 k tanh(1/T )
is an involution: T SG (T SG (T )) = T . In particular, T SG (0) = T cr , SG (T cr ) = 0, T SG (T SG ) = T SG .

Let us prove the last two Propositions in the vicinity of the point T 1 = T 2 = 0. After the analysis of the section III.B and using its technique, it is almost immediate.

Let us fix x to be the root 0 of our tree. Informally speaking, the magnetization at the root 0 in the state * T 2 σ D,+ is defined by the few bonds of the (rare) bond configuration D, which are in some proximity to 0. Moreover, this magnetization will take different values when these few bonds happen to be different. Since that happens with positive probability, our claim follows.

To be more formal, let B R be the ball of radius R centered at 0. Let b be a bond in B R . The proofs in general case involve the recursion relations, and can be obtained from the results of the papers [START_REF] Carlson | Bethe lattice spin glass: the effects of a ferromagnetic bias and external fields. I. Bifurcation analysis[END_REF], [START_REF] Carlson | Bethe lattice spin glass: the effects of a ferromagnetic bias and external fields. II. Magnetized spin-glass phase and the de Almeida-Thouless line[END_REF]. These results are summarized graphically by the following pictures.

FIG. 3: The first moment.

B

  . At low temperatures T the states • T σ D,+ are mutually singular. Theorem III.2. Pick two independent configurations σ D 1 ,+ , σ D 2 ,+ from the distribution • T ∅ . Then the two limiting states • T σ D 1 ,+ , • T σ D 2 ,+ exist and are mutually singular a.s. with respect to • T ∅ × • T ∅ , provided the temperature T is low enough. Proof. We denote by D a random configuration of bonds in T k , each bond picked independently with probability p, with the parameter p being fixed and small enough. We will study the Gibbs states • T σ D,+ at low temperatures, with the goal of showing their a.s. mutual singularity. Let b = (x, y) ∈ T k be a bond, and B N (b) ⊂ T k be a ball of radius N centered at b. Consider the ground state (=zero temperature Gibbs) measure • T =0,B N (b) σ D,+ in the box B N (b) with boundary condition σ D,+ . We call the bond b frustrated in the state • T =0,B N (b) σ D,+ , if the event σ x σ y = -1 happens with probability one in the state • T =0,B N (b) σ D,+

FIG. 1 :

 1 FIG.1:The curve T SG (p).

2 SGFIG. 2 :

 22 FIG. 2: Phase diagram. Phases I, II, III correspond to the behavior in Proposition IV.1,

0 T 2 σ 2 σ

 022 Define the set D b as the family of all realisations D which has b as the only bond in B R . Clearly, the probability Pr (D b ) is positive for every value of the parameter p (T 1 ) . Now, let b , b be two such bonds, with dist (b , 0) > dist (b , 0) , and let D ∈ D b , D ∈ D b be two typical configurations. Then, using the technique of the section III.B and a little of cluster expansions, one sees that there exists a constant c = c (b , b ) , such that σ D ,+ > c > σ 0 T D ,+ , provided both T 1 and T 2 are small enough. That proves the positivity of the variance of the EA random variable σ 0 T σ D,+ (with randomness coming from D).

Acknowledgments Part of this work has been carried out in the framework of the Labex Archimede (ANR-11-LABX-0033) and of the A*MIDEX project (ANR-11-IDEX-0001-02), funded by the "Investissements d'Avenir" French Government programme managed by the French National Research Agency (ANR). Part of this work, concerning the EA parameter, has been carried out at IITP RAS. The support of Russian Foundation for Sciences (project No. 14-50-00150) is gratefully acknowledged. This work was partially supported by the CNRS PICS grant "Interfaces aléatoires discrètes et dynamiques de Glauber" and by the grant PRC No. 1556 CNRS-RFBR 2017-2019 'Multi-dimensional semi-classical problems of Condensed Matter Physics and Quantum Mechanics ". CM thanks the hospitality of the CPT-Luminy at Marseille.

FIG. 4: The second moment.