E Garcia-Caurel
email: enric.garcia-caurel@polytechnique.edu

A Lizana

G Ndong

B Al-Bugami

C Bernon

E Al-Qahtani

F Rengnez

A De Martino

A mid-infrared Mueller ellipsometer with pseudoachromatic optical elements

The purpose of this article is to present a new broadband Mueller ellipsometer designed to work in the midinfrared range, from 3 to 14 microns. The Mueller ellipsometer, which can be mounted in reflection or in transmission configuration, consists of a polarization state generator (PSG), a sample holder, and a polarization state analyzer (PSA). The PSG consists in one linear polarizer and a retarder sequentially rotated to generate a set of four optimal polarization states. The retarder consists in a bi-prism made of two identical Fresnel rhombs disposed symmetrically and joined by optical contact, giving the ensemble a "V" shape. Retardation is induced by the four total internal reflections that the beam undergoes when it propagates through the bi-prism. Total internal reflection allows to generate a quasi-achromatic retardation. The PSA is identical to the PSG, but with its optical elements mounted in reverse order. After a measurement run, the instrument yields a set of sixteen independent values, which is the minimum amount of data required to calculate the Mueller matrix of the sample. The design of the Mueller ellipsometer is based on the optimization of an objective criterion that allows minimizing the propagation of errors from raw data to the Mueller matrix of the sample. The pseudo-achromatic optical elements ensure a homogeneous quality of the measurements for all wavelengths. The performance of the Mueller ellipsometer in terms of precision, and accuracy, is discussed and illustrated with a few examples.

1.Introduction

Standard spectroscopic ellipsometry is an optical characterization technique sensitive to the polarization of light, which is also noninvasive and non-destructive. Standard ellipsometers have been used for more than one century to characterize the optical properties of materials and also to measure thickness of coatings made of single films or stacks. Former ellipsometers were limited to work in the visible spectral range, but decades of work in instrumentation together with the progressive availability of new optical components allowed to push the spectral limits to the vacuum ultraviolet and to the infrared [START_REF] Irene | A brief History and State of the Art of Ellipsometry[END_REF]. There are two main reasons to explore the infrared range. The first one is the infrared absorption due to free carriers, it extends over a wide spectral range and depends on free carrier volume density and mobility. This feature can be used to assess the electrical properties of thin films made either of metals, semiconductors, or transparent conductors, which are commonly used to build low-emissivity windows, transparent electrodes for optoelectronic devices such as solar cells, light emitting diodes or liquid crystal displays. The second characteristic of the infrared is the occurrence of welllocalized and spectrally narrow absorption bands due to the resonant excitation of vibrational modes of single chemical bonds in molecules or collective oscillations in networks of atoms. Those resonances are highly specific of the type of chemical bonds and allow to study the chemical composition of samples. The room for ellipsometry appears because the absorption due to chemical bonds is highly directional. Directional absorption creates linear dichroism that is easily measured in polarimetry. Linear dichroism can then be used to evaluate structural properties of samples such as the degree of molecular ordering or the degree of crystallinity. The ordering of molecules can also induce birefringence which is easily measured by polarimetry and can be used as an indicator to characterize sample properties. Nowadays, there are several relevant industrial and scientific activities involving the use of polymer thin films such as: pharmaceuticals, alimentary packing, organic light emitting diodes synthesis, or flexible panel display fabrication that can profit from polarization sensitive optical measurements in the infrared [START_REF] Röseler | Spectroscopic Infrared Ellipsometry[END_REF][START_REF] Korte | Infrared spectroscopic ellipsometry: A tool for characterizing nanometer layers[END_REF][START_REF] Schubert | Infrared Ellipsometry on Semiconductor Layer Structures[END_REF].

Despite of its great interest, standard ellipsometry, is limited to non-depolarizing samples, i.e. samples which do not introduce any randomness to the temporal or spatial vibration of the electromagnetic field associated to a radiation beam. Mueller ellipsometry, does not suffer from this limitation and can be used to study non-depolarizing as well as partially depolarizing samples. Mueller and standard ellipsometry provide equivalent results for non-depolarizing samples. Inhomogeneity in the samples is an important and rather common cause of depolarization. Inhomogeneity can appear either in sample chemical composition or structure (thickness, roughness or porosity). Because samples currently found either in everyday life or manufactured by the industry, may not be perfectly homogeneous, they may induce depolarization to a given extent. Therefore, the use of Mueller ellipsometry appears as an interesting alternative to extend the characterization capabilities of ellipsometry. In the past ten years a number of spectroscopic Mueller ellipsometers have been developed. However, the vast majority of them were designed to work either in the visible [START_REF] Goldstein | Mueller matrix dual-rotating retarder polarimeter[END_REF][START_REF] Delplancke | Automated high-speed Mueller matrix scatterometer[END_REF][START_REF] Compain | Complete high-frequency measurement of Mueller matrices based on a new coupledphase modulator[END_REF][START_REF] Garcia-Caurel | Spectroscopic Mueller polarimeter based on liquid crystal devices[END_REF] or in the visible-near infrared spectral range [START_REF] Ladstein | Fast near-infra-red spectroscopic Mueller matrix ellipsometer based on ferroelectric liquid crystal retarders[END_REF]. A few of them can make measurements in the infrared spectrum [START_REF] Jones | A Simple Infrared Polarimeter[END_REF][START_REF] Chenault | Measurements of linear diatenuation and linear retardation spectra with a rotating sample spectropolarimeter[END_REF][START_REF] Vap | Optimization of a dual-rotating-retarder polarimeter as applied to a tunable infrared Mueller-matrix scatterometer[END_REF][START_REF] Stanislavchuk | Synchrotron radiation-based far-infrared spectroscopic ellipsometer with full Mueller-matrix capability[END_REF][START_REF] Deibler | Measurement of the complex refractive index of isotropic materials with Mueller matrix polarimetry[END_REF][START_REF] Thomson | Infrared ellipsometer/polarimeter system, method of calibration and use thereof[END_REF]. Our motivation was to develop a wide-spectral Mueller ellipsometer optimized to work in the mid-infrared to combine the sensitivity of ellipsometry to the physical structure of samples, with the chemical sensitivity characteristic of the infrared range. Such a Mueller ellipsometer will probably satisfy the growing demand of academic and industrial scientists and technicians, who need new information to characterize their samples. The instrument was patented a few years ago [START_REF] Garcia-Caurel | Bradband ellipsometer / polarimeter system[END_REF] and it was used in several studies [START_REF] Lizana | Enhanced sensitivity to dielectric function and thickness of absorbing thin films by combining total internal reflection ellipsometry with standard ellipsometry and reflectometry[END_REF][START_REF] Hu | Maximum likelihood method for calibration of Mueller polarimeters in reflection configuration[END_REF]. The aim of this article is to present a detailed description of the instrument design and operation in a format different from that of a patent, understandable to those who are not familiar with legal documents. The characteristics that make this apparatus different to others, is the simple operation protocol, based in the non-continuous rotation of two retarders, and the use of achromatic optical elements to guaranty an optimal and uniform quality of measurement over the measured spectral range. In the first part of this paper we describe the design principle and the practical implementation of the Mueller ellipsometer. In the second part we show the results of some tests to assess the performances of the apparatus in terms of precision and accuracy. Finally, in the third part we provide two examples to illustrate representative results that can be obtained with the Mueller ellipsometer.

A. Optical configuration and instrument operation 1. Optical configuration

The basic configuration of a Mueller ellipsometer consists of a light source, an input arm, a sample holder, an exit arm, and an acquisition system. In order to preserve the achromaticity of the system, gold coated mirrors have been used to define the optical path. The input arm includes an infrared source, a polarization state generator (PSG) and a retractable sample-holder for calibration samples. The exit arm includes also a mobile sampleholder for calibration samples, a polarization state analyser (PSA) and a detector system. The illumination source is a commercial FTIR spectrometer providing an infrared beam in a continuous spectral range from 2.5 to 14 µm. The PSG comprises a fixed grid type linear polarizer, and an achromatic retarder mounted in a rotatable holder. The rotatable holder allows orienting the retarder axis at four consecutive positions respect to the transmission direction of the fixed polarizer. Each one of the four orientations gives rise to a particular polarization state. The PSA is identical to the PSG, but with the optical elements, (retarder and polarizer) positioned in reverse order. Each one of the four polarization states created by the PSG, after being transformed by the sample, is projected to four optical configurations sequentially created by the PSA. Each analysis configuration is made by orienting the PSA retarder respect to the transmission axis of the polarizer at a given angle. Similarly to the PSG, the polarizer in the PSA is kept still during all the measurement process to avoid artifacts created by polarization sensitive detectors. Although measurements can in principle be performed either in reflection or in transmission, we have chosen to work in reflection because it is the most sensitive to measure thin films on isotropic substrates. In Fig. 1 we have represented two possible variants of the Mueller ellipsometer. In one configuration the beam incident on the sample is collimated, whereas in the other configuration the beam is focused. The collimated beam configuration is adequate to measure gratings or samples showing diffraction or scattering. This configuration provides a well-defined angle of incidence with controllable divergence (usually it is set to be around 1mrad) thanks to the aperture diaphragm (D2) placed at 100 mm from the second mirror (M2). The spot size on the sample, which is half of that of the aperture diaphragm because of the conjugation of focal lengths of mirrors M3 and M4, is 2x(2/cos(AoI)) mm. Here, AoI stands for angle of incidence. The spot size is not circular because of the oblique incidence. The focused beam configuration is well adapted to measure specular samples. It has the advantage of providing a smaller spot size than the collimated configuration, and the possibility to collect more light because the aperture diaphragm can be reasonably opened. The disadvantage of the focussed configuration is that the angle of incidence is less well-defined than in the collimated configuration. Indeed, focalization implies the existence of a distribution of angles of incidence across the focused beam waist, estimated to vary from 0.5 to 0.1 degrees depending on the diameter of the aperture diaphragm used. These values are related to the size of the beam and the focusing mirrors used in our set-up. Uncertainty related to the angle of incidence can be taken into account when modelling experimental data [START_REF] Irene | A brief History and State of the Art of Ellipsometry[END_REF].

The angle of incidence is not fixed and can be varied. Our first attempt was to use the classical 'goniometric' method which consists of rotating the sample and the exit arm by an angle α and 2α respectively. However, we disregarded this method because it was very difficult to keep the instrument properly aligned. Our second attempt was inspired but not identical to the one discussed by Matsumoto et al. [START_REF] Matsumoto | Measurement of the dielectric constant of thin films by terahertz time-domain spectroscopic ellipsometry[END_REF]. To change the angle of incidence by ω deg., the parabolic mirrors placed just before and after the sample must be pivoted by an angle ω and -ω respectively. After pivoting the mirrors, the sample must be shifted by a given distance Δz to compensate for the beam displacement. In practice, Δz is the distance required to keep the beam impinging on the same position of the sample surface, and to maximize the signal level measured by the detector. The method is sketched in Fig. 2. The main difference of our method and the one discussed in [START_REF] Matsumoto | Measurement of the dielectric constant of thin films by terahertz time-domain spectroscopic ellipsometry[END_REF] is that we worked with a collimated beam instead of a focused beam. In our case, the advantage of using a collimated beam is that the variation of the distance, Δz, does not affect either the alignment nor the divergence of the beam across the exit arm of the Mueller ellipsometer.

