Chen Xie

Vytautas Jukna

Carles Milian

Remo Giust

Ismail Ouadghiri Idrissi

Tatiana Itina

John Michaël Dudley

Arnaud Couairon

François Courvoisier

Tubular filamentation for laser material processing

teaching and research institutions in France or abroad, or from public or private research centers.

Tubular filamentation for laser material processing

Chen Xie [START_REF] Allen | Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes[END_REF] , Vytautas Jukna [START_REF] Grier | A revolution in optical manipulation[END_REF][START_REF] Fürhapter | Spiral interferometry[END_REF] , Carles Milián [START_REF] Grier | A revolution in optical manipulation[END_REF] , Remo Giust [START_REF] Allen | Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes[END_REF] , Ismail Ouadghiri [START_REF] Allen | Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes[END_REF] , Tatiana Itina [START_REF] Fürhapter | Spiral interferometry[END_REF] , John M. Dudley [START_REF] Allen | Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes[END_REF] , Arnaud Couairon [START_REF] Grier | A revolution in optical manipulation[END_REF] and Francois Courvoisier (1,*) In 1992, Allen et al demonstrated that light beams with a phase singularity or vortex charge (i.e. orbital angular momentum) can transfer mechanical torque [START_REF] Allen | Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes[END_REF] . These waves attracted intense interest for optical manipulation of micro and nano-objects [START_REF] Grier | A revolution in optical manipulation[END_REF] . The vortex singularity was identified as a supplementary degree of freedom for optical encoding and applications emerged in a diversity of fields such as optical metrology [START_REF] Fürhapter | Spiral interferometry[END_REF] , nonlinear optics [START_REF] Firth | Optical solitons carrying orbital angular momentum[END_REF][START_REF] Courtial | Second harmonic generation and the conservation of orbital angular momentum with high-order Laguerre-Gaussian modes[END_REF][START_REF] Soljačić | Integer and fractional angular momentum borne on self-trapped necklace-ring beams[END_REF] , and quantum information [START_REF] Mair | Entanglement of the orbital angular momentum states of photons[END_REF][START_REF] Allen | Optical Angular Momentum[END_REF] .

In this paper, we report on an important novel application of vortex beams to control nonlinear propagation of intense femtosecond light beams within dielectrics, and demonstrate its relevance for laser material processing. In particular, for the first time to our knowledge, ultra-intense light pulses are shown to be able to propagate as a "light tube" in dielectrics without deformation, generating a tubular plasma distribution that reaches optical breakdown densities. This novel tubular geometry for delivery of energy from femtosecond pulses in a propagation-invariant way is expected to generate novel breakthroughs in femtosecond laser material processing.

The complex nonlinear dynamics of femtosecond laser pulse propagation in dielectrics generally inhibits uniform energy deposition within the material along the beam path.

The addition of a vortex singularity on a Gaussian beam in the form of a doughnut has been shown to introduce some additional propagation control [START_REF] Desyatnikov | Optical Vortices and Vortex Solitons[END_REF] , associated with azimuthal modulation instability, filament formation and soliton dynamics [START_REF] Vuong | Collapse of optical vortices[END_REF][START_REF] Maleshkov | Filamentation and supercontinuum generation by singular beams in self-focusing nonlinear media[END_REF][START_REF] Polynkin | Self-Focusing of Ultraintense Femtosecond Optical Vortices in Air[END_REF] . Several solitonic regimes have been identified, where the beam can either remain doughnut shaped or split into several rotating filaments [START_REF] Desyatnikov | Optical Vortices and Vortex Solitons[END_REF] . However, this solitonic behavior which would be essential for material processing stops at higher intensity where strong nonlinear losses occur.

In contrast to Gaussian beams, however, nondiffracting Bessel beams have been successfully shown to sustain invariant propagation even in presence of nonlinear losses and Kerr effect [START_REF] Porras | Nonlinear unbalanced Bessel beams: stationary conical waves supported by nonlinear losses[END_REF] . This arises physically because of the particular property of Bessel beams where conical inward energy flow can compensate for losses and "selfheal" the beam during propagation [START_REF] Polesana | Filamentation in Kerr media from pulsed Bessel beams[END_REF] . Experimental results have applied this propagation-invariance of zero-order femtosecond Bessel beams to single shot laser processing of nanochannels in glass with high aspect ratios from 100:1 to 1000:1 [START_REF] Bhuyan | High aspect ratio nanochannel machining using single shot femtosecond Bessel beams[END_REF][START_REF] Bhuyan | Single shot high aspect ratio bulk nanostructuring of fused silica using chirped controlled ultrafast laser Bessel beams[END_REF] .

In this paper, we investigate nonlinear femtosecond Bessel beams carrying a vortex charge. The concept is shown in Figure 1. The beam structure consists of a primary "light tube" of high intensity surrounded by several other concentric tubes at lower intensity. This is shown in Figure 1(a). With linear propagation, the intensity distribution is quasi-invariant along z [START_REF] Paterson | Higher-order Bessel waves produced by axicon-type computer-generated holograms[END_REF] . The radial profile of Bessel beams carrying a vortex charge m is well approximated by the Bessel function of order m:

2 ( ) ( sin ) m I r J k r  
where k is the wavevector and  is the conical angle. As will be discussed below, the conical angle is a parameter of primary importance. For 0 m  , the tube radius is null and the conical angle is the angle that geometrical rays make with the optical axis. For 1 m  , geometrical analysis shows that light propagates tangentially to the tube and no light crosses the inner part of the tubular main intensity region [START_REF] Berry | Exact and geometrical optics energy trajectories in twisted beams[END_REF] . The experimental set up used to generate the beams and image their propagation is shown in Figure 1(b) and this is discussed in more detail later. This paper is organized as follows. We first describe propagation-invariant vortex solutions to the nonlinear Schrödinger propagation equation, including the Kerr effect and nonlinear losses. We then analyze with numerical simulations the domain of existence of such solutions that we refer to as propagation-invariant conical vortex waves, and we see how nonlinear propagation is attracted toward these propagationinvariant solutions. To experimentally confirm our numerical results, we needed an accurate technique to record beam propagation. Spatiotemporal measurements of nonlinear propagation have already been developed within liquids by using a cuvette with variable thickness [START_REF] Jarnac | Whole life cycle of femtosecond ultraviolet filaments in water[END_REF] . Here, for propagation within solids, we have developed a novel experimental approach. Based on beam scanning, our setup allows for a direct imaging of the 3D fluence distribution within transparent solids in near and far field.

Quantitative agreement is found with numerical simulations. We discuss our results in terms of three broad regimes of propagation: (i) Tubular propagation-invariant where the intensity distribution pattern does not vary with the propagation distance, (ii) rotating, where several filaments rotate around the optical axis during the propagation and (iii) speckle-like, where hot-spots are non-rotating, and appear and disappear in space.

Results

Propagation-invariant solutions to nonlinear propagation in dielectrics

The key physical effects in the physics of stationary conical vortex waves are diffraction, the optical Kerr effect and nonlinear losses. The nonlinear Schrödinger equation (NLSE) describing the propagation of a field envelope E in dielectrics is:

2 2 2 2 2 0 0 0 2 2 K K n E i E ik E E E E z k n          (1)
where 0 k is the wavevector, 2   is the transverse Laplacian operator, 0 n and 2 n are respectively the refractive and Kerr indices of the medium, and K  is the cross section for multiphoton ionization.

We seek propagation-invariant monochromatic solutions in the form of a vortex of charge m , such that the electric field amplitude is: ( , , ) ( )exp( ( )

) z E r z a r ik z i r im        .
( ) a r and ( ) r  are to be determined while the longitudinal component of the wavevector,

0 cos z k k  
, is a parameter that can be arbitrarily chosen. Our mathematical procedure is summarized in the methods section and more details are provided in reference [20]. It follows a similar approach as in references [13,21] We found numerically that there is a limit in maximal intensity for each cone angle and vortex charge, over which no propagation-invariant solution is found. In Figure 2(b), we show in white (resp. white and light grey) the domain of existence of propagation invariant solutions in the parameter space determined by peak intensity and conical angle for a vortex charge m=3 (resp. m=1) in fused silica. The location of the frontier (red for m=3 and blue for m=1) obviously depends on material parameters and vortex order but remains approximately linear with positive slope.

The error bars shown correspond to the accuracy to which numerical integration could locate the frontiers.

This shows that increasing cone angle allows stationary solutions to exist in a broader range of intensities where propagation-invariant solutions can be found. More importantly, this result is qualitatively valid for any kind of nonlinear losses, including those occurring by plasma absorption. This is a major result for the effectiveness of propagation-invariant conical vortex waves to applications requiring high peak powers such as femtosecond laser micromachining.

As with Bessel beams, ideal propagation-invariant conical vortex waves are only weakly localized, i.e., their amplitude tails decay as 1 / r and carry infinite power. As will be shown below, a finite power beam carrying a suitably designed spatial phase will reshape into apodized versions of propagation-invariant conical vortex waves.

Finite energy solutions and experimental results

Here, we compare experimental and numerical propagation of finite energy nonlinear Bessel vortices with the propagation-invariant solutions found above.

Numerical model

We produce higher order Bessel beams from a Gaussian beam by using a phase mask

( , ) sin r k r m       17
. This is equivalent to the phase applied by an axicon and a vortex phase plate of order m . The numerical propagation model is detailed in the methods section. Briefly, it takes the canonical form of a unidirectional envelope propagation equation (UPPE) written in the spectral domain [START_REF] Couairon | Practitioner's guide to laser pulse propagation models and simulation[END_REF] , coupled to a plasma equation including photoionization, avalanche and recombination (see methods section). In our conditions, temporal variations are negligible and for computational efficiency in 3D, we consider the electric field ( , , ) E x y z as independent of time.