Instrument operation and optimization

As previously stated, the PSG generates a sequence of four independent polarization states. After reflection on the sample, the polarization states are projected on the PSA, which analyses with a sequence of four independent polarization configurations each polarization state arriving to it. As a result, a complete measurement yields a set of 16 independent intensity values that allow the calculation of the Mueller matrix. The four states of polarization generated by the PSG, as well of the optical configurations of the PSA, are not obvious or intuitive and they have been selected according to an objective criterion, called the maximization of the condition number [START_REF] Taylor | Theory and Applications of Numerical Analysis[END_REF], which allows the minimization of the propagation of errors from the measured raw-data to the computed Mueller matrix.

In the following we review the physical model used to design the optical setup of the PSG and the PSA and to define the condition number. It is common to use two 4x4 matrices, W, and A, to represent the PSG and the PSA respectively. The columns of the W matrix correspond to the four Stokes vectors generated by the PSG. Reciprocally, the analysis configurations used by the PSA are represented by four transposed Stokes vectors written as rows in the A matrix. The result of a polarimetric measurement consists of sixteen intensity measurements that can be expressed in terms of a matrix:

, B AMW = ( 1
)
where M is the Mueller matrix of the sample.

If A and W are known, then M can be extracted from the raw data B as:

1 1 . M A BW - - = (2)
The determination of A and W is called calibration of the Mueller ellipsometer. Obviously, the matrices A and W must be non-singular. Moreover, in order to minimize error propagation from B to M, the analysis and modulation matrices, A and W respectively, must be "as close as possible" to unitary matrices. The best criterion to preserve A and W from singularity, is to optimize their condition numbers c(A) and c(W), which we define as the ratio of the smallest over the largest singular value of matrices A and W respectively.

In order to study the dependence of the condition number with the representative parameters of a given optical setup, a good approach is to create a functional model. Our PSG consists of a fixed linear polarizer followed by a retarder that can be rotated and placed at four different orientations respect to the polarizer. To simulate the optical behaviour of the PSG we treat the polarizer and the retarder as ideal components. The ideal Mueller matrix of the polarizer, P, and that of a retarder with a phase-shift "δ", and azimuth "θ", C(θ, δ), can be found in many publications and books such as "Handbook of Optics" [START_REF] Chipman | Polarimetry[END_REF] or [START_REF] Goldstein | Polarized Light[END_REF] which focusses on polarimetric instrumentation. Assuming that the light beam entering the PSG is completely depolarized, the Stokes vector of the vector leaving the PSG is given by the expression:

( ) ( ) ( ) , , 1 0 0 0 . T out S C P θ δ θ δ = ⋅ ⋅ (3)
The W matrix can then be built by taking four Stokes vectors characterized, each one, by a given azimuth θi for i=1 to 4, and a common value for the retardation "δ". Under the above representation, the condition number of the matrix W can be understood as a function of five parameters, the four orientations and the retardance. As shown by J. S. Tyo [START_REF] Tyo | Design of optimal polarimeters: maximization of signal-to-noise ratio and minimization of systematic error[END_REF], the theoretical value of the condition number for a PSG or a PSA consisting in one polarizer followed by at least one retarder, can vary between "0" and "3 -1/2 ". Accordingly, the values of the five parameters (retardation and azimuths), used to design the PSG, must be those that optimize the value of the condition number.

Finding the value of the parameters that maximize the condition number is a numerical problem that can be solved by means of standard numerical optimization packages. The resulting values for the retardation are: either 132°+n180° or 228°+n180°. The best sets of azimuths are: {θ1, θ2, θ3, θ4} = {38°+n180°, 74°+n180°, 106°+n180°, 142°+n180°}, "n" being a positive integer.

The tolerance of the absolute value to the retardation, as well as the four orientation angles of the retarder can be estimated with the help of the condition number. It is conventionally accepted that a minimum value of 0.25 is needed to obtain measurements with acceptable quality. Accordingly, the tolerances of the phase-shift, and the four azimuths of the retarder around their respective optimum value, are specified in order to keep the condition number higher than 0.25. The resulting tolerance range for the retardation and the orientation angles is ±15° and ±20° for respectively.

In view of the constraints imposed to the value of the retardation by the optimization of the condition number, and the wide spectral range intended to be measured, an achromatic retarder is needed. To achieve this objective, a total internal reflection (TIR) based achromatic retarder has been designed and built. The retarder consisted of a bi-prism made of two identical Fresnel rhombs disposed symmetrically and cemented together giving to the ensemble a "V" shape. An interesting review of alternative prism shapes has been written by J. M. Bennet [START_REF] Bennett | A Critical Evaluation of Rhomb-Type Quarterwave Retarders[END_REF]. The "V" shape bi-prism induces four total internal reflections and prevents the beam emerging from the bi-prism to be deviated from its original direction. The total retardance created by such "V" shaped bi-prism in contact with the air is given by: ( )

1 2 2 2 1 2 cos sin 1 8 tan , sin bi prism n n φ φ δ φ - -   -   =       ( 4
)
where φ is the internal reflection angle of incidence, and "n" is the refraction index of the bi-prism. The prisms have been built in ZnSe because this material is transparent over a wide spectral range in the infrared. Because of the refractive index of ZnSe in the infrared, 2.429, the prisms must be cut to an angle of 58.25° to get the optimal retardation of 228°.

The strategy chosen to calibrate the broadband Mueller ellipsometer is based on the eigenvalue method (ECM) described by E. Compain et al. [START_REF] Compain | General and selfconsistent method for the calibration of polarization modulators, polarimeters, and Mueller-matrix ellipsometers[END_REF] because of its robustness. The ECM requires the measurement of a set of reference samples, called by some authors calibration samples. According to Compain's method we used the following calibration samples. A linear polarizer oriented respectively at 0° and 90° to the plane of incidence, and an achromatic retarder providing a total retardation of 110° oriented at 30° to the plane of incidence. We used as calibration retarder another "V" shaped prism. The calibration samples were placed between the sample holder and mirror M4 (see Fig. 1(a)) to calibrate the PSG, and between the sample holder and mirror M5 to calibrate the PSA (see Fig. 1(a)).

B. Experimental validation of the optimized optical configuration

In this section, we show some results to discuss the conformity of the real Mueller ellipsometer to the theoretical model used to optimize the optical configuration. Moreover, based on the variance-covariance matrix formalism we provide a detailed evaluation of the noise (random error) propagation to each element of experimental Mueller matrices. This approach is useful to evaluate the uncertainty affecting each Mueller matrix element, in other words, the quality of measurements provided by the system. Fig. 3 shows the condition number of two typical experimental W and A matrices determined after calibration. As expected, the condition number of both, W and A, is quasiachromatic. The slight spectral dependence observed in Fig. 3, can be attributed, according to Eq.( 4), to the refractive index dispersion of the retarders. The spectral dependence is slightly higher for c(W) than c(A), but in practice those differences have a negligible impact on the quality of data. In our case, c(W) values vary between 0.40-0.45 and C(A) varies between 0.44 and 0.46, which are lower than the optimal cmax = 3-1/2 = 0.58 predicted by the ideal model and targeted in the design. This fact is explained because the real shape of the prisms does not correspond exactly to the optimal. Indeed, the prisms were cut to an angle of 60° instead of an ideal angle of 58.25° (Eq.( 4)). Accordingly, the real retardance provided by both retarders was measured to be 215°. However, as it will be discussed in the following, even though the real condition number is lower than the ideal one, the Mueller ellipsometer still provides measurements of high quality. The condition number formalism gives an overall estimation of noise propagation from the intensity matrix, B, to the Mueller matrix, M. It fails, however, to provide a particular estimation of noise propagation for each Mueller matrix element. This drawback can be solved using an alternative approach based on the variance-covariance matrix shown in [START_REF] Anna | Optimal Mueller matrix estimation in the presence of Poisson shot noise[END_REF]. In the following, it will be assumed that the data of each matrix element in B, is characterized by a Gaussian noise with null mean and constant standard deviation σ0. It is also assumed that the variance is independent of the absolute intensity measured by the detector, thus it is possible to assume the same value, σ0, for the all the elements of matrix B. Accordingly the variance of each Mueller matrix element can be expressed by: [ ]

1 1 2 2 , 1 ,16 ,
.

T T i o i i i AA WW σ σ - -         ∀ ∈ = ⊗             (5)
In Eq.( 5) the matrices A and W, have been rewritten as 16 component vectors, according to the dactylographic convention. Noise propagation can be evaluated as the square root of Eq. 5. In Eq. 5 the factor following σ0 determines by how much the variance characterizing the noise in matrix B will be amplified by the effect of matrices A and W. Therefore, we called the square root of this factor, 'noise propagation factor'. In the following example, we show an evaluation of the noise propagation factor from intensity matrices to Mueller matrices in three different cases represented by a Mueller ellipsometer working with: i) W and A matrices with an ideal optimal condition number equal to 0.58. ii) our experimental A and W matrices, and iii) W and A with a poor condition number equal to 0.25. In all three cases σ0 is assumed to be equal to 1 for simplicity. Results are shown in Fig. 4. To evaluate noise propagation in practical situations it is necessary to use an estimation of σ0, obtained from real data, instead of the value 1 used for this particular example.

As expected, the experimental noise propagation factors (case ii) are very close to the optimal ones (case i). The small differences among them are due to the fact that the condition number of the experimental A and W is not strictly optimal for the reasons discussed previously. Fig. 4 clearly shows that the accuracy obtained with our Mueller ellipsometer is close to the optimal one and definitely better than that obtained with a poorly conditioned system (case iii). Data in Fig. 4 also shows that even for the optimal configuration, the variance is not equally distributed for all Mueller matrix elements. Accordingly, noise propagation affecting the elements of the first column and row, is lower than that for the remaining matrix elements. This information may be of interest for instance in practical applications for which experimental data have to be fitted to a model. In the latter cases, being able to evaluate the uncertainty for each Mueller matrix element helps to decide whether the model accurately fits the experimental data or not. Fig. 4. Noise propagation factors for each Mueller matrix coefficient for the three situations described in the text: dashed line (optimum condition number), grey line (experimental matrices) and black line (poorly conditioned system). Scales are kept identical to highlight the differences among the three cases and also to make clear that noise propagation is not identical for all the elements.

As previously discussed, to quantify the accuracy of the Mueller ellipsometer, it is necessary to evaluate σ0 from real data. This can be achieved easily by measuring a set of N independent intensity B matrices and then evaluating, for each matrix element, the corresponding standard deviation. Since we worked with spectroscopic data we had to repeat the evaluation for each spectral point. To study the dependence of σ0 with N, we calculated σ0 for different sets of independent measurements with N=1, 2, 4, 16, 32, 64 and 128 respectively. Data shown in Fig. 5 corresponds to the standard deviation at one particular and representative wavenumber, 1500 cm -1 . We repeated the same experience with two different spectral resolutions, 8 cm -1 and 16 cm -1 respectively. In both cases, a 1 N dependence is observed, being characteristic of additive Gaussian noise. The inset in Fig. 5 shows a typical intensity spectrum. As it can be seen, the beam intensity is not constant across the measured spectral range. It is quite feeble at the extrema and shows a maximum in the central part, around 2500 cm -1 . The peak (at 2382.66 cm -1 ) and the two bands (at 1545.06 cm -1 and 3781.03 cm -1 ) are due, respectively, to the carbon dioxide and water vapor present in the atmosphere. If the instrument was operated in a vacuum chamber, those absorption bands would disappear. Because the signal level is not constant, the noise level, σ0, is higher at the ends than at the center of the measured spectral range. The error bars in Fig. 5 account for the uncertainty in the determination of σ0. The dependence of the standard deviation with N, shows a clear advantage of averaging 16 to 32 independent measurements because the noise level can be considerably reduced to an acceptable value without requiring and excessive acquisition time. Increasing the number of accumulated measurements beyond 64, does not allow to improve signal quality significantly and makes measurements considerably long. In practice, the information from Fig. 5 is used to set the number N required to obtain a target value of σ0. In our case we used to work with a target value lower than 1•10 -3 . The time required to measure a full intensity matrix can vary from 4 to 8 minutes depending on the selected spectral resolution and the reflectivity of the sample. Once the value of σ0 is known, it is just necessary to apply Eq.( 5) to evaluate the uncertainty of each Mueller matrix element due to random noise. According to the results of Fig. 4, the minimum experimental uncertainty due to random noise propagation is estimated to be around 2•10 -3 (for Mueller matrix elements in the first and second rows), and about 5•10 -3 for the remaining matrix elements. The effect of systematic errors is not taken into account in our evaluation of σ0 because it depends on a number of empirical parameters that may change from experience to experience or from day to day. There are many sources of systematic errors such as room temperature fluctuations that can affect the retardance created by the biprisms, the composition of the atmosphere (% of water vapour or CO2 for instance) that may change its absorption, or drifts in the radiance of the infrared source among others. We have found that the main factor that influences the systematic errors in our experiences is the misalignment occurring when samples are mounted and removed from the sample holder during routine operation. To circumvent this problem a (He:Ne) laser diode is inserted in the optical path prior to each measurement to check that the sample is well positioned and that the beam follows the same path that was defined in the calibration step. As a result of a careful operation of the Mueller ellipsometer in a temperature stabilized room and in a controlled atmosphere, allow to bring the impact of the systematic errors to a level comparable to that of random noise previously discussed. Thus, we can estimate that the maximum average accuracy that can be reached with our experimental setup is about 0.005 for Mueller matrix elements in the first row and column; and around 0.01 for the remaining matrix elements. In the following we provide several examples.

C. Results and discussion

1.A Measurement of a linear polarizer

To evaluate the overall system performance we measured the Mueller matrix of a linear polarizer (a commercial grid polarizer) rotated at various orientations respect to the plane of incidence. The polarizer was placed after the PSG. Since the Mueller ellipsometer was set in the reflection mode, a gold mirror was placed on the sample-holder. The Mueller matrix of the polarizer was obtained multiplying the matrix of the system polarizer-mirror by the inverse of the matrix of the mirror measured without the polarizer. The results corresponding to a wavenumber of 1500 cm -1 are shown in Fig. 6. Spectra was recorded with 16 cm -1 resolution and 32 integrations. In those conditions, according to data in Fig. 5 and Fig. 6 the uncertainty due to random noise in Mueller matrix elements is expected to be approximately 0.0025. The absolute value of the rotation angle was known with an accuracy of about 0.3° for each measurement. To interpret the dependence of Mueller matrix elements with the polarizer azimuth, θ, we assumed that the polarizer behaved as an imperfect dichroic retarder with a small amount of depolarization. The imperfections caused a slight isotropic depolarization were represented by the parameter D. The model matrix could be written as: Apart from the azimuth and the depolarization, the matrix elements depend on the ellipsometric functions Ic', Is and Is which are, in turn, trigonometric functions of the ellipsometric angles Ψ and ∆. Experimental data was fitted to Eq. 6 model using as parameters Ψ, ∆ and D giving as a result: Ψ=87.4±0.5°, ∆=294.2±0.5°, D=0.971±0.001. We found from the analysis of residuals that the average value of the differences among fitted and experimental data was 0.008. According to Eq. 6, there are some pairs of Mueller elements that have identical respective values. Therefore, calculating the differences among them is an heuristic way to evaluate the magnitude of systematic and nonsystematic errors in experimental data. The lower panel of Fig. 6. shows the differences among three symmetric pairs of matrix elements. As can be seen in the lower panel of Fig. 6, there is a strong dependence of the values shown in there with the polarizer orientation. For azimuths comprised from -180° to 110° the values of the differences are constant and around 0.005, however, for azimuths comprised between -110° and -45°, the differences become higher, up to 0.01. In our opinion the parallelism between the faces of the polarizer was not perfect, thus causing beam wandering when the azimuth was modified. Beam wandering can create systematic errors because the beam is deviated from the path that initially followed during the calibration procedure. In our opinion for azimuths from -180° to -110° the magnitude of systematic errors is comparable or lower than that of random noise, but for azimuths between -110° and -45° systematic errors become dominant. Thus, in ideal conditions, a nominal accuracy of the instrument can be estimated to be about 0.005, which is the sum of two balanced contributions, noise bounded to 0.0025 and systematic errors bounded to about 0.0025. In non-ideal conditions, systematic errors can surpass random noise thus increasing the uncertainty of the measurements.

( ) ( ) ( ) ( )
′ ′ - -     ′ - + - -   =   ′ - - +     + -   (6)

2.Characterization of a thermally oxidized c-Si wafer

The second example consists of the measurement of a crystal silicon (c-Si) wafer covered with a thin film of amorphous silicon oxide (SiO2) thermally grown on it, with a nominal thickness of 1050 nm. The Mueller matrix of this kind of samples (isotropic and non-depolarizing) can be expressed in terms of the ellipsometric angles Ψ and ∆ as [START_REF] Goldstein | Polarized Light, 2 nd edition Revised and Expanded[END_REF]: As a representative example, Fig. 7 shows the experimental spectroscopic Mueller matrix at an angle of incidence of 68.2°. It can be seen that the measured matrix has the same structure as Eq. 7, with null block-diagonal elements. A common practice to extract information from ellipsometric data, is to simulate the optical response of the sample expressed in terms of Mueller matrix elements or any other related observable, by a theoretical model. To do so, typical models allow to calculate the reflectivity of the sample in polarization parallel "p" and perpendicular "s" to the plane of incidence respectively, when it has been illuminated by beams polarized according to "p" and "s" directions respectively. Once these 4 reflectivity coefficients are calculated, it is straightforward to calculate the Jones and the Mueller matrix of the considered sample [START_REF] Goldstein | Polarized Light, 2 nd edition Revised and Expanded[END_REF][START_REF] Jellison | Data Analysis for Spectroscopic Ellipsometry[END_REF]. There are well stablished methods to calculate the optical response of samples consisting of a number of layers deposited on a substrate [START_REF] Yeh | Optical Waves in Layered Media[END_REF]. In general, theoretical models depend on a number of parameters characteristic of the sample, which must be adjusted to make the theoretical data to ''fit'' the measurements. A typical model for a stack of layers includes the angle of incidence, the thickness and complex refractive indexes of all the layers. The dispersion of the complex refractive index can be taken from reference tables or represented by a parameterized mathematical expression. Further discussion about data interpretation and modelling in the case of either, standard ellipsometry and Mueller ellipsometry, is given in specialized monographs [START_REF] Goldstein | Polarized Light, 2 nd edition Revised and Expanded[END_REF][START_REF] Jellison | Data Analysis for Spectroscopic Ellipsometry[END_REF][START_REF] Fujiwara | Spectroscopic Ellipsometry. Principles and Applications[END_REF]. Ellipsometric information taken at different angles of incidence may be useful when it comes to fit measurements with a model that depends on many parameters. As discussed in previous works [START_REF] Irene | A brief History and State of the Art of Ellipsometry[END_REF][START_REF] Jellison | Data Analysis for Spectroscopic Ellipsometry[END_REF][START_REF] Jenkins | Multiple-angle-of-incidence ellipsometry[END_REF], simultaneously fitting data taken at different angles of incidence allows to decrease linear correlations among the fitted parameters and thus, to improve the accuracy of the fitted values. We took measurements at different incident angles, (from 75° to 65° with a step of 2°). Experimental data was fitted with a simple mathematical model consisting of a semi-infinite c-Si substrate covered by a homogeneous SiO2 layer. The dielectric functions of c-Si and SiO2 were taken from a reference [START_REF] Palik | Handbook of Optical Constants of Solids[END_REF]. The fitted parameters were the SiO2 layer thickness and angle of incidence. As previously mentioned, in conventional ellipsometry is common to use parameterized dielectric functions to represent the optical properties of the materials analyzed. We have preferred to work with reference tabulated values for two reasons: i) to reduce the number of parameters to be adjusted, and ii) because both, c-Si and SiO2 are extensively studied materials, and their dielectric functions are well-known. Based on our experience, we knew that that the tabulated dielectric functions for c-Si and SiO2 corresponded to those of the materials present in the sample. A detailed discussion of the parameterization of dielectric functions for c-Si and SiO2 goes beyond the objectives of this paper, but this topic has been discussed in detail by Ossikovski et al. in [START_REF] Ossikovski | Infrared ellipsometry study of the thickness-dependent vibration frequency shifts in silicon dioxide films[END_REF]. The resulting bestfitted values for the SiO2 layer thickness was 1046 ± 10 nm, in good agreement with the manufacturer specifications. The best fitted angles of incidence were 74.8°, 73.0°, 71.2°, 69.7°, 68.2° and 66.5°. Correlation between angle of incidence and thickness was low, about 0.3. To illustrate the dependence of the polarimetric data with the angle of incidence we plotted in Fig. 8 the experimental and best-fitted values. Because the block diagonal elements are null, and some of them are redundant, only the elements [START_REF] Irene | A brief History and State of the Art of Ellipsometry[END_REF][START_REF] Röseler | Spectroscopic Infrared Ellipsometry[END_REF], [START_REF] Korte | Infrared spectroscopic ellipsometry: A tool for characterizing nanometer layers[END_REF][START_REF] Schubert | Infrared Ellipsometry on Semiconductor Layer Structures[END_REF] and [START_REF] Schubert | Infrared Ellipsometry on Semiconductor Layer Structures[END_REF][START_REF] Schubert | Infrared Ellipsometry on Semiconductor Layer Structures[END_REF] are shown. An excellent agreement is obtained between experiment and theoretical values, providing the usefulness of the model and the accuracy of data measured with the ellipsometer. In the plots two types of spectral features can be seen. The first one, around 2500 cm -1 , is related to a Bragg type interference between beams partially reflected by the first and second faces of the SiO2 layer. The spectral position of this feature is mainly determined by the SiO2 layer thickness, but changing the angle of incidence makes the position of the interference to slightly shift by 75 cm -1 . Changing the angle of incidence modifies the path length of the beam inside the layer and thus the wavenumber for which the Bragg condition is fulfilled. Wavenumber (cm -1 ) Fig. 8. Spectroscopic Mueller elements [START_REF] Irene | A brief History and State of the Art of Ellipsometry[END_REF][START_REF] Röseler | Spectroscopic Infrared Ellipsometry[END_REF], [START_REF] Korte | Infrared spectroscopic ellipsometry: A tool for characterizing nanometer layers[END_REF][START_REF] Schubert | Infrared Ellipsometry on Semiconductor Layer Structures[END_REF] and [START_REF] Schubert | Infrared Ellipsometry on Semiconductor Layer Structures[END_REF][START_REF] Schubert | Infrared Ellipsometry on Semiconductor Layer Structures[END_REF] as a function of the angle of incidence. Arrows indicate the dependence of spectroscopic values with the angle of incidence.