Experimental setup

Our experimental setup is shown in figure 1 Experimentally recorded fluences exceed 3 J/cm 2 , well above material damage threshold of dielectrics. (We note that material removal or material modification occurs at a timescale much larger than the light propagation scale so that it does not affect the intensity distribution [START_REF] Gamaly | The physics of ultra-short laser interaction with solids at nonrelativistic intensities[END_REF] )

We have numerically analyzed the quasi-propagation-invariant regime in terms of peak intensity and radial position of the main intensity ring for each propagation distance. We find that the actual propagation of finite energy pulses follows a family of propagation-invariant conical vortex waves, such as those shown in figure 2. This property is general for conical waves [START_REF] Polesana | Observation of conical waves in focusing, dispersive, and dissipative Kerr Media[END_REF] . This family is characterized by the same conical angle as the one of the input beam. In other terms, high intensity Bessel vortices reshape into a set of propagation-invariant solutions during the nonlinear propagation. Importantly, unlike solitons, this propagation-invariant solution sustains high losses occurring in the main ring. This is because a power flux from the external rings toward the main ring compensates the losses. For instance, energy loss at pulse energy of 5 µJ in Figure 3 exceeds 20% (see Methods section).

We observed both numerically and experimentally that the deviation from a quasipropagation-invariant regime occurs because of the appearance of nonlinear wave mixing, which generates novel spatial frequencies in the beam. As observed in the case of zeroth order Bessel beams [START_REF] Polesana | Filamentation in Kerr media from pulsed Bessel beams[END_REF][START_REF] Faccio | Nonlinear light-matter interaction with femtosecond highangle Bessel beams[END_REF] , the growth of modulation instability or fourwave mixing is inefficient at high conical angle, such as for the results shown in figure 3. The value of the conical angle 6.8

   was chosen to observe a clear quasi- stationary propagation regime at high pulse energy over a long distance. To observe the departure from propagation-invariant regime, we investigate the nonlinear propagation at lower cone angle, where nonlinear wave mixing growth is reached well before dense plasma generation limits the peak intensity.

Transitions from propagation-invariant to rotating and speckle-like regimes For a pulse energy of 5 µJ, the main ring splits azimuthally into several peaks that rotate around the optical axis. This seems to correspond to the regime observed by another group in water [START_REF] Shiffler | Self-focusing of femtosecond diffraction-resistant vortex beams in water[END_REF] . The sense of rotation is determined by the sign of the vortex charge. A second movie (Supplementary Material, movie 2), based on experimental results, shows the high degree of symmetry of the rotation during nonlinear propagation at 5 µJ for m=3 and m=-3. We note that in the linear case, no rotation was observed. Over 20 µJ, although one would expect multifilamentation, a novel regime appears, characterized by the fact that no continuous light channel is produced. Instead, multiple hot spots appear and disappear during the propagation with very limited rotation, resembling a speckle structure. Importantly, although very complex, this "speckle-like" regime is highly reproducible from shot to shot. Moreover, the same fluence distribution pattern was reproduced between identical experiments even when separated by several hours. This suggests that the speckle structure arises from inhomogeneity in the initial conditions rather than noise during propagation.

Spatial spectrum analysis

Spatial spectrum provides an efficient tool to understand the propagation dynamics in the low conical angle case. Experimentally, the Fourier spectrum was measured at the focal plane of a lens inserted in the imaging path of the previous setup (see methods).

In Fig. 5 For input pulse energies higher than 5 µJ, we observe the generation of spectral components with transverse wavevectors smaller than 0 r k in a complex pattern. These waves are generated by nonlinear wave mixing. They were removed when Kerr index was numerically set to zero and they disappeared in experiments when the pulse duration was temporally stretched. These spectral components arise from a combination of four-wave mixing and self-phase modulation. Simulations showed that the nonlinear propagation radially modulates the beam resulting in the generation of several concentric rings in the far field, as for the case of zeroth order nonlinear Bessel beams [START_REF] Gadonas | Self-action of Bessel beam in nonlinear medium[END_REF] . When intensity perturbations are added to the input beam, modulation instability splits these rings azimuthally. This effect is more and more visible as the input pulse energy increases.

We note that in the rotating regime, between ~ 5 and ~ 20 µJ, the number of hot spots does not necessarily correspond to the vortex order: it evolves with propagation. The actual number of peaks depends on the growth rate of each azimuthal mode [START_REF] Vuong | Collapse of optical vortices[END_REF] . At energies higher than 20 µJ, a number of azimuthally modulated rings are observed in the far field. They interfere in the near field and generate a highly complex, "speckle-like" pattern. The high degree of nonlinear cascading has already been observed to be deterministic [START_REF] Majus | Generation of periodic filament arrays by self-focusing of highly elliptical ultrashort pulsed laser beams[END_REF] . Besides, it makes the propagation extremely sensitive to nonlinear coefficients and a relatively small quantitative discrepancy is therefore observed here between simulations and experimental results, while the qualitative behavior of speckle-like propagation is still present in both.

Discussion. Applications to laser material processing

The propagation-invariant regime of Bessel Vortices opens numerous novel possibilities in the field of laser material processing. Indeed, we observed that this regime can be reached over a very wide range of parameters. Figure 6(a) shows a typical result of the longitudinal distribution of free-electron plasma in the propagation-invariant regime. Here, it corresponds to the case of the beam represented in figure 3 with a pulse energy of 5 µJ (i). The plasma distribution is longitudinally quasi-invariant over a range of 150 µm. The plasma cross section, shown in inset, peaks on the main intensity ring of the propagation-invariant Bessel vortex beam.

The tubular regime allows plasma generation on a very wide range of plasma densities and we note that density and geometrical parameters of the plasma tube can be controlled by six independent parameters (pulse energy, pulse duration, input Gaussian beam waist, cone angle, vortex order, laser central wavelength). for Corning 0211 and 0.5×10 -6 K -1 for fused silica).

The damage produced in glass is highly homogeneous with the shape of a cylinder extending from one side of the sample to the other. A pulse duration of 1 ps was chosen since the damage observed through optical transmission is much more visible than the one left by 100 fs, as already observed in previous work [START_REF] Bhuyan | Single shot high aspect ratio bulk nanostructuring of fused silica using chirped controlled ultrafast laser Bessel beams[END_REF] . The damage produced in fused silica under identical conditions as those for Corning 0211 was qualitatively similar but the index modification is smaller. Figure 6 (c) shows the cross section of a similar damage, obtained at pulse energy of 20 µJ. The picture is obtained under transmission microscopy and it is apparent that the central core of the damage is brighter than the other part of the transparent sample, thus indicating that the refractive index has increased by material compression and melting. Importantly, no light from the laser pulse has illuminated this central core. This core actually guides visible light and the near field image is shown in Figure 6(d). From the numerical aperture of the coupled light, we estimate the refractive index change to be ~8x10 -4 .

By using a similar approach as reference 23, we estimate the stopping distance of the shockwave to ~100 nm inside and outside the plasma tube [START_REF] Gamaly | The physics of ultra-short laser interaction with solids at nonrelativistic intensities[END_REF] . A model of cylindrical heat diffusion shows that laser deposited energy on the plasma tube is enough to overcome the melting temperature over the tube diameter. We note that the much high thermal expansion coefficient of 0211 glass compared to fused silica generates a higher thermal stress. This explains the higher material modification observed in Corning 0211 glass.

Our approach provides means to write waveguides with a single shot in media where waveguide writing by Gaussian beams usually generates negative index modification.

Hollow beams were previously used for this application in a regime with much smaller pulse energy, multishot illumination regime with continuous translation [START_REF] Long | Writing of stressed waveguides with tubular depressed cladding using femtosecond hollow beams[END_REF] . Here, an important novel step is that the high value of deposited energy in single shot allows the generation of material waves. Figure 6(e) shows schematically the propagation directions of inward and outward propagating cylindrical material waves generated from the tubular plasma. These waves are shockwaves, pressure waves and heat waves. The exact thermodynamical pathway yielding mechanical compression and material modification highly varies with radial energy distribution and material constants for heating, diffusion, thermal expansion and cooling after the plasma excitation. We found a variety of different damages depending on material and illumination parameters.

The tubular, propagation-invariant regime effectively de-couples light propagation from energy deposition within dielectric media. Our results are of primary importance because they open the way to control thermodynamics of materials under extreme conditions. Specifically, materials can be compressed and/or heated inside the hollow part of the beam, over arbitrary long distances, even if the material in the central core is of different composition as the dielectric where light propagates. In addition to very novel perspectives for laser material stress modification, high speed laser drilling and cutting of transparent materials, the multimegabar pressures generated by the explosion of dense plasmas from femtosecond pulses [START_REF] Juodkazis | Laser-induced microexplosion confined in the bulk of a sapphire crystal: evidence of multimegabar pressures[END_REF] can be applied to the synthesis of novel material phases in non-negligible amounts.

In conclusion, we have shown numerically and experimentally the existence of tubular propagation-invariant waves at high power in dielectric solids based on femtosecond Bessel beams with topological vortex singularity. Accurate numerical simulations and novel experimental approach of beam imaging allowed us to exhibit the transition from propagation-invariant to rotating and "speckle-like" regimes by a combination of four wave mixing and modulation instability. We anticipate that the impact of our results will also be important for filamentation in gases and specifically for long-range filamentation and microwave guiding by filaments in atmosphere [START_REF] Châteauneuf | Microwave guiding in air by a cylindrical filament array waveguide[END_REF] . We expect that propagation-invariant conical vortex waves and the subsequent novel degree of control on plasma geometry will have a dramatic impact in several new fields in physics such as plasma guiding with plasma tubes, high-aspect ratio laser material processing, femtosecond laser waveguide writing and material compression.