The second spectral feature, between at 1000 and 1300 cm -1 , is related to the absorption due to the different stretching modes of the Si-O bond. In particular, the signatures of a transversal optical mode at 1130 cm -1 and a longitudinal optical mode at 1260 cm -1 are well visible [START_REF] Ossikovski | Infrared ellipsometry study of the thickness-dependent vibration frequency shifts in silicon dioxide films[END_REF]. The spectral position of this band, contrarily to that of the Bragg interference, does not depend on the angle of incidence.

D. Conclusions

We have built a new spectroscopic infrared Mueller ellipsometer based on achromatic optical elements. The optical set-up consists of a polarization states generator, a PSG, and a polarization states analyzer, a PSA, both being identical. The optical design is based on an objective criterion, the maximization of the condition number of the matrices representing the PSG and the PSA, which allows minimizing the propagation of random and systematic errors to the measurements. The eigenvalue method has been applied to calibrate the PSG and the PSA without ambiguities because it does not need to model the optical elements forming the Mueller ellipsometer. We have included in this article several examples to illustrate the quality of the polarimetric measurements, compatible with that of commercial ellipsometers.

Fig. 1 .

 1 Fig. 1. Schematic representation of the optical configuration of the infrared Mueller ellipsometer: (a) collimated beam configuration and (b) focused beam configuration.

Fig. 2 .

 2 Fig. 2. Schematic representation of the different movements needed to modify the angle of incidence.

Fig. 3 .

 3 Fig.3. Spectroscopic dependence of the condition number for experimental matrices W (blue dots) and A (red dots) in the range from 700 to 5000 cm -1 . Spectral resolution: 16 cm -1 . Number of averaged spectra: 16.

Fig. 5

 5 Fig. 5 Standard deviation as a function of the accumulation number, N. Spectral resolution was set to 8 cm -1 (blue) and to 16 cm -1 (red). Error bars represent the uncertainty in the determination of the standard deviation. The inset shows the intensity dependence of the intensity in a typical spectra with 8 cm -1 resolution.

Fig. 6 .

 6 Top panel Mueller matrix of a grid linear polarizer oriented at various azimuths. Dots correspond to experimental data at 1500 cm -1 and solid line to the best-fitted model curve. Lower panel. Differences among some symmetric Mueller matrix elements indicated in the legend. Lines are just a guide to the eye.

Fig. 7 .

 7 Fig. 7. Experimental (blue dots) and best-fitted (red line) spectroscopic Mueller matrices of a silicon wafer coated with a thermal SiO2 film of nominal thickness 1050 nm.

Acknowledgements

The authors wish to express their acknowledgment to Dr. Y. Bonnassieux, and Dr. P. Bradu for their help in finding funds and providing part of the components. Moreover, the authors wish to express their gratitude to Mr. D. Clément, Mr. C. Jadaud and Mr. F. Farci for their help in designing and building the nonstandard opto-mechanical pieces used. The development of the broadband Mueller ellipsometer was a project included in the Master's Program at Ecole Polytechnique. The authors, wish to thank Mr. R. Landon, supervisor of one the programs, for his help in communication and integration of students in the project.


Alphonse Lovamamy Rasoloniaina

Vincent Huet

Monique Thual

Stéphane Balac

Patrice Féron

Yannick Dumeige
email: yannick.dumeige@univ-rennes1.fr

Analysis of third-order nonlinearity effects in very high-Q WGM resonator cavity ringdown spectroscopy

Introduction

High Quality (Q) factor whispering gallery mode (WGM) resonators have drawn lots of interest in the last decades in the aim of optical and photonics applications [START_REF] Braginsky | Quality-factor and nonlinear properties of optical whispering-gallery modes[END_REF][START_REF] Vahala | Optical microcavities[END_REF][START_REF] Ilchenko | Optical resonators with whispering-gallery modes-part II: applications[END_REF][START_REF] Chiasera | Spherical whispering-gallery-mode microresonators[END_REF][START_REF] Righini | Whispering gallery mode microresonators: Fundamentals and applications[END_REF]. In the linear regime, their high spectral selectivity has been exploited to obtain high sensitivity biological or chemical sensors [START_REF] Vollmer | Protein detection by optical shift of a resonant microcavity[END_REF][START_REF] Ward | WGM microresonators: sensing, lasing and fundamental optics with microspheres[END_REF][START_REF] Zhu | On-chip single nanoparticle detection and sizing by mode splitting in an ultrahigh-Q microresonator[END_REF]. Their ultra-high-Q factor associated with a small size lead to a very high finesse F (up to F = 10 7 [START_REF] Savchenkov | Optical resonators with ten million finesse[END_REF][START_REF] Rasoloniaina | Controling the coupling properties of active ultrahigh-Q WGM microcavities from undercoupling to selective amplification[END_REF]) for relatively small mode volumes (V m ). This unique property induces very high intracavity optical intensity and makes WGM resonators very suitable for nonlinear photonics applications. In this framework, several miniaturized nonlinear optical devices whose operation relies on second or third order nonlinearities have been proposed and demonstrated: Raman and Brillouin lasers [START_REF] Spillane | Ultralow-threshold Raman laser using a spherical dielectric microcavity[END_REF][START_REF] Grudinin | Brillouin lasing with a CaF2 whispering gallery mode resonator[END_REF], hyperparametric generators [START_REF] Savchenkov | Low threshold optical oscillations in a whispering gallery mode CaF2 resonator[END_REF][START_REF] Matsko | Optical hyperparametric oscillations in a whispering-gallery-mode resonator: Threshold and phase diffusion[END_REF], frequency combs [START_REF] Del'haye | Optical frequency comb generation from a monolithic microresonator[END_REF][START_REF] Savchenkov | Tunable optical frequency comb with a crystalline whispering gallery mode resonator[END_REF][START_REF] Saha | Broadband parametric frequency comb generation with a 1 µm pump source[END_REF][START_REF] Del'haye | Self-injection locking and phase-locked states in microresonator-based optical frequency combs[END_REF], frequency converters [START_REF] Dumeige | Whispering-gallery-mode analysis of phase-matched doubly resonant secondharmonic generation[END_REF][START_REF] Andronico | Difference frequency generation in GaAs microdisks[END_REF][START_REF] Moore | Continuous-wave ultraviolet emission through fourthharmonic generation in a whispering-gallery resonator[END_REF][START_REF] Kuo | Secondharmonic generation using -quasi-phasematching in a GaAs whispering-gallery-mode microcavity[END_REF][START_REF] Farnesi | Optical frequency conversion in silica-whispering-gallery-mode microspherical resonators[END_REF][START_REF] Strekalov | Optical sum-frequency generation in a whispering-gallery-mode resonator[END_REF], optical parametric oscillators [START_REF] Fürst | Lowthreshold optical parametric oscillations in a whispering gallery mode resonator[END_REF][START_REF] Beckmann | Highly tunable low-threshold optical parametric oscillation in radially poled whispering gallery resonators[END_REF], quantum light generators [START_REF] Fürst | Quantum light from a whispering-gallery-mode disk resonator[END_REF] or optical switches and memories [START_REF] Treussart | Evidence for intrinsic Kerr bistability of high-Q microsphere resonators in superfluid helium[END_REF][START_REF] Ibrahim | Photonic logic NOR gate based on two symmetric microring resonators[END_REF][START_REF] Pöllinger | All-optical signal processing at ultra-lowpowers in bottle microresonators using the Kerr effect[END_REF]. Moreover, in some cases, WGM resonators can be fabricated without complex technological steps using melting and surface coating techniques [START_REF] Collot | Very high-Q whispering-gallery mode resonances observed on fused silica microspheres[END_REF][START_REF] Dong | High-Q silica microsphere by poly(methyl methacrylate) coating and modifying[END_REF][START_REF] Murzina | Kerr versus thermal nonlinear effects studied by hybrid whispering gallery mode resonators[END_REF] and may therefore be used as a simple benchmark for novel nonlinear materials. In this purpose, all the optical properties of the WGM resonator must be thoroughly measured to deduce accurately its nonlinear properties. To detail this statement, we consider the WGM resonator sketched in Fig. 1 with its single mode access line consisting of a tapered fiber [START_REF] Knight | Phase-matched excitation of whispering-gallery-mode resonances by a fiber taper[END_REF]. The input and output signals are respectively s in (t) and s out (t); u(t) is the resonator mode amplitude which is normalized in such a way that |u(t)| 2 is the energy stored in the cavity. The mode amplitude

WGM microcavity

Input signal s in (t) s out (t)

u(t)
Output signal τ 0 τ e

Tapered fiber

Fig. 1. Sketch of a WGM optical resonator coupled to a single access line. The input and output fields are sin(t) and sout(t), u(t) is the resonator mode amplitude. τ0 is the intrinsic photon lifetime and τe the coupling photon lifetime.

lifetime τ is related to the intrinsic (or optical loss limited) lifetime τ 0 and to the external (or coupling) lifetime τ e by: τ -1 = τ -1 0 +τ -1 e . Note that with this definition, τ is related to the overall Q-factor Q by Q = ω 0 τ /2 where ω 0 is the resonance angular frequency of the cavity. We consider a purely Kerr material with a refractive index given by:

N = N 0 + N 2 I (1)
where N 0 is the linear refractive index, I the optical intensity and N 2 the nonlinear refractive index given in m 2 /W. This resonator can have a bistable behavior for a threshold input power P th given by (see Appendix A):