Methods

Experimental setup

Bessel vortex beam. The experimental setup is based on shaping the beam of a 130 fs Ti:Sa laser with a spatial light modulator (SLM). The SLM plane is imaged through a telescope with a first lens L 1 of focal length f 1 =1m and a microscope objective (MO 1 , x20) placed in a 4f configuration. In the image of the SLM, light amplitude distribution is demagnified by a factor 110. As described in reference 32, the SLM phase mask applies a combination of the target phase and a reference phase so that the diffraction orders are spatially separated in the focal plane of lens L 1 . Spatial filtering allows us to filter out all undesired diffraction orders. The Bessel vortex beam onsets from the focal plane of the microscope objective.

Sample imaging and positioning. The samples (1 mm thick) were mounted on a 3 axis motorized XYZ stage (bidirectional repeatability 0.14 µm) and a two axis tilt clear-aperture piezo-mount. To image both the sample and the femtosecond beam, we placed a x50, NA 0.8 microscope objective (MO 2 ) after the sample, on a separate translation stage (bidirectional repeatability 0.14 µm) to follow the sample z position along the optical axis. A CCD camera is placed in the focal plane of a f 3 =200 mm lens to image the sample with a magnification factor 55. Under white light illumination, the depth of field of sample imaging is less than 0.5 µm. The sample tilt compensation was then adjusted over the whole sample area (20x10 mm 2 ) yielding a planarity better than 1 µm over 20 mm (50 µrad).

Fluence distribution imaging.

For beam imaging, we use the same imaging system as for sample imaging. A high dynamical range (16 bits) CCD camera was electronically synchronized with the laser and imaged a single laser shot. Light intensity was attenuated by neutral densities placed before the camera and the camera gain was maintained to zero. Careful calibration of neutral densities, imaging setup transmission and camera response provides camera signal conversion in physical fluence (J/cm 2 ).

Imaging a plane inside the sample would be an indirect measurement of the fluence distribution in this plane since amplitude and phase can be affected by nonlinear propagation after this plane. Our imaging procedure makes use of the fact that the nonlinear propagation regime in fused silica has a finite extent, smaller than the sample thickness. In addition, we limit ourselves to pulse energies where we The spatial integration of fluence provides the pulse energy at each propagation distance, after correction of Fresnel losses. The ratio between the pulse energy before and after the high-intensity region provides the energy loss.

Far field imaging. Far field images were obtained by inserting a supplementary lens (f 2 =125 mm) between MO 2 and L 3 to perform an optical Fourier-transform.

Samples. We used high purity grade synthetic fused silica Lithosil ® Q1 from Schott.

The samples were 1 mm thick with a total thickness variation (TTV) less than 10 µm.

For laser processing and results shown on figure 6, we used Corning 0211 glass, with thickness 150 µm.

Modelling

Direct numerical simulations for beam propagation. The propagation model takes the canonical form of a unidirectional envelope propagation equation written in the spectral domain for the electric field envelope ˆ( , , )

x y E k k z : 22

2 2 2 0 0 2 0 ˆ( ) | | (1 / )| | (1 ) 2 2 K K z n t c E i k k E i n E E E E i z c                  (2)
where In addition, the effects of high intensity are of interest for potential applications of conical vortex waves to laser energy deposition over a tubular structure. We therefore considered the generation of an electron-hole plasma, with density  described by the rate equation:

2 0 | | (1 / ) / K nt r g E t K U                         (3)
where  denotes the cross section for inverse carrying a vortex charge m and a linear spatial phase equivalent to that induced by an axicon:

2 2 0 0 0 ( , , 0) exp / ( sin ) E r z E r r i k r m             (4)
Eigenvalue problem for propagation-invariant conical vortex waves. From previous works on nonlinear Bessel beams [START_REF] Porras | Nonlinear unbalanced Bessel beams: stationary conical waves supported by nonlinear losses[END_REF][START_REF] Polesana | Filamentation in Kerr media from pulsed Bessel beams[END_REF] , we know that the key ingredients in the propagation equation for finding propagation-invariant states in the form of nonlinear Bessel beams are the nonlinear losses and the focusing/defocusing nonlinearity. We note that, for the present study, UPPE-based model provides the same results as Nonlinear Schrödinger equation (NLS) -based numerical model, which justifies the use of the NLS equation here. We rewrite our model in the form of a NLS equation suitable for the analysis of stationary solutions, i.e., with a diffraction operator expanded in cylindrical coordinates.

2 2 2 2 2 2 2 0 1 1 [ (| | ) (| | )] 2 E i E if E g E E z k r r r r                     (5)
where

2 2 0 2 0 0 (| | ) ( / ) | | / 2 c f E k n n E       and 2 2 2 (| | ) | | /2 / 2 K K g E E       .
Without loss of generality, we looked for propagation-invariant solutions in the case 0   , i.e., we considered only one physical effect of each type, namely, the optical Kerr effect and multiphoton absorption losses. We seek propagation-invariant beams carrying angular momentum, in the form

  ( , , ) ( )exp ( ( ) ) E r z a r i r m z        , with 0  
resulting in the set of ordinary differential equations for ( ) a r and ( ) / q r d dr

  2 2 2 3 2 0 0 2 0 2 2 0 2 2 2 K K n a m a q a k a k a a r n r a q q q k a a r                (6)
where dots means differentiation with respect to r. Equations ( 5) and ( 6), together with boundary conditions: ( ) 0 a r  as 0 r  ; ( ) 0 a r  as r   ; ( ) 0 q r  as 0 r  , represent an eigenvalue problem for the eigenvector ( ) a r and eigenvalue  . In practice, a continuous spectrum of solutions exists for 0   , allowing us to equivalently fix the axial phase shift  from the cone angle  , through 2 0 ( / 2)sin k    . We then solve numerically equations ( 5) and ( 6) with boundary conditions compatible with the behavior of a Bessel beam:

( 1 ) 2 )

 12 Département d'Optique P. M. Duffieux, Institut FEMTO-ST, UMR 6174 CNRS Université de Franche-Comté, F-25030 Besançon cedex, France (Centre de Physique Théorique, CNRS, Ecole Polytechnique, F-91128 Palaiseau, France (3) Laboratoire Hubert Curien, UMR CNRS 5516, Université de Lyon, Université Jean Monnet, F-42000 Saint-Etienne, France francois.courvoisier@femto-st.fr An open challenge in the important field of femtosecond laser material processing is the controlled internal structuring of dielectric materials. Although the availability of high energy high repetition rate femtosecond lasers has led to many advances in this field, writing structures within transparent dielectrics at intensities exceeding 10 13 W/cm 2 has remained difficult as it is associated with significant nonlinear spatial distortion. This letter reports the existence of a new propagation regime for femtosecond pulses at high power that overcomes this challenge, associated with the generation of a hollow uniform and intense light tube that remains propagation invariant even at intensities associated with dense plasma formation. This regime is seeded from higher order nondiffracting Bessel beams, which carry an optical vortex charge. Numerical simulations are quantitatively confirmed by experiments where a novel experimental approach allows direct imaging of the 3D fluence distribution within transparent solids. We also analyze the transitions to other propagation regimes in near and far fields. We demonstrate how the generation of plasma in this tubular geometry can lead to applications in ultrafast laser material processing in terms of single shot index writing, and discuss how it opens important perspectives for material compression and filamentation guiding in atmosphere.

8 

 8 (b). The vortex Bessel beam is generated from the image of the SLM phase mask, placed at the focal plane of a high numerical aperture microscope objective. This allows us to generate the beam within the bulk of solid dielectrics. We have developed a novel beam procedure to image the fluence distribution in 3D for quantitative comparison with simulations. It is based on progressive precise scanning of the beam within the sample and single shot imaging of the sample exit side. The procedure is detailed in the methods section.Quasi-invariant propagation of finite energy beamsA first set of experiments and numerical simulations investigated propagation of high conical angle Bessel vortices, with 6.  in glass. Figure3shows the results for several input pulse energies. The subfigures in the first and second rows show the numerical and experimental longitudinal fluence distributions. The third row compares the transverse cross-sections for a propagation distance of 300 µm from the phase mask image in the sample. Aside from the high degree of agreement between simulation and experiment, it is apparent that both numerically and experimentally, no distortion of the central main ring appears even at high input energies: only the contrast between the main ring and the surrounding rings decreases. The transverse profile changes only smoothly along the propagation, which shows the propagation is quasi-propagation-invariant.

Figure 4 8 

 48 Figure 4 compares numerically and experimentally the fluence distribution with increasing pulse energies for a vortex of order m=3 and

Figures 5

 5 Figures 5 (c,d) show the spatial spectrum for an input pulse energy of 1 µJ. A slight increase of the ring radius is characteristic of the propagation-invariant regime at high intensity, where, in the near-field, lobe compression occurs (see Fig 2.a). As the input pulse energy further increases, self-phase modulation and four-wave mixing modify the spatial spectrum and generate one or more rings with radii close to the initial one.