P th = N 2 0 ω 0 V m 3 √ 3cQ 2 N 2 • τ e τ , (2)
we have assumed here for simplicity that V m ≈ V Kerr where V Kerr represents the effective mode volume characterizing the Kerr effect. As already stated, Eq. [START_REF] Vahala | Optical microcavities[END_REF] shows that low mode volumes and high Q-factors enhance the nonlinear effects and reduce the threshold power. Moreover, as it is the case for the linear dispersive properties of the resonator, its coupling regime is also very important. If we deduce the nonlinear refractive value from the measurement of the bistability threshold power, the value of the coupling strength between the resonator and the access line (through the parameter τ e ) play a crucial role as shown in Eq. ( 2). Finally, all the linear properties of a resonator: Q, τ e and V m are important to describe accurately its nonlinear operation. The Q-factor is usually deduced by measuring the full-width at half-maximum (FWHM) of the intensity transmission of the resonator or by cavity-ring-down techniques. The coupling properties or the relative value of τ 0 and τ e can be obtained by interferometric measurements [START_REF] Slagmolen | Phase-sensitive reflection technique for characterization of a Fabry-Perot interferometer[END_REF][START_REF] Vanier | Intrinsic quality factor determination in whispering gallery mode microcavities using a single Stokes parameters measurement[END_REF]. Dispersive and coupling properties can also be inferred by analyzing the ring-down signal [START_REF] Grudinin | Ultrahigh optical Q factors of crystalline resonators in the linear regime[END_REF] thanks to a simple linear coupled-mode model of the resonator when the probe frequency is rapidly swept across the resonance [START_REF] Dumeige | Determination of coupling regime of high-Q resonators and optical gain of highly selective amplifiers[END_REF]. Recently, it has been shown that this model can be developed to take into account the thermal effects [START_REF] Dong | Ringing phenomenon in silica microspheres[END_REF]. The mode volume can be deduced from the angular and radial structure of the WGM. This can be measured by mapping the light intensity using near field techniques [START_REF] Knight | Mapping whispering-gallery modes in microspheres with a near-field probe[END_REF][START_REF] Knight | Characterizing whispering-gallery modes in microspheres by direct observation of the optical standingwave pattern in the near field[END_REF][START_REF] Lin | Excitation mapping of whispering gallery modes in silica microcavities[END_REF] or the directional emission of asymmetric resonant WGM microspheres [START_REF] Gorodetsky | High-Q optical whispering-gallery microresonators: precession approach for spherical mode analysis and emission patterns with prism couplers[END_REF][START_REF] Dong | Directly mapping whispering gallery modes in a microsphere through modal coupling and directional emission[END_REF]. Note that in active WGM microspheres, the structure of the electric field can be inferred from the analysis of the fluorescence spectrum [START_REF] Lissillour | Whispering-gallery-mode laser at 1.56 µm excited by a fiber taper[END_REF][START_REF] Zhi | Structure of whispering gallery mode spectrum of microspheres coated with fluorescent silicon quantum dots[END_REF]. Furthermore, it has been already extensively reported that Rayleigh backscattering is enhanced in high finesse (and low mode volume) WGM resonators [START_REF] Ilchenko | Thermal nonlinear effects in optical whispering gallery microresonators[END_REF][START_REF] Weiss | Splitting of high-Q Mie modes induced by light backscattering in silica microspheres[END_REF][START_REF] Gorodetsky | Rayleigh scattering in high-Q microspheres[END_REF][START_REF] Kippenberg | Modal coupling in traveling-wave resonators[END_REF]. This phenomenon leads to a resonant coupling of the two co-and counter-propagating modes (often referred as modal-coupling) within the cavity which can be detected by measuring a reflected signal form the WGM resonator and a frequency splitting in the transmission spectrum [START_REF] Mazzei | Controlled coupling of counterpropagating whispering-gallery modes by a single Rayleigh scatterer: A classical problem in a quantum optical light[END_REF][START_REF] Li | Unified approach to mode splitting and scattering loss in high-Q whispering-gallery-mode microresonators[END_REF]. Cavity ring-down has been extended to take into account the beating between the two resonances.

In this case a thorough analyze of the ringing profile allows not only the Q-factor and the coupling regime but also the modal-coupling strength to be measured [START_REF] Trebaol | Ringing phenomenon in coupled cavities: Application to modal coupling in whispering-gallery-mode resonators[END_REF].

In this paper, we theoretically and experimentally show that by taking into account third order nonlinearities, the actual input power and the mode volume can also be inferred from the analysis of the transient response of the resonator when excited by a linearly chirped probe. Consequently we show that the cavityring-down spectroscopy method gives a full characterization of high-finesse WGM resonators under modalcoupling without the use of near field techniques or laser frequency locking to the cavity resonance [START_REF] Savchenkov | Low threshold optical oscillations in a whispering gallery mode CaF2 resonator[END_REF]. The paper is organized as follows. In section 2 we describe the studied microresonator configuration and the coupled mode theory (CMT) used to model its linear and nonlinear properties. The experimental setup is depicted in section 3. In a first step (section 4.A), the resonator is studied in the linear regime, which gives its dispersive properties and modal-coupling strength. The second step (section 4.B) consists in increasing the input power to reach the nonlinear regime of the resonator. Using the linear parameters deduced in the first step and comparing the calculations to the experimental results we show that it is possible to infer the actual input power P in = |s in | 2 and the effective mode volume V Kerr . Section 5 is devoted to a discussion regarding the consistency of the experimental results.

2. Structural description -Theory 2.A. CMT in the third order nonlinear regime Figure 2 is a sketch of the experimental configuration. The evanescent tail of the mode of a tapered fiber is side coupled to an high finesse microcavity WGM whose amplitude is u 2 (t). Due to the intracavity Rayleigh backscattering, the mode u 2 (t) is coupled to the counterpropagating mode u 1 (t) with a rate γ/2 [START_REF] Weiss | Splitting of high-Q Mie modes induced by light backscattering in silica microspheres[END_REF][START_REF] Kippenberg | Modal coupling in traveling-wave resonators[END_REF]. We assume that the two modes share the same cavity with the same resonance frequency (ω 0 = 2πc/λ 0 where λ 0 is the resonance wavelength) and same photon lifetimes (τ e , τ 0 ). For an input signal s in (t), the system have now two output ports giving a transmitted signal s out (t) and a reflected signal r out (t) with [START_REF] Trebaol | Ringing phenomenon in coupled cavities: Application to modal coupling in whispering-gallery-mode resonators[END_REF][START_REF] Haus | Coupled-mode theory[END_REF]:

       s out (t) = -s in (t) + 2 τ e u 2 (t) r out (t) = 2 τ e u 1 (t). (3)
In this model, we assume a good fidelity of the coupler which is the case for our experimental conditions (diam-eter of the tapered fiber ≈ 2 µm and Q-factor > 10 8 , see section 3) [START_REF] Spillane | Ideality in a fiber-taper-coupled microresonator system for application to cavity quantum electrodynamics[END_REF]. The output amplitudes can be obtained from the following differential system [START_REF] Trebaol | Ringing phenomenon in coupled cavities: Application to modal coupling in whispering-gallery-mode resonators[END_REF][START_REF] Dumeige | Stability and time-domain analysis of the dispersive tristability in microresonators under modal coupling[END_REF]:

         du 1 dt = j [ω 0 + ∆ω 1 (t)] - 1 τ u 1 + jγ 2 u 2 du 2 dt = j [ω 0 + ∆ω 2 (t)] - 1 τ u 2 + jγ 2 u 1 + 2 τ e s in (t) (4
) where ∆ω 1 and ∆ω 2 are the resonance frequency shifts due to nonlinear effects. These frequency shifts have two contributions: ∆ω therm due to thermal effects within the microcavity and ∆ω Kerr,i (for i ∈ {1, 2}) due to the Kerr effect:

∆ω i = ∆ω therm + ∆ω Kerr,i . (5)
Denoting by L the perimeter of the resonator, we have N ω 0 L/c = 2pπ where p ∈ N is the interference order. For a small change ∆L of the size L and ∆N of the refractive index N of the resonator the frequency shift ∆ω obeys:

- ∆ω ω 0 = ∆N N + ∆L L . (6)
Considering an increase ∆Θ of the temperature Θ of the resonator in the vicinity of the WGM, we can write the shift of the resonance frequency as [START_REF] Park | Regenerative pulsation in silica microspheres[END_REF]:

- ∆ω therm ω 0 = 1 N 0 ∂N 0 ∂Θ + 1 L ∂L ∂Θ ∆Θ = β∆Θ. (7)
Assuming that the intrinsic photon lifetime (τ 0 ) of the resonator is limited by a residual absorption α, the time variation of the temperature increase of the effective mode volume V ef f is such that [START_REF] Johnson | Self-induced optical modulation of the transmission through a high-Q silicon microdisk resonator[END_REF]:

ρV ef f C p d∆Θ dt source = αc N 0 |u 1 (t)| 2 + |u 2 (t)| 2 , (8)
where, ρ and C p are respectively the density and the heat capacity of the material of the resonator. Since τ 0 = 2N0 cα , introducing the effective thermal relaxation time 1/γ therm obtained from the thermal diffusion equation [START_REF] Fomin | Nonstationary nonlinear effects in optical microspheres[END_REF] and using the source term given in Eq. ( 8), we deduce the differential equation governing the time evolution of the temperature in the vicinity of the WGM [START_REF] Park | Regenerative pulsation in silica microspheres[END_REF]:

d∆Θ dt = -γ therm ∆Θ + 2 |u 1 (t)| 2 + |u 2 (t)| 2 τ 0 ρC p V Kerr . (9)
where we have assumed that V ef f ≈ V Kerr [START_REF] Schmidt | Nonlinear thermal effects in optical microspheres at different wavelength sweeping speeds[END_REF]. The thermal relaxation rate γ therm can be expressed using the thermal diffusivity D and the half thickness b of the mode in the largest field gradient by γ therm = D/b 2 [START_REF] Ilchenko | Thermal nonlinear effects in optical whispering gallery microresonators[END_REF]. This effective thickness can be estimated using the analytical results available for a well confined WGM for which b ≈ λ 0 ℓ 1/3 /(πN 0 ) [START_REF] Ilchenko | Thermal nonlinear effects in optical whispering gallery microresonators[END_REF] where ℓ ≈ N 0 L/λ 0 is the angular momentum number of the mode. For a single mode resonator, the Kerr effect induces a change of the refractive index ∆N = N 2 I. The optical intensity I is related to the electromagnetic energy density w = dU dV by:

I = c N 0 dU dV = c N 0 |u| 2 V Kerr . ( 10
)
Taking into account the cross Kerr effect between modes u 1 and u 2 , the resonance frequency shift due to the Kerr effect can be written as [START_REF] Matsko | Optical hyperparametric oscillations in a whispering-gallery-mode resonator: Threshold and phase diffusion[END_REF][START_REF] Yanik | Submicrometer all-optical digital memory and integration of nanoscale photonic devices without isolators[END_REF]:

- ∆ω Kerr,i ω 0 = N 2 c N 2 0 V Kerr |u i | 2 + 2 |u 3-i | 2 . ( 11
)
For a given input signal s in (t), integration of Eq. ( 4) and ( 9) using Eq. ( 11) and ( 7) provides the instantaneous values of the mode u 1 (t) and u 2 (t). The output signals are then deduced from Eq. ( 3). The amplitude transmission and reflection coefficients defined by x = s out /s in and y = r out /s in are related to the intensity transmission and reflection coefficients by T = |x| 2 and R = |y| 2 .

2.B. Linear approximation

Assuming that the input power |s in (t)| 2 is sufficiently weak, the response of the resonator can be assumed to be linear and thus ∆ω 1 = ∆ω 2 = 0. Denoting u ± = u 2 ±u 1 , the differential system (4) reads:

du ± dt = Λ ± u ± + 2 τ e s in (t) (12)
where Λ ± = j ω 0 ± γ 2 -1 τ . Integrating Eq. ( 12) we deduce the amplitude of the two modes u 1 (t) and u 2 (t) by the relations u 1 = (u + +u -)/2 and u 2 = (u + -u -)/2.

2.C. Stationary linear regime

In the stationary regime, the input signal can be written as s in (t) = s 0 e jωt . Solving Eq. ( 12) for this source term, we obtain the stationary modes:

u ± (t) = 2 τe s 0 e jωt j δ ∓ γ 2 + 1 τ , (13)
where we have introduced the quantity δ = ω -ω 0 . The amplitude transmission coefficient is then