Figure 6 ( 8 

 68 Figure 6(b) shows an image by optical transmission microscopy of the damage induced in 150 µm-thick glass slide (Corning 0211, see Methods section) by a single pulse with cone angle

  measured the propagation in air was linear. We image the fluence distribution only at the sample exit side: we reconstruct the 3D fluence distribution by translating the beam within the sample and always image the fluence at the sample exit side. The image of the SLM is first placed at the exit side of the sample, by translating the sample toward the laser source. The sample is illuminated by a single pulse and an image is recorded by our imaging system. Then, the sample is transversally shifted in plane by x  =50 µm in order to illuminate a fresh part of the sample and moved along the optical axis by a distance z  =2 µm for the high cone angle case and 5 µm for the low angle one. The imaging microscope objective MO 2 is translated axially by the same amount to image the sample exit side. The fluence distribution measured at the exit side then corresponds to the one at a propagation distance n z  within the sample, where n is the sample index of refraction.

  we record the maximum intensity of the conical vortex wave for all solutions that decay back to zero as a Bessel beam. In this way, we map the region of existence of stationary conical Bessel waves in the plane ( , )

Figure 1 Figure 1 .

 11 Figure 1

Figure 2 Figure 2 .

 22 Figure 2

Figure 3 Figure 3 .

 33 Figure 3

Figure 4 .

 4 Figure 4.

Figure 4 .

 4 Figure 4. Linear, stationary, rotating and "speckle-like" regimes. They are reached during the propagation for low conical angle (2.8°) and vortex charge m=3, when energy is progressively increased. (Top row) Numerical simulations starting from slightly distorted beam transverse profile. (Central row) Experimental results. (Bottom row) Comparison of the transverse numerical and experimental fluence profiles at a propagation distance z=0.82 mm, where each regime is fully developed.

Figure 5 Figure 5 .

 55 Figure 5

Acknowledgments

We acknowledge Luca Furfaro, Roland Salut, Amaury Mathis and Pierre-Ambroise Lacourt for their very valuable support. We acknowledge funding from Region Franche-Comte, French ANR, contract 2011-BS04-010-01 NANOFLAM. This work has been performed in cooperation with the Labex ACTION program (contract ANR-11-LABX-01-01).

Author contribution statement CX and VJ equally contributed to this work. VJ, CM and AC developed the semianalytical approach. VJ, CM, AC and TI developed numerical codes and performed the numerical simulations. AC and FC developed the original idea and conducted the research along with TI and JD. CX, RG, JD, FC developed the experimental setup. CX and IO acquired the experimental data that were processed by CX. The manuscript was jointly written by all co-authors. Supplementary Information: Movie 1: Experimental propagation in linear, propagation invariant, rotating and "speckle-like" regimes. This data is the same as the one shown in figure 4. Movie 2 : Comparison of the linear and rotating regimes for a vortex charge m=3 and m=-3, for pulse energy 5µJ.

The authors declare no competing financial interests.


Christophe Fluhr

Serge Grop

Yann Kersalé

Enrico Rubiola

Vincent Giordano

Benoît Dubois

Characterization of the individual frequency stability of Cryogenic Sapphire Oscillators at the 1e-16 level

Keywords: Ultra-stable Oscillators, frequency stability, phase noise

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

T he Cryogenic Sapphire Oscillator (CSO) based on a sapphire whispering gallery mode resonator cooled near 6 K is currently the microwave signal source presenting the highest short-term frequency stability for integration time τ = 1 . . . 10, 000 s. With a fractional frequency stability better than 1 × 10 -15 , the CSO allows the operation of the laser-cooled microwave atomic clocks at the quantum limit [START_REF] Santarelli | Quantum projection noise in an atomic fountain: A high stability cesium frequency standard[END_REF]. It provides the means to improve the resolution of the space vehicles ranging and Doppler tracking provided by Deep Space Networks as well as those of Very Long Baseline Interferometry (VLBI) Observatories [START_REF] Dick | Temperature compensated sapphire resonator for ultra-stable oscillator capability at temperatures above 77K[END_REF]- [START_REF] Rioja | The impact of frequency standards on coherence in VLBI at the highest frequencies[END_REF]. The CSO can also enhance the calibration capability of Metrological Institutes or help the qualification of high performances clocks or oscillators [START_REF] Santarelli | Quantum projection noise in an atomic fountain: A high stability cesium frequency standard[END_REF], [START_REF] Dolgovskiy | Ultra-stable microwave generation with a diode-pumped solid-state laser in the 1.5-µm range[END_REF], [START_REF] Takamizawa | Atomic fountain clock with very high frequency stability employing a pulse-tube-cryocooled sapphire oscillator[END_REF]. The CSO performances are only concurrenced at short integration times (τ = 0.1...10 s) by laser stabilized to a room temperature Fabry-Perot (FP) cavity made with low-expansion vitro-ceramic materials [START_REF] Millo | Ultrastable lasers based on vibration insensitive cavities[END_REF]- [START_REF] Häfner | 8×10 -17 fractional laser frequency instability with a long room-temperature cavity[END_REF]. More recently fractional frequency stabilities better than 1 × 10 -16 have been reported with laser stabilized on cryogenic Silicium FP cavity [START_REF] Kessler | A sub-40-mHz-linewidth laser based on a silicon single-crystal optical cavity[END_REF]. Conversly to the CSO, the FP-stabilized laser suffers from a large drift (typically 10 -12 /day, and requires a femtosecond (FS) laser to deliver the microwave output. In turn, the FS laser is expensive and cumbersome, it can hardly work longer than weeks without loosing internal frequency locking, and suffers from the low C. Fluhr, Y. Kersalé, E. Rubiola and V. Giordano are with the Institute FEMTO-ST (Franche-Comté Electronique, Mécanique, Thermique et Optique -Sciences et Technologies), Unité Mixte de Recherche 6174, CNRS (Centre National de la Recherche Scientifique), 25000 Besanc ¸on France, (e-mail: giordano@femto-st.fr) S. Grop is with Alemnis GmbH, Feuerwerkerstrasse 39, 3602 Thun Switzerland.

B. Dubois is with the center for technological development FEMTO Engineering, 32 avenue de l'Observatoire 25000 Besanc ¸on France.

Manuscript received XX, YY 2015; revised XX, YY 2015.

signal-to-noise ratio per mode inherent in the wide bandwidth.

Fractional frequency stability better than 1 × 10 -15 has already been demonstrated by beating two nearly identical CSOs [START_REF] Giordano | New generation of cryogenic sapphire microwave oscillator for space, metrology and scientific applications[END_REF], [START_REF] Hartnett | Ultra-low-phase-noise cryocooled microwave dielectric-sapphire-resonator oscillators[END_REF]. The fractional frequency stability of one unit is simply obtained by substracting 3 dB to the actual result, assuming the two oscillator noises are equivalent and uncorrelated. This assumption is not generally true especially when developing a new instrument at the state-of-the-art. In our lab, we recently achieved the implentation of three nearly identical CSOs in the frame of the Oscillator IMP project. This was the opportunity to test the three-cornered-hat method [START_REF] Gray | A method for estimating the frequency stability of an individual oscillator[END_REF] to extract the individual frequency stabilities. Very preliminary measurements based on these method have been presented in [START_REF] Fluhr | Characterization of a set of cryocooled sapphire oscillators at the 10 -16 level with the three-cornered hat method[END_REF]. Since one of the CSO has been improved and more data have been accumulated. The current results demonstrate the capability of the method, which has been simply implemented with commercial components and counters. Althought based on the same general configuration, the three CSOs lightly differ from each other. Different thermal configurations have been tested and only one resonator is completely optimized, the two others recently implemented still need some adjustments. The three-cornered-hat method gives us information about each CSO and thus will help us to optimize its functioning.

II. CSO DESCRIPTION

Our most advanced CSOs incorporate a 54 mm diameter and 30 mm height cylindrical sapphire resonator. It operates on the quasi-transverse magnetic whispering gallery mode W GH 15,0,0 near 10 GHz [START_REF] Grop | ELISA: A cryocooled 10 GHz oscillator with 10 -15 frequency stability[END_REF]. The Q factor can achieve 1 × 10 9 at the liquid-He temperature depending on the sapphire crystal quality and on the resonator adjustment and cleanning. In an autonomous cryocooled CSO, the sapphire resonator is placed into a cryostat and in thermal contact with the second stage of a pulse-tube (PT) cryocooler delivering typically 0.5 W of cooling power at 4 K (see Fig. 1). The gas flow in the cryocooler induces mechanical vibrations and a temperature modulation at about 1 Hz, which need to be filtered. In our cryostats the heat-links between the PT 2 nd stage and the flange supporting the 4 K thermal shield and the resonator is made with copper braids or foils. The same mechanical decoupling is implemented between the PT 1 st stage and the 50 K thermal shield. This simple arrangement is sufficient to limit the resonator displacement below 1 µm at the PT cycle frequency [START_REF] Grop | 10 GHz cryocooled sapphire oscillator with extremely low phase noise[END_REF]. The thermal filtering is obtained by combining the heat-link thermal resistance and the thermal mass of the 4 K flange and its ballast that could be added. Eventually, the resonator is stabilized at its turnover temperature T 0 , where its thermal sensitivity nulls at first order. T 0 depends on the residual paramagnetic impurities present in the sapphire crystal [START_REF] Giordano | Tests of sapphire crystals manufactured with different growth processes for ultra-stable microwave oscillators[END_REF], and thus is specific to each resonator. T 0 is typically between 5 and 8 K for a high-quality sapphire crystal.

The CSO is a Pound-Galani oscillator [START_REF] Galani | Analysis and design of a single-resonator GaAs FET oscillator with noise degeneration[END_REF]. In short, the resonator is used in transmission mode in a regular oscillator loop, and in reflection mode as the discriminator of the classical Pound servo [START_REF] Pound | Electronic frequency stabilization of microwave oscillators[END_REF]. The sustaining stage and the control electronics are placed at room temperature. The insertion loss through the cryostat is ≈ -30 dB. The sustaining stage is made up of commercial components. Two low noise microwave amplifiers provide a small signal gain of ca 54 dB. A Voltage Controlled Attenuator (VCA) allows the control of the power injected in the resonator. Two Voltage Controlled Phase Shifters (VCPS) are used for the Pound servo. A 70 kHz phase modulation is applied through the first one. The correction is applied through the second VCPS. A 80 MHz bandwidth filter and some isolators complement the circuit. The error signals needed for the Pound and the power servos are derived from the low frequency voltages generated by two tunnel diodes placed near the resonator input port. The Pound detector is directly connected to a lockin amplifier (Stanford Research Systems SR 810).

Three oscillators were assembled successively since 2012. They are identical in the principle but show however some differences (see table I).

The thermal ballast is a piece of stainless steel placed between the PT 2 nd stage and the resonator support. Associated with the thermal resistance of the mechanical link, it is equivalent to a first order filter with a time constant τ B .

D is the Pound servo frequency discriminator sensitivity in (V/Hz). It depends on P R : the power injected in the resonator. Formally, D is taken at the demodulator output (see Fig. 1). For slow frequency fluctuations ∆ν(t), the error signal at the demodulator output is ǫ(t) = D∆ν(t). D is experimentaly determined by applying an offset at the lockin amplifier output and measuring the resulting frequency shift.

CSO-1 is an optimized copy of our first demonstrator, i.e. ELISA developed for the European Space Agency (ESA) and implemented in the Deep Space Network station DSA-3 in Malargue Argentina [START_REF] Grop | ELISA: A cryocooled 10 GHz oscillator with 10 -15 frequency stability[END_REF]. The second one ULISS is a transportable unit, which has been travelled since 2012 in some european laboratories to be tested in real field applications [START_REF] Giordano | New generation of cryogenic sapphire microwave oscillator for space, metrology and scientific applications[END_REF]. The third CSO has a new designed cryostat and was put into operation in october 2014. No completely optimized, it incorporates a crystal manufactured by the Kyropoulos growth method instead of a HEMEX crystal.

III. RELATIVE FREQUENCY STABILITY MEASUREMENTS

The measurement set-up is schematized in figure 2. The CSO output signals are mixed to obtain the three beatnotes: ν 12 = 7 MHz, ν 13 = 0.9 MHz and ν 23 = 7.9 MHz. They are simply counted using a multi-channels K&K-FXE SCR counter [START_REF] Kramer | Multi-channel synchronous digital phase recorder[END_REF]. The three channels work in parallel thereby the data acquisitions are synchronous. All datas are processed by one second averaging time using a Λ windowing [START_REF] Rubiola | On the measurement of frequency and of its sample variance with high-resolution counters[END_REF]. The computed two-sample deviation σ Λ (τ ) differs from the true Allan deviation σ y (τ ). The correspondence between σ Λ (τ ) and σ y (τ ) can be found in [START_REF] Dawkins | Considerations on the measurement of the stability of oscillators with frequency counters[END_REF]. For white and flicker FM noise these two representations of the fractional frequency stability are almost identical: σ Λ (τ ) ≈ 1.14 × σ y (τ ).

A. Beatnotes and phase noise

Figure 3 shows σ Λ (τ ) computed from the three beatnotes and compared to the fractional stability of a typical highperformance hydrogen Maser (H.M.) The bold lines represent σ Λ (τ ) without any post data processing. The thin lines are computed from the data after a linear drift removing. The measured frequency stabilities are better than 1×10 -15 for τ ≤ 2, 000 s. If the oscillator noises are equivalent and uncorrelated, that means that the fractional frequency stability of one unit is better than 7 × 10 -16 , which is a conservative value as we will see in the next section. This number is coherent with the phase noise measurement between CSO-2 and CSO-3 shown in figure 4.

This result corresponds to the phase noise of one unit assuming equivalent and uncorrelated the two CSOs: 3 dB has been substracted from the measured spectrum. For Fourier frequency f ≤ 10 Hz, the phase noise spectrum can be approximated by a white frequency noise (f -2 slope) and a flicker frequency noise (f -3 slope):

S ϕ (f ) ≈ 1 f -2 + 0.2 f -3 × 10 -10 rad 2 Hz -1 (1)
This phase noise spectrum is equivalent in the time domain to a deviation σ Λ (τ ) such as [START_REF] Rubiola | Phase noise and frequency stability in oscillators[END_REF]:

σ Λ (τ ) = 8.2 × 10 -16 τ -1/2 + 6.0 × 10 -16 (2)
It should be noted that the CSO white PM noise S ϕ (f ) ≈ 10 -14 rad 2 Hz -1 is filtered by the Λ-averaging [START_REF] Dawkins | Considerations on the measurement of the stability of oscillators with frequency counters[END_REF]. This will not be the case with a traditional Π-counter. Thus the Allan deviation will appear limited at short integration times to 9 × 10 -15 τ -1 . due to the integration of the white PM noise over the counter input bandwidth of 10.9 MHz.

CSO-1 and CSO-3 present the best performances at short term, whereas a perturbation near τ = 10 s affects CSO-2. The long term behaviors also differ: CSO-1 and CSO-2 do not show a frequency drift but seems limited by a random walk process. Conversely CSO-3 is drifting with a rate of 2 -3 × 10 -14 /day.

B. Three Cornered Hat Method

The individual frequency stabilities have been computed using the three-cornered-hat method implemented in the Sta-ble32 software. The results are given the figure 5. At short term (τ < 50 s), all the individual σ Λ (τ ) improve with the integration time but do not follow the expected τ -1/2 slope. Indeed it is expected that the CSO shortterm frequency stability is ultimatly limited by the Pound discriminator white frequency noise. For CSO-1 and CSO-2 this white frequency noise is completely masked by another process leading to a hump in σ Λ (τ ). In section IV-A we show this perturbation can be attributed to residual resonator temperature fluctuations. At longer integration times, the frequency stabilities reach an apparent flicker floor expanding over approximatly two decades. Thus for the best oscillator, i.e. CSO-1, σ Λ (τ ) ≈ 1.5 × 10 -16 for 100 s ≤ τ ≤ 5, 000 s.

The humps appearing at τ ∼ 400 s, can be attributed to an oscillation in the air conditionning system of the lab. Thus they reveal the residual sensitivity to the room temperature of each CSO. Nevertheless the resulting CSO frequency variations are obviously correlated, which is in conflict with the requirement of the three-cornered-hat method. The calculated σ Λ (τ ) near 400 s cannot be considered as the actual frequency stability. Correlated noises indeed induces false results: the inversed hump in the CSO-1 stability is symptomatic of this situation. The same caution must be taken in the analysis of the long term fluctuations resulting also for a large part from the room temperature variations.

IV. SHORT TERM FREQUENCY STABILITY ANALYSIS

In the figure 6 the fractional frequency stability of CSO-1 is compared to the most advanced CSO performances previously published [START_REF] Giordano | New generation of cryogenic sapphire microwave oscillator for space, metrology and scientific applications[END_REF], [START_REF] Hartnett | Ultra-low-phase-noise cryocooled microwave dielectric-sapphire-resonator oscillators[END_REF], [START_REF] Locke | Invited Article: Design techniques and noise properties of ultrastable cryogenically cooled sapphire-dielectric resonator oscillators[END_REF]. These ones have been evaluated using the beatnote between two assumed equivalent CSOs. In [START_REF] Locke | Invited Article: Design techniques and noise properties of ultrastable cryogenically cooled sapphire-dielectric resonator oscillators[END_REF] Locke describes a liquid helium cooled CSO. The power injected in the resonator is higher than those used in our devices. Depending on the input resonator port coupling, he evaluates the optimimum injected power between 6 and 60 mW. We do not share this approach. In our CSOs the injected power is much lower, i.e P R ∼ 100 µW. It is set to the value for which the resonator power sensitivity is minimized [START_REF] Giordano | Influence of the electron spin resonance saturation on the power sensitivity of cryogenic sapphire resonators[END_REF], and thus is specific to each resonator. Locke measured the intrinsic noise of the Pound frequency discriminator. The resulting frequency stability limitation was evaluated to 2 × 10 -16 τ -1/2 , which is well below the observed value. It also demonstrated that the CSO frequency stability is affected by AM-index fluctuations of the interrogation signal. Such a sensitivity makes mandatory a supplementary control loop to suppress the spurious AM. In our CSOs we do not implement an AM suppression lock loop as we have never highlighted such a strong sensitivity.

In [START_REF] Hartnett | Ultra-low-phase-noise cryocooled microwave dielectric-sapphire-resonator oscillators[END_REF] Hartnett describes a cryocooled CSO where the injected power is set to 100 µW. He claims a fractional frequency stability limited by a white frequency noise process. However it is obvious from the figure 6 that σ Λ (τ ) does not follow the expected τ -1/2 slope for τ ≈ 1 s.

From the data presented in the figure 6, it is clear that there still exist some sources of fluctuation responsible for a frequency stability degradation with respect to the Pound servo intrinsic noise.

A. Noise in the resonator temperature control

The resonator is stabilized at its turnover temperature T 0 using a LakeShore 340 Temperature controller in the Proportional-Integral (PI) mode. The temperature sensor is a Cernox type CX-1050 with a sensitivity of approximately 3 kΩ/K. The controller drives a current I H in a 25 Ω heater. The LakeShore controller makes 10 readings per second and offers a variety of digital processing that can be done to the raw sensor data before applying the PI control equation. Thus the input information can be in sensor unit (Ω) or converted in temperature unit (K).

The measurement resolution deduced from the LakeShore 340 User's Manual is δR = 1 Ω equivalent to δT = 300 µK for the sensor we used. The same datasheet indicates a control stability equals to ±2 Ω (or equivalently ±600 µK) without any information about the considered bandwidth and the type of noise (white and/or flicker?) affecting the temperature control. Near T 0 the residual thermal sensitivity of the CSO frequency is:

1 ν 0 ∆ν ∆T = 1.9 × 10 -9 (T -T 0 ) (3)
A rough estimation of the noise floor imposed by the temperature controler can be done by assuming (T -T 0 ) ∼ δT and a rms temperature fluctuation ∆T ∼ 600 µK. The resulting fractional frequency stability is ∼ 3.4 × 10 -16 . This shows that below a fractional frequency stability of 1 × 10 -15 the noise of the temperature controller has to be considered. A temperature stability of some hundreds of microKelvin is a common performance for such a cryogenic configuration based on a commercial controller and a thermistor [START_REF] Hartnett | Cryogenic sapphire oscillator with exceptionally high long-term frequency stability[END_REF]- [START_REF] Hasegawa | High-precision temperature control and stabilization using a cryocooler[END_REF].

Better temperature stability of a few µK have been reported on systems using ac-resistance bridge and custom electronic design [START_REF] Luiten | Latest Results of the UWA Cryogenic Sapphire Oscillator[END_REF], [START_REF] Li | Cryogenic precision digital temperature control with peaked frequency response[END_REF]. Such a solution should be envisaged to improve the CSO short frequency stability and reach the limit imposed by the noise of the Pound frequency discriminator.

We do not know the details of the algorithms implemented in the controller. We checked that for small departures from the temperature set-point, the Model 340 behaves like a traditional PI controler. If ǫ(t) is the error signal expressed in Kelvin, the output current I H can be written as:

I H (t) = G P ǫ(t) + 1 τ I ǫ(t)dt (4)
G is an internal gain which depends on the controller configuration and on the sensor sensitivity. In our case we measured G = 22.5 mA/K and G = 1.7 mA/K when the controller works in sensor unit and temperature unit respectively. P is the dimensionless proportional gain. τ I is the integral time [START_REF] Åström | Advanced PID Control. ISA -The Instrumentation, Systems, and Automation Society[END_REF].

Both P and τ I can be adjusted from the controller front panel. The input data processing chosen, the setting of PI-controller gains is done using the autotuning procedure of the LakeShore controller. To start the autotuning, we adjusted the inital values of P and τ I to those determined by a manual Ziegler-Nichols procedure. Eventually, P and τ I can be slightly varied by checking the short term frequency stability and searching for the best result. Figure 7 shows σ Λ (τ ) calculated from the beatnote between CSO-1 and CSO-2 for two different configurations of the CSO-2 temperature controller: reading in sensor unit (Ω) or in temperature unit (K). The same P and τ I are used for both. The effect on the frequency stability is obvious: the hump maximum is shifted to longer integration times when the reading is made in Kelvin with all other parameters being kept constant. In a second time, we intentionaly operate CSO-1 far from its turnover temperature, i.e at T = T 0 + 100 mK. At this point the first order resonator temperature sensitivity is 1 ∆T ∆ν ν0 = 1.9×10 -10 /K. Thus the CSO frequency fluctuations will follow those of the resonator temperature. Figure 8 shows the resonator temperature deviation, i.e. σ T (τ ), for different temperature controller configurations.

The curve 4 is the computed temperature deviation, assuming a resonator temperature modulated around its mean value T with a period θ m = 20 s and an amplitude ∆T = 4 × 10 -5 K, which represents well the actual behavior of the configuration 2. The curve 3 corresponds to the best configuration we found and that has been used for the measurements presented in the previous section. The curve 3 is just above the resolution measurement limited by the frequency noise of CSO-2. The temperature modulation is still there leading to the hump in the CSO-1 fractional frequency stability curve, (see figure 5). In spite of all our efforts we did not manage to find a better tuning for CSO-1 and CSO-2. CSO-3 behaves better as no hump has been observed in its σ Λ (τ ) curve.

The cause of this temperature modulation as well as its impact on the CSO frequency stability, are still not well understood. A modulation amplitude as ∆T = 4×10 -5 K will lead to a frequency instability of 5 × 10 -16 if the resonator temperature is 6 mK above or bellow its turnover value T 0 . Such an error on the temperature setpoint seems not realistic as it is much higher than the controller measurement resolution of 300 µK. Moreover we tried to adjust the temperature setpoint by step of 1 mK around the expected T 0 without finding any better tuning. The temperature modulation could result from an unexpected time lag in the thermal system making the control loop unstable. CSO-1 and CSO-2 are operating for a long time. They have been subjected to several stops, and both have been transported by car in the frame of the ULISS project. Thermal anchorages and bold tightenings into the cryostat could have been degraded by the resulting mechanical perturbations. That is coherent with the fact that CSO-3 is immune to these thermal perturbation as it was recently assembled and has not been transported.

B. Noise in the Frequency discriminator

As explained in the previous section, CSO-3 is the only one not limited at short term by the resonator temperature fluctuations. Its fractional frequency stability is shown in figure 9 for several values of the Pound frequency discriminator gain D. The latter has been simply varied by changing the power P R injected into the resonator, all the other parameters being kept constant.

The measured frequency stabilities follow a τ -1/2 slope until τ ≈ 30s and then a flicker floor. Both vary with the inci- The frequency discriminator gain D being 1) 1.3×10 -4 V/Hz, 2) 3.5×10 -4 V/Hz, 3) 1 × 10 -3 V/Hz.