x(δ) = jδ + 1 τ 1 τe -1 τ0 -jδ -γ 2 4 j δ + γ 2 + 1 τ j δ -γ 2 + 1 τ ( 14
)
whereas the amplitude reflected coefficient is:

y(δ) = j γ τe j δ + γ 2 + 1 τ j δ -γ 2 + 1 τ . ( 15
)
2.D. Frequency chirped input signal in the linear approximation We consider now an input signal with a carrier frequency linearly swept in time such that s in (t) = s 0 e jϕ(t) with:

ϕ(t) = ω i t + πV S t 2 , ( 16
)
where ω i is the angular frequency at t = 0 and V S is the frequency scanning speed. These quantities are related to the instantaneous frequency by:

ν(t) = 1 2π dϕ dt = ω i 2π + V S t. (17)
Integration of Eq. ( 12) requires to evaluate the indefinite integral

t 0 e [j(ωi-ω0∓ γ 2 +2πVS t ′ )+ 1 τ ]t ′ dt ′ = g ± (t) -g ± (0). ( 18
)
The function g ± can be analytically expressed using the complex error function erf [START_REF] Trebaol | Ringing phenomenon in coupled cavities: Application to modal coupling in whispering-gallery-mode resonators[END_REF]:

g ± (t) = - j 4V S exp -j 4πV S ω 0 -ω i ± γ 2 + j τ 2 ×erf ω 0 -ω i ± γ 2 -2πV S t + j τ √ 4πjV S ( 19
)
We obtain

u ± (t) = s 0 2 τ e 1 j ω i -ω 0 ∓ γ 2 + 1 τ + g ± (t) -g ± (0) × exp jω 0 ± j γ 2 - 1 τ t . (20)
Using Eq. ( 3), [START_REF] Dumeige | Whispering-gallery-mode analysis of phase-matched doubly resonant secondharmonic generation[END_REF] and [START_REF] Andronico | Difference frequency generation in GaAs microdisks[END_REF] we are in position to compute the dynamic responses of the cavity s out (t) and r out (t) for any frequency sweeping speed V S .

Experimental setup -Method 3.A. Experimental setup

We have used ultra-pure silica WGM microspheres as model systems to test the impact of the third order nonon the transient time domain response of high-Q microresonators. These microspheres are obtained by melting the tip of a first tapered fiber (Fig. 3.a) using a fiber fusion splicer (ERICSSON FSU 925). Using this process it is possible to obtain silica microspheres with diameters between 40 µm and 200 µm (Fig. 3.b). The experimental setup is sketched in Fig. 4.a). A second tapered fiber is used to inject the optical signal in the microsphere. This fiber taper has a diameter smaller than 2 µm over a length much longer than the microsphere diameter. Once the position of the microsphere is set, the fiber taper is mounted on 3-axis micro-positioning stages in order to control their relative distance (typically less than the wavelength). The probe laser is an external cavity laser diode in the C-band (wavelength λ ≈ 1550 nm) with a long-term linewidth of 150 kHz. The laser frequency is linearly swept (with a frequency scanning speed V S ) over 6 GHz using a 100 Hz triangle 18 V peak to peak signal. The signal power at the input of the tapered fiber is denoted P 0 (Fig. 4.b). Both transmitted and reflected signal are measured thanks to two amplified photodiodes (PD1 and PD2) and an optical circulator (C). The bandwidth of PD1 is 80 MHz and that one of PD2 is 14 MHz.

3.B. Method

For very low sweeping speed: V S ≪ V 0 where V 0 = 2/(πτ ) 2 corresponds to a frequency range 1/(πτ ) equal to the FWHM of the resonance scanned during a cavity lifetime τ /2, the response of the resonator is stationary and the recording of the transmission in time gives the spectral profile of the resonance as shown in Fig. 4.b).

In the linear regime, the Q-factor is deduced from i) the measurement of the FWHM of the experimental curve and ii) a calibration of the frequency sweeping speed. However for sweeping speeds such as V S V 0 , the dynamic response of the resonator can not be neglected. We have already demonstrated that in this case, using the model described in section 2.D it is possible to infer the Q-factor and the dispersive properties of a resonator without any frequency speed calibration [START_REF] Dumeige | Determination of coupling regime of high-Q resonators and optical gain of highly selective amplifiers[END_REF]. The experimental transmission signal T exp (t) is normalized using the off-resonance value. In the linear regime, the experimental signal is then compared to the theoretical model T theo (t) depending on τ 0 , τ e , γ and V S calculated from Eq. ( 3), ( 20) and [START_REF] Dumeige | Whispering-gallery-mode analysis of phase-matched doubly resonant secondharmonic generation[END_REF]. A nonlinear optimization procedure on τ 0 , τ e , γ and V S is used to minimize the value of σ 2 defined as

σ 2 (τ 0 , τ e , γ, V S ) = N k=1 [T exp (t k ) -T theo (t k , τ 0 , τ e , γ, V S )] 2 (21)
where N is the number of experimental sampling points. In the general case we have checked that the error on the measurements of τ 0 , τ e and γ was less than ±15 %. Note that the method also enables us to deduce the values of V S from the experiments. This value can be compared to its nominal value (depending on the peak to peak and the frequency of the electrical command signal) in order to check the consistency of our approach.

A first validation of the method is now presented using a same resonance of a 145 µm diameter microsphere for several sweeping speeds. τ 0 , τ e and γ are inferred using the procedure previously described from the transmission profiles. From We can notice a good agreement between the experiments and the calculations carried out from Eq. ( 14) and ( 15) using the linear parameters inferred from Fig. 5.c).

Full Characterization of a high-Q WGM resonator in the nonlinear regime

In this section we show that the nonlinear behavior of the resonator can be used to deduce other important parameters such as the resonator input power (P in ) and the mode volume. For this experimental demonstration we have used a silica microsphere with a diameter φ = 135 µm. For the same microsphere, we first deduce its linear properties using the method previously described. Then we increase the input power to characterize its nonlinear behavior.

4.A. WGM resonator linear property measurements

The first step of the method consists in measuring the linear parameters of the microsphere according to the method described in section 3.B. This is made using a backward frequency sweeping (V S < 0) and a forward frequency sweeping (V S > 0). In the nonlinear regime, the response signals obtained for V S < 0 and V S > 0 are different due to the nonlinear effects [START_REF] Savchenkov | Kilohertz optical resonances in dielectric crystal cavities[END_REF][START_REF] Upham | Enhancing optical field intensities in gaussian-profile fiber bragg gratings[END_REF]: i) for a V S < 0 the nonlinear frequency shift pushes the resonance away from the probe laser which leads to a broadening of the resonance; ii) for V S > 0, the nonlinear shift pulls the resonance towards the probe laser when the intracavity energy increases giving a resonance width smaller than in the linear case. In the linear regime, the frequency sweeping does not affect the response of the resonator whereas in the nonlinear regime this parameter strongly influences the shape of the transmission and reflection signals as it will be shown in the next section. Figure 6 shows the transmission and reflection signals for the 135 µm diameter microsphere and a power P 0 = 26.5 µW at the input of the taper. The two signal profiles are comparable for V S > 0 and Table 1. Linear parameters deduced from the fits of the transmission signals given in Fig. 6 V S < 0 which indicates that the input power is such that the linear regime is reached. The measurements deduced from the two fits in Fig. 6.a) and 6.b) are given in Tab. 1. The two sets of data are quite similar and the sweeping speeds are consistent with the nominal values (≈ 5 MHz/µs). The linear parameters deduced from Fig. 6.a) will be used in the nonlinear result analysis of section 4.B.

4.B. Experiments in the nonlinear regime

In a second step, we increase the input power to reach the nonlinear regime. In Fig. 7 we show the reflected and transmitted signals for two high nominal input powers P 0 = 210 µW and P 0 = 420 µW using both backward and forward frequency sweepings. The transmission and reflection signals can be simulated using the model described in section 2.A. An explicit embedded Runge-Kutta method of orders 4 and 5 (namely the Dormand-Prince pair) is used to numerically integrate Eqs. [START_REF] Chiasera | Spherical whispering-gallery-mode microresonators[END_REF]. All the parameters used for the calculations are given in Tab. 2. The only free parameters are the mode volume V Kerr (depending of the spatial structure of the excited WGM) and the actual input power P in at the waist of the tapered fiber. Simulation results given in Fig. 7 have been carried out by manually changing the values of V Kerr and P in in order to obtain a good Table 2. Physical parameters used in the calculations given in Fig. 7 and Fig. 8.

Parameter Value References N0 1.46 [57] β 8.83 × 10 -6 K -1 [64] N2 2.5 × 10 -16 cm 2 • W -1 [57] ρ 2.2 g • cm -3 [57, 59] Cp 670 J • kg -1 • K [59] D 9.5 × 10 -3 cm 2 • s -1 [59, 64]
agreement with the experimental results. We deduced from this procedure that V Kerr = (2.5 ± 0.2) × 10 4 µm 3 and that P in = (90 ± 8) µW for Fig. 7.a) and P in = (180 ± 15) µW for Fig. 7.b). Taking into account the taper losses (3.2 dB), the power at the input of the tapered fiber is respectively: P 0,calc = (190 ± 15) µW and P 0,calc = (380 ± 30) µW depending on the values of P in . It is in good agreement with the experimental data. The input power has been still increased to higher values (P 0 = 1.2 mW) as shown in Fig. 8. For V S < 0, the transmission signal is strongly broadened and the stiff part of the signal (obtained for low frequencies) is characteristic of a bistable operation. The model well reproduces the time domain variations of the transmitted and reflected signals for P in = (480 ± 40) µW or P 0,calc = (1.03 ± 0.8) mW close to the measured value. For V S > 0, an input power of P in = (1.2 ± 0.1) mW (or P 0,calc = (2.5 ± 0.2) mW) is required in the simulations to obtain calculated signals similar to experimental results. For moderate input powers, the model taking into account Kerr and thermal effects is well suited to describe the experimental results. For higher power, the discrepancy (observed in the reflected signal) between the measured and the input power deduced from calculations shows that other nonlinear effects (such as thermal effects within the taper, high-order nonlinearities in the microsphere, . . . ) have to be taken into account in the model.

Remarks and discussion

In Fig. 8.i) we also give calculations carried out for N 2 = 0, which are quite different from those taking into account the Kerr effect. This shows that both Kerr and thermal effects must be considered simultaneously as already reported for the stationary regime in Ref. [START_REF] Schmidt | Nonlinear thermal effects in optical microspheres at different wavelength sweeping speeds[END_REF].

We discuss now the value of the effective mode volume (V Kerr ) deduced from calculations. WGM propagating within a dielectric microsphere by total internal reflection can be understood as high angular momentum electromagnetic modes. They are labeled using a set of three quantum numbers (n, ℓ, m): n is the number of maxima of the field envelope in the radial direction, ℓ is the azimutal number (or the interference order) and m is the projection of the angular momentum on the quantiza- tion axis. For an ideal sphere, the resonance positions (λ 0 ) depend only on n and ℓ whereas ℓ-|m|+1 gives the number of maxima of the field envelope in the quantization direction. Thus the effective mode volume V Kerr defined by [START_REF] Braginsky | Quality-factor and nonlinear properties of optical whispering-gallery modes[END_REF][START_REF] Schmidt | Nonlinear thermal effects in optical microspheres at different wavelength sweeping speeds[END_REF] V

Kerr = R 3 |E(r)| 2 d 3 r 2 V |E(r)| 4 d 3 r (22)
(where E is the electric field and V is the physical volume of the microsphere) strongly depends on the values of these three numbers. The electric field E can be expressed in terms of Bessel and spherical harmonic functions [START_REF] Schiller | High-resolution spectroscopy of whispering gallery modes in large dielectric spheres[END_REF]. For an ideal sphere with φ = 135 µm and N 0 = 1.46, comparing the numerical values of the effective volume calculated using Eq. ( 22) to the value V Kerr = (2.5 ± 0.2) × 10 4 µm 3 inferred from the measurements we estimate that in the experiments reported in section 4.B the value of the three quantum numbers were: n = 2, ℓ = 377 and m = 371 for a resonant wavelength λ 0 ≈ 1550 nm. We identificated other sets of values for (n, ℓ, m) but they correspond to n > 2 and are experimentally less feasible using a tapered fiber as coupling system [START_REF] Knight | Characterizing whispering-gallery modes in microspheres by direct observation of the optical standingwave pattern in the near field[END_REF]. With these values, we estimate the mode volume by the following integral [START_REF] Collot | Very high-Q whispering-gallery mode resonances observed on fused silica microspheres[END_REF]:

V m = R 3 w(r)d 3 r max [w(r)] ( 23
)
where w is the electromagnetic energy density and find for the above mentioned set of quantum parameters V m = 9.6 × 10 3 µm 3 . Note that all the numerical results are given for the TE polarization, for TM polarization the difference in the mode volume is smaller than the uncertainty on the value inferred from the measurement. The frequency splitting can be estimated from the mode volume value by [START_REF] Weiss | Splitting of high-Q Mie modes induced by light backscattering in silica microspheres[END_REF]:

γ 2π = c λ 0 ρ sc α 2 V m , (24)
where ρ sc is the scatterer number density and α their average linear polarisability. For silica we have ρ sc α 2 = 2 × 10 -12 µm 3 [START_REF] Weiss | Splitting of high-Q Mie modes induced by light backscattering in silica microspheres[END_REF] and thus the frequency splitting is γ/(2π) = 2.8 MHz in very good agreement with the data inferred from the fitting curves carried out in the linear regime (see Tab. 1). Another expression of the frequency splitting can be found in Ref. [START_REF] Gorodetsky | Rayleigh scattering in high-Q microspheres[END_REF]. It is obtained by replacing V m by V Kerr in Eq. [START_REF] Strekalov | Optical sum-frequency generation in a whispering-gallery-mode resonator[END_REF]. In this case, with our data we find γ/(2π) = 1.7 MHz which is slightly different from our measurement. Nevertheless Note that in i) we have also added the simulations obtained for N2 = 0 in dash line.

we have shown that the method is self consistent giving a mode volume estimation compatible with both the modal coupling rate and the nonlinear third order effect power threshold.

Conclusion

We have proposed a simple model using the CMT to take into account third order optical nonlinear effects in high-finesse WGM resonators. This method allows cavity-ring-down under modal-coupling analysis since it is directly formulated in the time domain. In the linear regime we are able to deduce from cavity-ring-down spectroscopy signals the coupling regime, the Q-factor and the mode-coupling strength. We have shown that the method can be extended to the nonlinear regime and that enables the input power and the mode volume to be estimated for silica microspheres. The mode volume deduced from nonlinear experiments is consistent with estimations made from the frequency splitting due to the mode coupling resulting from enhanced Rayleigh backscattering. Consequently, this shows that for a given material the combination of measurement of backscattering strength and nonlinear effects power threshold can be used to deduce the nonlinear index of this material at the cost of the knowledge of several physical parameters (thermal and thermo-optic properties; scatterers density and polarisability). This could be applied to determine the nonlinear properties of novel optical materials [START_REF] Murzina | Kerr versus thermal nonlinear effects studied by hybrid whispering gallery mode resonators[END_REF].

Fig. 3 .

 3 Fig. 3. Microsphere fabrication process. a) A 125 µm silica wire is drawn to obtain a tapered fiber. b) The tip of the fiber is melted using a fiber fusion splicer.

Fig. 4 .

 4 Fig. 4. a) Experimental setup. ISO: optical isolator, VA: variable attenuator, PC: polarization controller, C: optical circulator, PD1 and PD2: amplified photodiodes. P0: injected optical power at the input of the tapered fiber, Pin: optical power at the input of the WGM resonator, T (t) and R(t): respectively the transmitted and reflected signals. b) Sketch of a typical transmission signal and the corresponding instantaneous frequency variations of the probe laser.

Fig. 5 .

 5 Fig. 5. Transmission and reflection signals measured for a 145 µm diameter silica microsphere for several nominal frequency sweeping speeds: a) and b) VS = 5 MHz/µs; c) and d) VS = 1 MHz/µs; e) and f) VS = 0.3 MHz/µs. The calculation given in b) has been carried out using the parameters deduced from the fit of a). The theoretical curves shown in d)-f) have been calculated using the parameters deduced from the fit of c).

Fig. 5 .

 5 a) we deduce the following parameters values: τ 0 = 360 ns, τ e = 3.2 µs, γ/(2π) = 3.2 MHz and V S = 4.7 MHz/µs whereas from Fig. 5.c) we deduce: τ 0 = 480 ns, τ e = 4.3 µs, γ/(2π) = 3.3 MHz and V S = 1.2 MHz/µs. First, we can observe that the values of the sweeping speed deduced from the fits are very close from the nominal values. These two measurements give a mean value of the cavity lifetime τ = (380 ± 50) ns leading to V 0 = 4.4 MHz/µs and a Q-factor around 2 × 10 8 . Using these parameters we can simulate the reflected signals which are in good agreement with the measured signals as shown in Fig. 5.b) and 5.d). Finally, the transmission and reflection signals have been measured for a low sweeping speed V S = 0.3 MHz/µs verifying V S ≪ V 0 . In this case we obtain the stationary profile of the resonator response, moreover the abscissa axis can be converted in frequency using the value of V S as shown in Fig 5.e) and 5.f).

Fig. 6 .

 6 Fig. 6. Transmission and reflection signals in the linear regime measured for a 135 µm diameter silica microsphere for a nominal frequency sweeping speeds |VS| = 5 MHz/µs: a) and c) for VS < 0; b) and d) VS > 0. The theoretical reflection signals have been deduced from the fits of the transmission signals.

  .a) and Fig. 6.b). Sweeping τ0 τe γ/(2π) |VS| Figure VS < 0 405 ns 1.01 µs 2.5 MHz 4.2 MHz/µs 6.a) VS > 0 370 ns 1.04 µs 2.7 MHz 4.8 MHz/µs 6.b)

Fig. 7 .

 7 Fig. 7. Left columns: experimental transmission and reflection signals measured for backward and forward frequency sweepings in the nonlinear regime for two different optical powers P0 at the input of the tapered fiber. Right columns: theoretical calculations obtained for a) Pin = (90 ± 8) µW and b) Pin = (180 ± 15) µW with an effective mode volume VKerr = (2.5 ± 0.2) × 10 4 µm 3 .

Fig. 8 .

 8 Fig. 8. Left column: experimental transmission and reflection signals measured for backward and forward frequency sweepings for an optical power P0 = 1.2 mW at the input of the tapered fiber. Right column: theoretical calculations obtained for i) Pin = (480±40) µW and ii) Pin = (1.2±0.1) mW with an effective mode volume VKerr = (2.5±0.2)×10 4 µm 3 .Note that in i) we have also added the simulations obtained for N2 = 0 in dash line.

WGM microcavity + Rayleigh backscaterring Input signal

		τ 0
		γ
		u 1 (t)
		u 2 (t)	τ 0
	s in (t) r out (t)	τ e	s out (t)
	Reflected signal	Tapered fiber	Transmitted signal
	Fig. 2. Sketch of a WGM optical resonator under modal-
	coupling. The input signal sin(t) excites the mode u2(t)
	which is coupled to the counter-propagating mode u1(t)
	thanks to Rayleigh backscattering. The two modes have the
	same intrinsic and coupling photon lifetimes. The mode u1(t)
	generates a reflected signal rout(t) from the WGM resonator.

Acknowledgments

We acknowledge G. Nunzi Conti for providing us with the tapered fiber. This research was supported by the Agence Nationale de la Recherche through Projects ORA (ANR 2010 BLAN-0312), CALIN (ANR 2010 BLAN-1002) and by the Centre National d'Etudes Spatiales through Action R&T SHYRO R-S10/LN-0001-004. A.R. acknowledges support from the Conseil Régional de Bretagne (ARED). V.H. acknowledges support from the Direction Générale de l'Armement (DGA). Y.D. is member of the Institut Universitaire de France.

Appendix A: Bistability power threshold For a single mode resonator without modal coupling the mode amplitude dynamic can be deduced using Eqs. ( 4) and ( 11)

In the stationary regime at angular frequency ω, Eq. ( 25) can be analytically solved and we have with the notations introduced in section 2.C

Finally, by injecting the expression of |u th | 2 in Eq. ( 26) we find the threshold input power given by Eq. ( 2). This expression is exactly that of the hyperparametric threshold [START_REF] Matsko | Optical hyperparametric oscillations in a whispering-gallery-mode resonator: Threshold and phase diffusion[END_REF] which shows that the knowledge of the coupling properties is useful for the thorough analysis of several third order nonlinear processes.


Jérôme Fort
email: fort.jerome@gmail.com

Yves Cherel

Ann M A Harding

Carsten Egevang

Harald Steen

Grégoire Kuntz

Warren P Porter

David Grémillet

The feeding ecology of little auks raises questions about winter zooplankton stocks in North Atlantic surface waters

Keywords: diving behaviour, energetic modelling, food requirements, marine food web, stable isotope, seabird

Running head: seabirds and zooplankton stocks

Copepods are essential components of marine food webs worldwide. In the North Atlantic, they are thought to perform vertical migration and to remain at depths >500 m during winter. We challenge this concept through a study of the winter feeding ecology of little auks (Alle alle), a highly abundant planktivorous seabird from the North Atlantic. By combining stable isotope and behavioural analyses, we strongly suggest that swarms of copepods are still available to their predators in water surface layers (<50 m) during winter, even during short daylight periods. Using a new bioenergetic model, we estimate that the huge number (20-40 million birds) of little auks wintering off southwest Greenland consume 3600-7200 tonnes of copepods daily, strongly suggesting substantial zooplankton stocks in surface waters of the North Atlantic in the middle of the boreal winter.

INTRODUCTION

Copepods are essential components of aquatic food webs and play an important role in marine ecosystems. Their life strategies have been extensively studied (e.g. Falk-Petersen 2009) and their ecology is considered to be well known. At high latitudes and during winter, copepods are thought to perform vertical migration to deeper water to undergo diapause, and to become unavailable to most predators (Falk-Petersen 2009). However, a recent study has shown that in high Arctic waters copepods maintain their synchronized diel vertical migration throughout the polar night [START_REF] Berge | Diel vertical migration of Arctic zooplankton during the polar night[END_REF], suggesting that further investigations are necessary to yield a complete view of copepod winter ecology.