dent power indicating that the CSO-3 performance is limited by the intrinsic noise of the Pound frequency discriminator. In presence of a white voltage noise at the demodulator output, the CSO short term fractional frequency stability is [START_REF] Dawkins | Considerations on the measurement of the stability of oscillators with frequency counters[END_REF]:

σ Λ (τ ) = 2 3 e n D ν 0 τ -1/2 (5)
where e n (V/ √ Hz) is the demodulator output voltage spectral density. The computed values of e n are given in the table II. The constant value e n ≈ 8 nV/ √ Hz is the equivalent voltage noise of the Pound detector. The direct measurement of e n requires to duplicate the Pound detector and thus to place a second cryogenic diode receiving the signal reflected by the cavity. This has not been foreseen in our current CSO design. However this value is compatible with the expected noise contributions of the Lockin amplifier and of the diode detector itself.

Assuming an identical Pound detector noise for CSO-1 and CSO-2 leads to a short term frequency stability of 3 × 10 -16 τ -1/2 and 4.4 × 10 -16 τ -1/2 respectively. For this two CSOs the Pound discriminator noise is not the dominant process.

V. SUMMARY

We applied the three-cornered-hat method to measure the individual fractional short-term frequency stability of three Cryogenic Sapphire Oscillators. This method implemented with commercial instruments and softwares permits a comparison at the 10 -16 level. The method also reveals that there still exist technical sources of fluctuation responsible for a degradation of the oscillator frequency stability with respect to the Pound servo intrinsic noise. These perturbations could be minimized by a careful optimisation of the thermal system and of the resonator temperature stabilization. Despite these perturbations, all the tested CSOs present a short term frequency stability better than 7 × 10 -16 at 1 s and than 5 × 10 -16 between 30 s and 3, 000 s. The best CSO shows a frequency stability of 4.6 × 10 -16 at 1 s and a flicker floor below 2 × 10 -16 .

Fig. 1 .

 1 Fig. 1. Scheme of the Cryogenic Sapphire Oscillator showing the cryogenic part (cooled resonator) and the electronics placed at room temperature .

Fig. 2 .

 2 Fig.2. Measurement Set-Up. Each beatnote is low-pass filtered (10.9 MHz) and amplified (Minicircuits ZFL-1000+ Amplifier). The multichannel counter is referenced on the 10 MHz coming from an Hydrogen Maser (H.M.) .

Fig. 3 .

 3 Fig.3. Two-sample deviations σ Λ (τ ) calculated from the three beatnotes.

Fig. 4 .

 4 Fig. 4. One unit phase noise obtained by comparing CSO-2 and CSO-3. Measurement realized with the Symmetricom 5125A phase test set-up.

Fig. 5 .

 5 Fig.5. Individual fractional frequency stabilities obtained by applying the three-cornered-hat method.

2 Fig. 6 .

 26 Fig. 6. CSO-1 fractional frequency stability compared previous results: V. Giordano et al. [14]; • Locke et al. [26]; • J.G. Hartnett et al. [15] .

Fig. 7 .

 7 Fig. 7.σ Λ (τ ) calculated from the beatnote CSO-1 vs CSO-2, for two temperature controller configuration: 1) reading in sensor unit (Ω), 2) reading in temperature unit (K).

10 Fig. 8 .

 108 Fig.8. CSO-1 resonator temperature deviation σ T (τ ) for some temperature controller configurations, the proportional gain is kept constant P = 300: 1) reading in temperature unit (K) , τ I = 16 s; 2) reading in sensor unit (Ω), τ I = 16 s; 3) reading in sensor unit (Ω), τ I = 8 s, 4) Computed temperature deviation for a time varying temperature as: T (t) = T + ∆T sin(2π t θm ) with ∆T = 4 × 10 -5 K and θm = 20 s.