Seabirds are powerful indicators of marine food webs (e.g. [START_REF] Piatt | Seabirds as indicators of marine food supplies: Cairns revisited[END_REF]). Here we used little auks (Alle alle) as biosamplers of copepods, since this small Arctic species is the only seabird in the North Atlantic which feeds almost exclusively on zooplankton, essentially calanoid copepods [START_REF] Stempniewicz | Little Auk Alle alle[END_REF]. We used isotopic tracers [START_REF] Kelly | Stable isotopes of carbon and nitrogen in the study of avian and mammalian trophic ecology[END_REF] and Time-Depth-Recorders to investigate little auk diet and foraging depths during winter. It is the first study to record the winter diving behaviour of such a small seabird. We tested the hypothesis that copepods are still available in upper surface layers and consumed by little auks, even during the short winter daylight period.

MATERIAL AND METHODS

Little auks breed at high latitudes around the Arctic Ocean. The North-East Atlantic populations winter south of the ice edge, mainly off southwest Greenland and Newfoundland [START_REF] Stempniewicz | Little Auk Alle alle[END_REF]. Twenty breeding little auks were caught during summer 2007 at Kap Höegh (East Greenland; 70°44'N 21°35'W). Blood samples (~0.2ml) were collected from the brachial vein, as well as 16 chick meals transported in adult's sublingual pouch. Twenty wintering birds were legally shot at sea off Nuuk (southwest Greenland; 64°10′N 51°45′W) in January 2007 and blood samples collected from the cardiac clot. Stomach contents were removed, but digestion precluded identifying the few prey items.

To investigate copepod isotopic signatures, only chick meals containing >90% copepods were analysed. Stable isotope analyses were performed following [START_REF] Cherel | Stable isotopes document seasonal changes in trophic niches and winter foraging individual specialization in diving predators from the Southern Ocean[END_REF], with blood  15 N ( 15 N/ 14 N) values reflecting the trophic position / diet [START_REF] Kelly | Stable isotopes of carbon and nitrogen in the study of avian and mammalian trophic ecology[END_REF]. Values are mean ± S.D.

Eighteen breeding little auks from Longyearbyen (Spitsbergen; N78°13'N 15°20'E) and 22 from Kap Höegh were equipped with a Time-Depth-Recorder archival tag (G5, CEFAS Technology Limited, Lowestoft, UK) during summers 2007 and 2008, respectively. Recorders were implanted into the abdominal cavity following [START_REF] Grémillet | Cormorants dive through the Polar night[END_REF]. Each batch of loggers was programmed to record pressure during 24 hour periods at a sampling rate of 2 seconds; every 12 th day from 6 th Stored data were analysed using MultiTrace-Dive software (Jensen Software Systems). To avoid artefacts due to waves or tag accuracy, only dives ≥1.5 m were analyzed. To compare winter results of behavioural and isotopic analyses, only dive data recorded in December and January were used.

Individual daily energy and food requirements during December-January were calculated using the bioenergetics model Niche Mapper™ [START_REF] Fort | Thermodynamic modelling predicts energetic bottleneck for seabirds wintering in the northwest Atlantic[END_REF]). The principles of Niche Mapper™, its validation and all input data are detailed in the electronic supplementary material (ESM-1).

Given the isotopic results (see below), birds were assumed to feed on copepods during winter. We then estimated the daily food intake for the whole little auk population wintering off southwest Greenland. This population is composed to a large extent of birds breeding in northwest Greenland, east Greenland and Svalbard [START_REF] Stempniewicz | Little Auk Alle alle[END_REF]. The number of little auks wintering off southwest Greenland is therefore estimated as 20-40 million birds [START_REF] Barrett | Seabird numbers and prey consumption in the North Atlantic[END_REF][START_REF] Karnovsky | Estimation of carbon flux to dovekies (Alle alle) in the North Water[END_REF].

RESULTS

During summer, the δ 15 N value of chick meals (which correspond to copepod prey, [START_REF] Harding | Estimating prey capture rates of a planktivorous seabird, the little auk (Alle alle), using diet, diving behaviour, and energy consumption[END_REF]) averaged 8.4±0.5‰ (figure 1). Summer and winter bird blood δ 15 N values were almost identical (figure 1), yet statistically different (11.4±0.1 and 11.7±0.4‰, respectively; t-test, t=3.75, p<0.001).

During summer, adult bird blood was 3‰ enriched in δ 15 N compared to copepod prey, as expected for an increase of one trophic level (figure 1). During winter, measured bird blood was also 3‰ enriched in δ 15 N compared to copepods (figure 1, [START_REF] Sato | Stable isotopic compositions of overwintering copepods in the arctic and subarctic waters and implications to the feeding history[END_REF].

During winter (December-January), the mean depth of dives was 12±11 m, with a maximum dive depth of 50 m (n=6056 dives). The mean dive duration was 45±25 s. On average, birds spent 24% of their time diving, an activity which was equally performed during the daylight and night periods (48% and 52% of time diving, respectively) (figure 2). Dives were mostly concentrated during the daylight period and the first part of the night (figure 2; detailed results in ESM-2).

During January, Niche Mapper™ predicted an individual daily energy expenditure of 438±41 kJ•day -1 , and a food requirement necessary to meet this expenditure of 180±17g zooplankton•day -1 (wet mass). Consequently, we estimate that the southwest Greenland wintering population of 20-40 millions little auks consumes 3600-7200 tonnes of zooplankton per day (wet mass).

DISCUSSION

Using a multidisciplinary-approach combining stable isotope analysis, biotelemetry and energetic modelling, this study presents indirect evidence that during winter daylight periods copepods are still numerous in surface waters <50m off southwest Greenland. Little auks are known to feed almost exclusively on copepods in summer (e.g. [START_REF] Karnovsky | Seasonal changes in diets of seabirds in the North Water Polynya: a multiple-indicator approach[END_REF], and blood δ 15 N values found here are in agreement with birds preying upon copepods during this period [START_REF] Karnovsky | Seasonal changes in diets of seabirds in the North Water Polynya: a multiple-indicator approach[END_REF]. Isotopic δ 15 N signatures are also similar for summer and winter copepods [START_REF] Sato | Stable isotopic compositions of overwintering copepods in the arctic and subarctic waters and implications to the feeding history[END_REF], figure 1), as well as for summer and winter adult blood samples (figure 1). This strongly suggests that little auks feed predominantly on copepods in winter. During this period, birds may also consume other prey, but in a minor proportion since amphipods, euphausiids or fish larvae have higher isotopic signatures than copepods [START_REF] Tamelander | Pelagic and sympagic contribution of organic matter to zooplankton and vertical export in the Barents Sea marginal ice zone[END_REF]) and would thus lead to higher blood δ 15 N values. Winter blood values slightly higher than those measured in summer might reflect that amphipod/euphausiid effect, as well as the general food web enrichment from summer to winter [START_REF] Rolff | Seasonal variation in δ 13 C and δ 15 N of size-fractionated plankton at a coastal station in the northern Baltic proper[END_REF].

Further, we infer from little auk diving behaviour that copepods remain abundant in winter surface waters, even during the daylight period when they are supposed to perform vertical migration to greater depths (figure 2, [START_REF] Fortier | Visual predators and the diel vertical migration of copepods under Arctic sea ice during the midnight sun[END_REF][START_REF] Berge | Diel vertical migration of Arctic zooplankton during the polar night[END_REF]. Even though our dive depth results are consistent with summer recordings [START_REF] Harding | Estimating prey capture rates of a planktivorous seabird, the little auk (Alle alle), using diet, diving behaviour, and energy consumption[END_REF], further studies are now needed to confirm the first winter dive data for little auks. Their performance is impressive, with hundreds of dives conducted daily, over 1 minute and to depths of up to 50 m (see ESM-2). However, this maximum dive depth is extremely shallow compared to the depths of >500 m at which copepods are assumed to occur in winter.

Copepods are classically thought to enter diapause in the fall and to spend all winter in deeper water (Falk-Petersen 2009). During early spring, and mostly synchronized with the phytoplankton spring bloom, they are supposed to migrate back to surface waters (Falk-Petersen 2009). However, a recent acoustic study has shown that even in winter and during the high-latitude polar night, some copepods perform a synchronized diel vertical migration in the epipelagic layer [START_REF] Berge | Diel vertical migration of Arctic zooplankton during the polar night[END_REF], whereby they migrate to the food-rich surface layers during darkness and move to deeper water during daylight where they are out of reach from predators [START_REF] Fortier | Visual predators and the diel vertical migration of copepods under Arctic sea ice during the midnight sun[END_REF]. Our study strongly suggests that during winter, swarms of copepods are still present in surface waters (sensu [START_REF] Berge | Diel vertical migration of Arctic zooplankton during the polar night[END_REF], where they are targeted by little auks both during the daylight and the dark period.

Using a bioenergetic model, we estimated that the amount of copepods available to little auks wintering off southwest Greenland is substantial, since the daily intake of the little auk population wintering off southwest Greenland is estimated to amount to 3600-7200 tonnes. This result is based on a diet composed only of copepods and it might be slightly overestimated since little auks probably also consume a few other prey items such as amphipods and euphausiids. The standing stock of copepods is currently unknown and it is impossible to estimate the proportion of the population caught per day by wintering birds. However, a range of 3600-7200 tonnes of copepods caught per day during the winter off southwest Greenland is larger than the daily amount of copepods consumed during summer by the breeding population of the North Water Polynya (i.e. 3450-6900 tonnes of wet food consumed per day, including 2340-4680 tonnes of copepods; [START_REF] Karnovsky | Estimation of carbon flux to dovekies (Alle alle) in the North Water[END_REF]. Therefore, this suggests that the daily presence of copepods in the upper layer of the water column throughout the winter is much more important than usually assumed, with far-ranging consequences for regional food webs and energy flow.

  November 2007 to 22 nd February 2008 and every 10 th day from 1 st October 2008 to 10 th May 2009. Five birds from Longyearbyen were resighted in July 2008 and three of those were recaptured. Only one logger successfully recorded pressure. Thirteen birds from Kap Hoëgh were resighted and recaptured in July 2009. Only two loggers successfully recorded pressure in December-January.

Figure 1 .

 1 Figure 1. Stable nitrogen isotope values for chick meals (n=9) and for little auks blood samples (n=20)

Figure 2 .

 2 Figure 2. Average number of dives ( ) and maximum depth reach for any dive ( ) recorded at each

Acknowledgments

We thank all assistants for their hard work in the field. This project was supported by the Institut Polaire Français Paul Emile Victor (Grant 388 to D.G.) and the National Science Foundation (Grant 0612504 to A.H. & N.Karnovsky). Fieldwork led in Greenland was conducted under permits of the IPEV Ethics Committee (MP/12/24/05) and with the permission of the Greenland Home Rule Government (Danish Polar Center permit 512-240). Fieldwork led in Spitsbergen was conducted under permits 234 Forsøksdyrutvalget, Mattilsynet, Norway. J.F. is supported by a grant from the University Louis Pasteur (Strasbourg, France).