Fig. 9 .

 9 Fig.9. CSO-3 frequency stability extracted from three-cornered-hat method. The frequency discriminator gain D being 1) 1.3×10 -4 V/Hz, 2) 3.5×10 -4 V/Hz, 3) 1 × 10 -3 V/Hz.

TABLE II COMPUTED

 II VOLTAGE WHITE NOISE SPECTRAL DENSITY en

	Curve in Fig. 8	σ Λ (1s) ×10 -15	P R (µW)	D ×10 -4 (V/Hz)	en (nV/ √	Hz)
	1	5	12,5	1.3	8.0
	2	2	33	3.5	8.6
	3	0.7	87	9.3	8.0

ACKNOWLEDGMENTS

The work has been realized in the frame of the ANR project Equipex Oscillator IMP. The Oscillator IMP project funded in the frame of the french national Projets d'Investivement d'Avenir (PIA) targets at being a facility dedicated to the measurement of noise and short-term stability of oscillators and devices in the whole radio spectrum (from MHz to THz), including microwave photonics. The authors would like to thank the Council of the Région de Franche-Comté for its support to the Projets d'Investissements d'Avenir and the FEDER for funding one CSO.

Enrico Rubiola is full professor at the Université de Franche Comté and deputy director of the Department of Time and Frequency of the CNRS FEMTO-ST Institute, Besanc ¸on, France. Formerly, he was a full professor at the Université Henri Poincaré, Nancy, France, a guest scientist at the NASA JPL, a professor at the Università di Parma, Italy, and an assistant professor at the Politecnico di Torino, Italy. He graduated in electronic engineering at the Politecnico di Torino in 1983, received a Ph.D. in Metrology from the Italian Minister of University and Research, Roma (1989), and a Sc.D. degree from the Université de Franche Comté in 1999. Prof. Rubiola has worked on various topics of electronics and metrology, navigation systems, time/frequency comparisons, and frequency standards. His main fields of interest are precision electronics form dc to microwaves, and time and frequency metrology. This includes phase and amplitude noise, analog and digital frequency synthesis, high spectral purity oscillators, photonic systems, sophisticated instrumentation, spectral analysis, and noise. He has developed innovative instruments for AM/PM noise measurement with ultimate sensitivity, and a variety of dedicated signal-processing methods. Since 2012, Enrico is the PI of Oscillator IMP, a platform dedicated to the measurement of AM/PM noise and short-term frequency stability. A wealth of articles, reports and conference presentations are available on the Enrico's home page http://rubiola.org.

Vincent Giordano was born in


Loic Vichard
email: loic.vichard@utbm.fr

Simon Morando

Alexandre Ravey

Fabien Harel

Pascal Venet

Serge Pelissier

Daniel Hissel

Battery modeling using real driving cycle and Big-Bang Big-Crunch algorithm

The purpose of this study is to model battery packs integrated in a ten fuel cell hybrid electric vehicle fleet developed among the European project Mobypost dedicated to postal delivery applications. This project led to create and feed a big database as the vehicles were deeply monitored. Thanks to this database and Big-Bang Big-Crunch optimization algorithm this paper proposes a method to model battery using real driving cycle data in few minutes with a NRMSE less than 0.02.

I. INTRODUCTION

With the growing popularity of Electric Vehicles (EVs) and fuel cell hybrid electric vehicles (FCHEVs), batteries used in those applications are now a famous topic [START_REF] Blunier | Fuel cells, energy management using fuel cells and supercapacitors[END_REF]. However it is frequently pointed out that the batteries are still too expensive and a very critical part among the entire electric power-train. Indeed their lifetime and their autonomy are still a limitation to a worldwide deployment as it cannot compete internal combustion engine (ICE) yet. Consequently in order to improve the energy management of the batteries and thus their autonomy and lifetime, a good understanding of the dynamic behavior is required. This is why mathematical models are needed in order to represents and predict this dynamic behavior. Generally battery models are determined from battery characterizations by using constant current or simplified current profile. Models can also be determined by using impedance spectroscopy. These processes require to set up experiments needing equipment, cost and time. The purpose of this paper is to present a fast way to model a battery cell through a Thevenin Electrical Equivalent Circuit (EEC) by using an optimization algorithm and field use experimental data. The particularity of this study is the use of a rich database created from an ten fuel cell hybrid electric vehicle fleet.

In a first part the Mobypost European project and the vehicles will be presented. In a second part the data used to model the batteries will be enumerated. In a third part Thevenin representation and mathematical equations will be presented. Then, in a fourth part, the optimization algorithm used to identified the model parameters will be explained. Finally in fifth part the results will be illustrated and discussed to go on conclusion and perspectives.

II. MOBYPOST PROJECT

A. Concept

The european project Mobypost aimed at developing and instrumenting a fleet of ten hybrid hydrogen vehicles powered by fuel cell for postal delivery applications (Fig. 1) [START_REF] Ravey | Energy management of fuel cell electric vehicle with hydrid tanks[END_REF], [START_REF] Vichard | Battery aging study using field use data[END_REF]. The design of such vehicles, lightweigth hybrid electric including fuel cell and batteries, includes a powertrain design and components sizing taking into account the type of application : postal delivering services. The powertrain was designed using serie architecture by operating the fuel cell as a range extender [START_REF] Higel | Part 1: Mobypost vehicles powertrain modeling, simulation and sizing[END_REF], [START_REF] Faivre | Part 2-mobypost vehicles powertrain design and experimental validation[END_REF].

The battery pack is composed of four modules (U24-12XP) connected in series and manufactured by Valence. Those modules, present a nominal voltage of 12.8V and a nominal capacity of 110Ah measured with a current C/5 at 23°C. Each module is also composed of four parallel 3.2V LiFeMgPO4 cell groups in series.

B. fleet monitoring and used data

This project led to create a big database by monitoring ten postal delivery vehicles during almost three years. Physical values from sensors are measured and saved on a server every second [START_REF] Vichard | Battery aging study using field use data[END_REF]. Now an access to data of more than 1500 real driving cycles can be used to model the battery cells and observe the aging. In this study the data used are the battery current (𝐼𝑏𝑎𝑡𝑡), the state of charge (𝑆𝑂𝐶), the battery cell voltage (𝑉 𝑐𝑒𝑙𝑙) and the battery temperature (𝑇 ).

The following figures Fig. 2

A. Thevenin model

Thevenin electrical equivalent circuit can be used to represent the temporal voltage response of a battery. Among all the equivalent circuits which can be found in the litterature, a most common structure is illustrated in the Fig. 6. The battery cell voltage response being composed of transient periods of many different duration, this is necessary to add several R//C circuits. The R//C circuits will get different time constants (𝜏 ) in order to make the model capable to reproduce larger interval of transient periods duration. In this structure, R//C series circuits are used to represent impedances linked to the charge transfers effects, double layer and diffusion effects. In the literature, models can be found integrating only one R//C circuit as in [START_REF] Tian | An experimental study and nonlinear modeling of discharge i-v behavior of valve-regulated leadacid batteries[END_REF]- [START_REF] Dubarry | Development of a universal modeling tool for rechargeable lithium batteries[END_REF], with two R//C circuits [START_REF] Chen | Accurate electrical battery model capable of predicting runtime and iv performance[END_REF]- [START_REF] Hu | A technique for dynamic battery model identification in automotive applications using linear parameter varying structures[END_REF], or again with three R//C circuits [START_REF] Kroeze | Electrical battery model for use in dynamic electric vehicle simulations[END_REF] [START_REF] Hu | Electro-thermal battery model identification for automotive applications[END_REF]. In this study the Thevenin structure is used with three R//C circuits in order to get more accuracy. Indeed these three R//C circuits are used to represent the charge/discharge phases, the diffusion phases and the relaxation phases.

B. Model equations

Using Kirchhoff's circuit laws from the Fig. 6 , the equation ( 1) can be used to determine Vcell, the estimated voltage response of the battery.

𝑒𝑠𝑡𝑖𝑚𝑎𝑡𝑒𝑑 𝑉 𝑐𝑒𝑙𝑙 = 𝑂𝐶𝑉 -𝑈 0 -𝑈 1 -... -𝑈𝑛 (1)
The next step is to resolve the differential equation ( 2) due to the capacitor. After having resolved (2), the equation (3) was deduced as solution in order to express Un. It has to be notified, since the acquisition frequency is 1𝐻𝑧, t is then discretized with a range of 1 second to solve the equation.

U𝑛 = 1 𝐶𝑛 .𝐼 - 1 𝑅𝑛.𝐶𝑛 .𝑈 𝑛 (2) 𝑈𝑛(𝑡) = 𝑈𝑛(𝑡 -1).𝑒 -1/(𝑅𝑛.𝐶𝑛) + 𝑅𝑛.(1 -𝑒 -1/(𝑅𝑛.𝐶𝑛) ).𝐼(𝑡) (3)
In order to improve both accuracy and reality meaning, some dependencies have to be taken into account. Open Circuit Voltage (OCV) will be dependent on the SOC and on the Temperature as in the function (4). This OCV expression has been deduced during this study after a long work to find out how to reproduce OCV evolution versus the SOC and the temperature. Internal resistances will be dependent only on the temperature as in the function [START_REF] Faivre | Part 2-mobypost vehicles powertrain design and experimental validation[END_REF].

a, b, c, d, e, f, g, h and i being constants to be identified

𝑅 𝑛 (𝑡) = 𝑎 𝑛 .𝑒𝑥𝑝(𝑏 𝑛 .𝑇 ) (5)
a 𝑛 and b 𝑛 being constants to be identified From these equations, an objective function can be deduced. Then an optimization algorithm is run using the objective function : Normalized Root Mean Square Error (NRMSE) [START_REF] Dubarry | Development of a universal modeling tool for rechargeable lithium batteries[END_REF] in order to determine the parameter values of the model. This NRMSE is obtained from the following equations : Error (6) and Root Mean Square Error (RMSE) [START_REF] Liaw | Modeling of lithium ion cellsa simple equivalent-circuit model approach[END_REF].

𝐸𝑟𝑟𝑜𝑟 = 𝑟𝑒𝑎𝑙 𝑉 𝑐𝑒𝑙𝑙 -𝑒𝑠𝑡𝑖𝑚𝑎𝑡𝑒𝑑 𝑉 𝑐𝑒𝑙𝑙 ( 6
)
𝑅𝑀 𝑆𝐸 = ⎷ 1 𝑛 . 𝑛 ∑ 𝑖=1 𝐸𝑟𝑟𝑜𝑟 2 (7) 𝑁𝑅𝑀𝑆𝐸 = 𝑅𝑀 𝑆𝐸 𝑚𝑎𝑥(𝑟𝑒𝑎𝑙 𝑉 𝑐𝑒𝑙𝑙 ) -𝑚𝑖𝑛(𝑟𝑒𝑎𝑙 𝑉 𝑐𝑒𝑙𝑙 ) (8)
IV. OPTIMIZATION ALGORITHM

A. Optimization concept

Optimization problem can be summarized as finding the set of input parameters (x 𝑛 ) minimizing a mathematical function f(x 𝑛 ). Nature is a main source for proposing optimization methods such as ant colony system (ACO), genetic algorithms (GA) or again big-bang big-crunch (BB-BC). In this paper, the objective is to identified the set of inputs parameters (parameters of the Thevenin structure), minimizing the NRMSE which means the 𝑒𝑠𝑡𝑖𝑚𝑎𝑡𝑒𝑑 𝑉 𝑐𝑒𝑙𝑙 is as close as possible than the 𝑟𝑒𝑎𝑙 𝑉 𝑐𝑒𝑙𝑙 .

B. Big-Bang Big-Crunch algorithm

In this study the BB-BC algorithm has been used. This optimization algorithm relies on one of the theories of the universe evolution [17] [18]. This theory, based on the general relativity of Einstein, explains that the expansion phase of the universe due to the big bang will end to be replaced by the big crunch phase which consists of the universe contraction. This algorithm is used in this study as it shows good results and efficiency in [START_REF] Morando | Reservoir computing optimisation for pem fuel cell fault diagnostic[END_REF]. It refers to the creation of the initial population randomly. This phase consists in spreading the candidates solutions all over the search space in an uniform manner.

The BB-BC optimization algorithm is based on the sequence of the following steps :

Step 1: Create an initial population of N 𝑝𝑜𝑝 candidates randomly, respecting the limits of the search space.

Step 2: Evaluate the objective function value of all candidate solutions from the initial population.

Step 3: Find the center of mass referring to the best candidate solution. (candidate giving the lower error).

Step 4: Calculate new candidates around the center of mass by adding or subtracting a normal random number which its value decreases as the iterations elapse (reduction of the search space). Also an exploration probability is integrated in order to get out of an eventual local minimum.

Step 5: Return to step 2 until stopping criteria is met.

The reduction of the search space is represented by a function inversely proportional to the number of generation. The equation ( 9) illustrates the research space division according to the number of generations.

𝑅𝑆 𝑑𝑖𝑣 (𝑔𝑒𝑛) = 1 𝑔𝑒𝑛 (9)
This function drastically reduces the search space for the first generations and this can lead to progress to a local minimum instead of the desired global minimum. In order to avoid this, another way to reduce the search space, proposed in [START_REF] Morando | Reservoir computing optimisation for pem fuel cell fault diagnostic[END_REF], has been used in this study. Indeed this way divides more slowly the search space during first generations and then accelerate to find the results more efficiently. The proposed algorithm replaces the linear function by a translated exponential function.

It has to be specified that this algorithm is stochastic and to be sure to achieve the global minimum, an infinite number of iterations is required. In this study what is desired is a set of inputs parameters which will be able to reproduce the voltage for the entire temperature range with a good accuracy. Also as a short computing duration, is desired, a big population is avoided, and not a lot of generations are proceeded so here finding the global minimum is not the target. Tables I and II display the three average time constant values and their capacitance values for the driving cycle which was used for parameters identification. These tables also show the values of the constants in equations ( 4) and ( 5). The time constants present three different values, 2 seconds to reproduce charge transfer and double layer effects, 66 seconds and 500 secondes to reproduce both diffusion phenomena and relaxation phases. Regarding relaxation phenomena, to be able to better reproduce it at best, others R//C circuits should be added with higher time constants up to several hours [START_REF] Li | Fast characterization method for modeling battery relaxation voltage[END_REF]. The Fig. 8 illustrates the evolution of the internal resistance 𝑅𝑖𝑛 (which is the sum of the four resistances) versus the temperature. These results can lead to the following conclusion. The model is able to reproduce the voltage response of one driving cycle by taking into account the dependencies on the temperature and on the state of charge with a good accuracy after having identified the parameters. Here results are illustrated with the driving cycle which has been used to identify the parameters, so the method is not validated and it must be verified that the parameters are able to reproduce all the driving cycles of this battery.

V. RESULTS

VI. VALIDATION

To validate this modeling and verify if the identified parameters are coherent, another driving cycle has to be reproduced to observe if the error remains low. To do that, a driving cycle where the temperature is between 5 ∘ 𝐶 and 9 ∘ 𝐶 was chosen. In this study a battery model for one driving cycle taking into account dependency on temperature and state of charge has been deduced by using a Thevenin electrical equivalent circuit and data directly extracted from real driving cycle through an optimization algorithm. This method presents a good accuracy and efficiency as the computing time is few minutes and the NRMSE is low. Future works should aim at optimizing again the dependency on the temperature and on the SOC.

VIII. ACKNOWLEDGMENT

This research work was funded by the region of Franche-Comté, France.

Fig. 1 :

 1 Fig. 1: 3 vehicles of the Mobypost fleet

Fig. 3 Fig. 4 Fig

 34 . 5 illustrate the data mentioned above recorded during one postal delivery. This driving cycle was realized using the vehicle only in pure electric mod (fuel cell shutdown).

Fig. 2 :Fig. 3 :

 23 Fig. 2: Battery current during the driving cycle

Fig. 4 :Fig. 5 :

 45 Fig. 4: State of charge evolution along the driving cycle

Fig. 6 :

 6 Fig. 6: Battery RC model

  𝑂𝐶𝑉 (𝑡) = 𝑎 + 𝑏.(25 -𝑇 (𝑡)) 𝑆𝑂𝐶(𝑡) + 𝑐 𝑆𝑂𝐶(𝑡) + 𝑑.𝑆𝑂𝐶(𝑡) + 𝑒.𝑙𝑛(𝑆𝑂𝐶(𝑡)) + 𝑓.𝑙𝑛(1.001 -𝑆𝑂𝐶(𝑡)) + 𝑔.𝑙𝑛(1.01 -𝑆𝑂𝐶(𝑡)) + ℎ.𝑒𝑥𝑝(𝑖.𝑇 (𝑡))

Figures 7

 7 Figures 7 and 10 compare the estimated battery voltage from the model regarding the real measured cell voltage from the database. The error that presents the model is 𝑁𝑅𝑀𝑆𝐸 = 0.0185 proving the effectiveness of the proposed method.

Fig. 7 :

 7 Fig. 7: Cell voltages comparison with driving cycle used for parameters identification

  Figures 9 and 11 can validate the identified parameters as the model is able to adjust the voltage in function of the temperature range.

Fig. 9 :

 9 Fig. 9: Cell voltages comparison with driving cycle used for validation

Fig. 10 :Fig. 11 :

 1011 Fig. 10: Cell voltages comparison with entire driving cycle used for parameters identification

TABLE I :

 I Model parameters values and time constants for 20.5°C<T<24.5°C

	Components	Values
	𝜏 1	66 𝑠
	𝜏 2	2 𝑠
	𝜏 3	500 𝑠
	𝐶1	72000𝐹
	𝐶2	2700𝐹
	𝐶3	500000𝐹
	𝑎	3.22 𝑉
	𝑏	-0.001 𝑉 / ∘ 𝐶
	𝑐	-0.0216 𝑉
	𝑑	0.182 𝑉
	𝑒	-0.0886 𝑉
	𝑓	-0.0632 𝑉
	𝑔	0.0778 𝑉
	ℎ	-0.001 𝑉
	𝑖	-0.5 ∘ 𝐶 -1

TABLE II :

 II Model resistance coefficient values

		Resistance	𝑎𝑛	𝑏𝑛
		𝑅0	0.0003 Ω	-0.0753 ∘ 𝐶 -1
		𝑅1	0.0010 Ω	-0.0052 ∘ 𝐶 -1
		𝑅2	0.0013 Ω	-0.0301 ∘ 𝐶 -1
		𝑅3	0.0010 Ω	-0.0025 ∘ 𝐶 -1
		-3
		2.64
	Internal Resistance Rin (Ohms)	2.52 2.54 2.56 2.58 2.6 2.62
		2.5
		2.48
		20 20.5 21 21.5 22 22.5 23 23.5 24 24.5
			T (°C)

Rin

Fig.

8

: Internal resistance evolution versus temperature