Thomas Laupretre

Lucas Groult

Bachir Achi

Michael Petersen

Yann Kersalé

AND Marion Delehaye

Clément Lacroute
email: clement.lacroute@femto-st.fr

Thomas Lauprêtre

Clément Lacroûte

Absolute frequency measurements of the 1S0 -&gt 1P1 transition in ytterbium

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The ytterbium atom has several advantages for modern applications in atomic physics. Cold ytterbium experiments have therefore spread worldwide with applications in metrology, quantum information and quantum gases. Neutral ytterbium offers a narrow optical clock transition used for state-of-the-art frequency metrology [START_REF] Nemitz | Frequency ratio of Yb and Sr clocks with 5 × 10 -17 uncertainty at 150 seconds averaging time[END_REF][START_REF] Schioppo | Ultrastable optical clock with two cold-atom ensembles[END_REF] and can be cooled down to ultra-low temperatures to study the physical properties of ultra-cold quantum gases [START_REF] Aguilera | Non-linear relaxation of interacting bosons coherently driven on a narrow optical transition[END_REF]. In addition, the Yb + ion offers two optical clock transitions at 436 nm and 467 nm [START_REF] Godun | Frequency ratio of two optical clock transitions in 171 Yb + and constraints on the time variation of fundamental constants[END_REF][START_REF] Huntemann | Single-ion atomic clock with 3 × 10 -18 systematic uncertainty[END_REF], and its microwave clock transition at 12.6 GHz can be used for q-bit manipulation in quantum information processing experiments [START_REF] Randall | Generation of high-fidelity quantum control methods for multilevel systems[END_REF].

The 1 S 0 → 1 P 1 transition is the first cooling transition of neutral ytterbium and the first isotopic-selective stage in 2-photon ionization of ytterbium. However, there is a large dispersion in the measurements of the 1 S 0 → 1 P 1 transition frequency [START_REF] Meggers | The first spectrum of ytterbium (Yb I)[END_REF][START_REF] Das | Absolute frequency measurements in Yb with 0.08 ppb uncertainty: Isotope shifts and hyperfine structure in the 399-nm 1 S 0 → 1 P 1 line[END_REF][START_REF] Nizamani | Doppler-free Yb spectroscopy with the fluorescence spot technique[END_REF][START_REF] Enomoto | Comparison of resonance frequencies of major atomic lines in 398-423 nm[END_REF][START_REF] Kleinert | Measurement of the Yb I 1 S 0 -1 P 1 transition frequency at 399 nm using an optical frequency comb[END_REF]. The first measurement was performed in 1977 by acquiring spectrograms [START_REF] Meggers | The first spectrum of ytterbium (Yb I)[END_REF] and listed the energy of the transition as 25 068.227 cm -1 (751.526 54 THz), which was used as the reference for the NIST database [START_REF] Kramida | NIST Atomic Spectra Database (Ver. 5.6)[END_REF]. More recent measurements started in 2005, when Das et al. used a ring-cavity resonator stabilized on a Rb transition to measure the absolute frequency in 174 Yb with high accuracy [START_REF] Das | Absolute frequency measurements in Yb with 0.08 ppb uncertainty: Isotope shifts and hyperfine structure in the 399-nm 1 S 0 → 1 P 1 line[END_REF]. However, the value they report, 751.525 987 761(60) THz, differs from the one measured in [START_REF] Meggers | The first spectrum of ytterbium (Yb I)[END_REF] by more than 550 MHz. Nizamani et al. measured the transition frequency in 2010 [START_REF] Nizamani | Doppler-free Yb spectroscopy with the fluorescence spot technique[END_REF] using a wavelength-meter specified with a 60 MHz accuracy. Their measurement, yielding the value 751.526 65(6) THz, is closer to [START_REF] Meggers | The first spectrum of ytterbium (Yb I)[END_REF] (110 MHz difference) but not compatible with [START_REF] Das | Absolute frequency measurements in Yb with 0.08 ppb uncertainty: Isotope shifts and hyperfine structure in the 399-nm 1 S 0 → 1 P 1 line[END_REF]. In 2016, Enomoto et al. [START_REF] Enomoto | Comparison of resonance frequencies of major atomic lines in 398-423 nm[END_REF] measured the F = 1/2 → F = 3/2 hyperfine transition frequency in 171 Yb to be 751.527 48 [START_REF] Enomoto | Comparison of resonance frequencies of major atomic lines in 398-423 nm[END_REF] THz by accurately calibrating resonance frequencies of an ultra-stable Fabry-Perot cavity with lasers referenced to transitions in rubidium and calcium. They did not measure the transition frequency in the 174 Yb isotope, but their result is in good agreement with [START_REF] Meggers | The first spectrum of ytterbium (Yb I)[END_REF] and with [START_REF] Nizamani | Doppler-free Yb spectroscopy with the fluorescence spot technique[END_REF], given their uncertainty and the existing isotopic shift measurements from literature [START_REF] Das | Absolute frequency measurements in Yb with 0.08 ppb uncertainty: Isotope shifts and hyperfine structure in the 399-nm 1 S 0 → 1 P 1 line[END_REF][START_REF] Nizamani | Doppler-free Yb spectroscopy with the fluorescence spot technique[END_REF][START_REF] Kleinert | Measurement of the Yb I 1 S 0 -1 P 1 transition frequency at 399 nm using an optical frequency comb[END_REF][START_REF] Deilamian | Isotope shifts and hyperfine splittings of the 398.8-nm Yb I line[END_REF][START_REF] Loftus | Optical double-resonance cooled-atom spectroscopy[END_REF][START_REF] Banerjee | Precise measurements of UV atomic lines: Hyperfine structure and isotope shifts in the 398.8 nm line of Yb[END_REF][START_REF] Wang | A novel method to measure the isotope shifts and hyperfine splittings of all ytterbium isotopes for a 399-nm transition[END_REF]. Using an optical frequency comb referenced to the D 2 line of Rb, Kleinert et al. [START_REF] Kleinert | Measurement of the Yb I 1 S 0 -1 P 1 transition frequency at 399 nm using an optical frequency comb[END_REF] made a very accurate measurement of the transition frequency in 2016, measuring 751.526 533 49(33) THz, disagreeing with the values from [START_REF] Das | Absolute frequency measurements in Yb with 0.08 ppb uncertainty: Isotope shifts and hyperfine structure in the 399-nm 1 S 0 → 1 P 1 line[END_REF] and [START_REF] Nizamani | Doppler-free Yb spectroscopy with the fluorescence spot technique[END_REF].

Given the non-negligible dispersion of the measurements published so far, this work aims to establish and report on additional independent measurements of the absolute frequency of the 1 S 0 → 1 P 1 transition in Yb. Two independent measurements of this transition for different isotopes of ytterbium are obtained by recording the fluorescence spectrum of an atomic beam. Care was taken in the experimental setup to mitigate Doppler shifts and power broadening effects, reducing systematic effects below 5 MHz. The weighted mean of our measurements results in the value 751 526 537 ± 27 MHz for the 174 Yb isotope transition frequency.

Experimental setup

Fig. 1. Sketch of the experimental setup. In an ultra-high vacuum chamber, a 399 nm laser beam crosses an atomic beam of neutral ytterbium atoms to excite the 1 S 0 → 1 P 1 transition. The resulting fluorescence is detected by a low-noise photodiode. An angle θ can be introduced with respect to the perfect orthogonal configuration and a mirror allows for back reflection of the laser beam for alignment purposes (see text for details).

The experimental set-up is sketched out in Fig. 1. A commercial dispenser of ytterbium with natural isotopic abundance [17] is driven at 5.5 A and thus heated close to 600 • C in order to produce an atomic beam inside an ultra-high vacuum chamber (background pressure < 10 -10 mbars). The 1 S 0 → 1 P 1 transition in ytterbium is excited by a 399 nm laser beam crossing the atomic beam nearly orthogonally. The laser light, produced by a commercial external-cavity diode laser (ECDL) [18] with a linewidth of a few MHz, is brought to the experiment chamber via a single-mode fiber. The fiber output coupler is mounted on a rotatable breadboard in order to adjust the angle θ between the probe laser and the direction perpendicular to the atomic beam. The average intensity of the probe beam is ∼ 8 mW.cm -2 corresponding to a saturation parameter of 0.14. The laser frequency is slowly scanned over a ∼ 3 GHz range using the piezoelectric actuator driving its grating (PZT) during approximately 100 s. High-apperture optics collect photons scattered by atoms excited to the upper state and decaying back to the ground state. The fluorescence signal is detected with a photodiode [19]. While the fluorescence spectrum of the transition is recorded, the frequency of the laser is sampled every ∼ 25 ms on two independent commercial wavemeters [20]. The first wavemeter is a WS/7 Super Precision specified by the manufacturer with an absolute accuracy of ±60 MHz when calibrated with its built-in neon lamp. Its frequency drift in absence of thermal isolation has been measured to be below 30 MHz/day [START_REF] Saleh | Frequency stability of a wavelength meter and applications to laser frequency stabilization[END_REF]. The second wavemeter is a WS8-2 calibrated to the F = 3 → F = 4 hyperfine transition of the D 1 line of caesium at 894 nm, which results in an absolute accuracy of ±30 MHz at 399 nm as specified by the manufacturer.

A mirror placed at the output of the chamber can be used to reflect the laser beam back onto itself. In that case, the two counter-propagating beams have perfectly opposite Doppler shifts with respect to the atomic beam for non-zero angles θ, allowing to tune θ to zero with mrad precision (see section 4).

Results

The fluorescence spectrum was recorded at θ = 0 • as a function of the frequency for each wavemeter. A Lorentzian fit is applied to the data to extract the frequency of the 1 S 0 → 1 P 1 transition for each of the ytterbium isotopes. Fig. 2 shows the spectrum obtained with the WS8-2. The common estimated full width at half-maximum (FWHM) of 37 MHz is larger than the 28 MHz natural linewidth. This minor broadening is mainly attributed to laser linewidth (∼ 4 MHz), power broadening (∼ 2 MHz) and Doppler broadening due to the transverse velocity distribution in the atomic beam. [START_REF] Kleinert | Measurement of the Yb I 1 S 0 -1 P 1 transition frequency at 399 nm using an optical frequency comb[END_REF]. In [START_REF] Enomoto | Comparison of resonance frequencies of major atomic lines in 398-423 nm[END_REF], the authors do not provide a frequency value for the 174 Yb isotope. For comparison purposes, we use the value of 751 526.65 ± 0.11 GHz calculated from their measurement for the 171 Yb isotope translated by an isotopic shift of 832.6 MHz, the weighted mean [START_REF]Using the inverse of the squared uncertainties as weights[END_REF] of the shifts measurements from [START_REF] Das | Absolute frequency measurements in Yb with 0.08 ppb uncertainty: Isotope shifts and hyperfine structure in the 399-nm 1 S 0 → 1 P 1 line[END_REF][START_REF] Nizamani | Doppler-free Yb spectroscopy with the fluorescence spot technique[END_REF][START_REF] Kleinert | Measurement of the Yb I 1 S 0 -1 P 1 transition frequency at 399 nm using an optical frequency comb[END_REF][START_REF] Deilamian | Isotope shifts and hyperfine splittings of the 398.8-nm Yb I line[END_REF][START_REF] Loftus | Optical double-resonance cooled-atom spectroscopy[END_REF][START_REF] Banerjee | Precise measurements of UV atomic lines: Hyperfine structure and isotope shifts in the 398.8 nm line of Yb[END_REF][START_REF] Wang | A novel method to measure the isotope shifts and hyperfine splittings of all ytterbium isotopes for a 399-nm transition[END_REF]]. An uncertainty of 10 MHz was added to account for the uncertainty on the isotopic shifts measurements. The reported values are spread out over 600 MHz. When restraining the analysis to values from [START_REF] Meggers | The first spectrum of ytterbium (Yb I)[END_REF][START_REF] Nizamani | Doppler-free Yb spectroscopy with the fluorescence spot technique[END_REF][START_REF] Enomoto | Comparison of resonance frequencies of major atomic lines in 398-423 nm[END_REF][START_REF] Kleinert | Measurement of the Yb I 1 S 0 -1 P 1 transition frequency at 399 nm using an optical frequency comb[END_REF] as well as ours, we find a reasonable agreement, with error bars still covering nearly 400 MHz. Within the error bars, our results are compatible with [START_REF] Meggers | The first spectrum of ytterbium (Yb I)[END_REF][START_REF] Enomoto | Comparison of resonance frequencies of major atomic lines in 398-423 nm[END_REF][START_REF] Kleinert | Measurement of the Yb I 1 S 0 -1 P 1 transition frequency at 399 nm using an optical frequency comb[END_REF], and our WS8-2 measurement differs from [START_REF] Nizamani | Doppler-free Yb spectroscopy with the fluorescence spot technique[END_REF] by a few MHz only. By doing a simple averaging of all 6 values, we find ν174 = 751 526 574 ± 35 MHz, while a weighted mean [START_REF]Using the inverse of the squared uncertainties as weights[END_REF] yields ν174 = 751 526 533.5 ± 0.3 MHz.

Other isotopes

Table 1 displays our measurements of the ytterbium 1 S 0 → 1 P 1 transition absolute frequencies for the different isotopes and hyperfine transitions. The uncertainty of the wavemeters when measuring a frequency difference over the ∼GHz range is specified by the manufacturer to be equal to the absolute accuracy. For this reason, we do not find necessary to present a thorough comparison of our results with the isotopic shifts measurements found in literature, which have

■ ■ ○ ○ ★ ★ ▲ ▲ ◇ ◇ ▽ ▽ □ □

Meggers and Tech

(1978)

Das et al.

(

) 2005
Nizamani et al.

(

) 2010
Enomoto et al.

(

) 2016
Kleinert et al.

(

) 2016
This work

(WS/7)
This work Fig. 3. Absolute measurements of the 174 Yb 1 S 0 → 1 P 1 transition frequency. Blue square: [START_REF] Meggers | The first spectrum of ytterbium (Yb I)[END_REF]; yellow circle: [START_REF] Das | Absolute frequency measurements in Yb with 0.08 ppb uncertainty: Isotope shifts and hyperfine structure in the 399-nm 1 S 0 → 1 P 1 line[END_REF]; green star: [START_REF] Nizamani | Doppler-free Yb spectroscopy with the fluorescence spot technique[END_REF]; red upward pointing triangle: calculated from [START_REF] Enomoto | Comparison of resonance frequencies of major atomic lines in 398-423 nm[END_REF]; purple diamond: [START_REF] Kleinert | Measurement of the Yb I 1 S 0 -1 P 1 transition frequency at 399 nm using an optical frequency comb[END_REF]; brown downward pointing triangle: this work with WS/7; light blue square: this work with WS8-2. Frequencies are plotted relatively to the atomic transition in the 174 Yb isotope as measured in [START_REF] Meggers | The first spectrum of ytterbium (Yb I)[END_REF].

a much smaller dispersion than this uncertainty. Our results are nevertheless in very good agreement with all the isotopic shifts measurements found in [START_REF] Das | Absolute frequency measurements in Yb with 0.08 ppb uncertainty: Isotope shifts and hyperfine structure in the 399-nm 1 S 0 → 1 P 1 line[END_REF][START_REF] Nizamani | Doppler-free Yb spectroscopy with the fluorescence spot technique[END_REF][START_REF] Kleinert | Measurement of the Yb I 1 S 0 -1 P 1 transition frequency at 399 nm using an optical frequency comb[END_REF][START_REF] Deilamian | Isotope shifts and hyperfine splittings of the 398.8-nm Yb I line[END_REF][START_REF] Loftus | Optical double-resonance cooled-atom spectroscopy[END_REF][START_REF] Banerjee | Precise measurements of UV atomic lines: Hyperfine structure and isotope shifts in the 398.8 nm line of Yb[END_REF][START_REF] Wang | A novel method to measure the isotope shifts and hyperfine splittings of all ytterbium isotopes for a 399-nm transition[END_REF].

It should be noted that even though the wavemeters manufacturer recommends using the absolute uncertainty for frequency difference measurements, this seems to be a too conservative estimate. To verify this, we successively locked an 870 nm ECDL to adjacent teeth of a frequency comb. The teeth had a fixed spacing of about 250 MHz and stable well below the MHz level. By reading the frequency-doubled output of this laser with the wavemeters after each 250 MHz step, we checked that the uncertainty on frequency difference measurements was below 6 MHz for both the WS/7 and WS8-2 wavemeters across a ∼7 GHz range around 435 nm. We are unfortunately unable to perform this test around 399 nm to verify that our measurement could be used to estimate the isotope shifts at the few MHz level.

Sources of uncertainty

Measurement method and frequency scale

We record fluorescence spectra by acquiring a time trace of the fluorescence signal with an oscilloscope while simultaneously monitoring the laser frequency with the wavemeters. During a measurement the 100 s-long frequency scan is made over a range of 3 GHz, which induces uncalibrated non-linearity when using the laser PZT. In practice, we sample the laser frequency every ∼ 25 ms using the wavemeters and interpolate between the measurements to convert the oscilloscope time trace to a frequency scale. The error on the frequency scale origin from this conversion results in an uncertainty corresponding to one time step at maximum, and is thus evaluated to be below 0.8 MHz.

Systematic errors

In addition to the uncertainty inherent to the measurement method, different sources of error of physical origin can be identified, including frequency shifts and broadenings due to the laser intensity and linewidth, the residual magnetic fields and the laser polarization, and the Doppler effect.

First, the resonances in the recorded spectrum have no detectable amount of asymmetry. The position of their maxima stays within the error bar returned by the fit, which is found to be sub-MHz for all of them. It is therefore very conservative to assume an error of ±1 MHz from the combined effects leading to a broadening of the resonances, including non-linearity of the frequency measurement over the linewidth range.

The experiment, not protected by a magnetic shield, is subject to the earth field and environmental background fields. The residual magnetic field around the chamber was measured to have an amplitude below 1 G leading to a splitting of the transitions of less than 1.4 MHz for both odd and even isotopes. Moreover, the laser polarization is random. In these conditions, Zeeman splitting leads to a minor asymmetric broadening of the recorded resonances (natural linewidth of 28 MHz) and sub-MHz shifts of the maxima. For odd isotopes, quantum interferences and polarization effects [START_REF] Kleinert | Measurement of the Yb I 1 S 0 -1 P 1 transition frequency at 399 nm using an optical frequency comb[END_REF][START_REF] Brown | Quantum interference and light polarization effects in unresolvable atomic lines: Application to a precise measurement of the 6,7 Li D 2 lines[END_REF] can also lead to sub-MHz shifts of the resonances. In the particular case of the F = 5/2 → F = 3/2 hyperfine transition in 173 Yb, very specific polarizations of the laser with respect to the residual magnetic field would result in an optical pumping effect [START_REF] Mcknight | Comment on "laser cooling of 173 Yb for isotope separation and precision hyperfine spectroscopy[END_REF] and forbid detection of fluorescence for this transition. We do not claim in this work to measure separately transition frequencies for 172 Yb and 173 Yb F=5/2→F =3/2 , and the presence or the absence of the 173 Yb F=5/2→F =3/2 peak in the fluorescence spectrum can only induce sub-MHz level uncertainty on the common measurement largely dominated by the 172 Yb transition frequency value. In conclusion, residual magnetic fields and polarization effects are conservatively considered to induce an overall error of ±3 MHz on the transition frequency.

Doppler shifts can be caused by the non-perpendicularity of the laser beam relatively to the atomic beam direction, which is not perfectly known inside the vacuum chamber. The longitudinal velocity distribution of the atomic beam is centered around a non-zero value, which leads to a Doppler shift depending on the angle θ. From the temperature of the dispenser T ∼ 875 K, we can roughly estimate this shift [START_REF]Boltzmann distribution with most probable speed v at = 2k B T M , with k B the Boltzmann constant and M the mass of an ytterbium atom[END_REF] to v at sin θ λ ∼ sin θ × 725 MHz, with v at the most probable speed and λ the laser wavelength. This amounts to ∼ 13 MHz/ • at small angles. This sensitivity can thus induce a significant contribution to the systematic error budget. In order to reduce the uncertainty due to this effect, we recorded the fluorescence spectrum at different angles θ using the previously described retro-reflection mirror (section 2 and Fig. 1).

Figure 4 shows recorded fluorescence spectra for different values of the angle θ. For each θ value, each resonance peak is split into a doublet of separate resonances corresponding to opposite Doppler shifts. For clarity, spectra have been normalized to the spectrum obtained with θ = 1 • . A horizontal average has also been applied in order to smooth out the curves for better readability. For analysis, the data have been fitted with the sum of two translated spectra which include all isotopes lines (one for each way of the laser propagation). In the fitting function, each line is considered a Voigt profile, which reduces processing to an acceptable time, and the symmetrical Doppler shift of each of the two spectra is left as a common free parameter. Error bars returned by the fits are all sub-MHz regarding the Doppler shift results. These shifts are then plotted as a function of the angle θ and finally linearly fitted as a function of sin θ in order to find the position for which the laser beam is perpendicular to the most probable velocity, ie θ = 0 • . These results are presented in the inset of Fig. 4. The 1 • data point has been discarded in this last analysis even though the fitting process of the recorded spectrum seems to return a reliable result. The fit indicates a sin θ × (787 ± 5) MHz dependence of the Doppler shift. Given the error on the perpendicular position from the analysis (±0.3 • ) and the angular resolution of the rotatable breadboard (5'∼ 0.083 • ), we estimate the uncertainty on frequency measurements caused by this parameter to be ±4.3 MHz. Note that an asymmetrical Doppler broadening due to the same effect is to be expected on the final recorded spectrum, but this has also been minimized by the above procedure.

• • • • • • • • • • • • 0

Overall uncertainty

All pre-mentionned uncertainties have been summarized in Table 2 as well as the accuracy for each wavemeter. From these results, we deduce an overall uncertainty of ±60 MHz for the WS/7 measurement and ±30 MHz for the WS8-2. With all other sources of uncertainties having a very small contribution, the overall uncertainty is completely determined by the wavemeters. Using these uncertainties and deriving the weighted mean of our two results, the absolute frequency of the 1 S 0 → 1 P 1 transition of 174 Yb is found to be 751 526 537 ± 27 MHz. For other isotopes, the weighted averages have been indicated in Table 1.

Source of error

Uncertainty (±MHz)

Frequency scale conversion 0.8

Line broadening effects 1

Residual magnetic field and laser polarization 3

Laser/atomic beams non orthogonality 4.3

Wavemeter WS/7 60

Wavemeter WS8-2 30

Table 2. Summary of the conservative estimates on the uncertainties affecting the 1 S 0 → 1 P 1 transition measurement.

Conclusion

In this article, we have reported on two independent measurements of the 1 S 0 → 1 P 1 transition frequency for different isotopes of neutral ytterbium. The measurements were done by recording fluorescence spectra using two independent commercial wavemeters. We implemented a method based on counterpropagating beams that minimizes the uncertainty due to the Doppler shift, by identifying the angle for which the laser beam is perpendicular to the atomic beam. By presenting these results, we hope to reduce the existing uncertainty on this transition frequency, suffering from a non-negligible dispersion in the literature. Our results confirm in particular the most recent accurate measurement by Kleinert et al. [START_REF] Kleinert | Measurement of the Yb I 1 S 0 -1 P 1 transition frequency at 399 nm using an optical frequency comb[END_REF] from 2016. This might be useful to future starting experiments based on neutral and ionized ytterbium where knowing accurately this transition frequency is crucial.

Fig. 4 .

 4 Fig.[START_REF] Godun | Frequency ratio of two optical clock transitions in 171 Yb + and constraints on the time variation of fundamental constants[END_REF]. Fluorescence spectra of ytterbium atoms as a function of the laser frequency while retro-reflecting the probing beam at different angles:θ ∼ 1 • (blue), θ ∼ 3 • (yellow), θ ∼ 5 • (green), θ ∼ 7 • (red).Frequencies are plotted relatively to the atomic transition in the 174 Yb isotope. Inset: Doppler shift as a function of the angle θ and associated linear fit indicating the zero-shift angle for perpendicular configuration.

Table 1 .

 1 Absolute measurement of the frequency of the 1 S 0 → 1 P 1 transition depending on the isotope and the hyperfine transition. A single measurement is provided for 172 Yb and F = 5/2 → F = 3/2 in 173 Yb as the two resonances are not distinguishable.

			Transition frequency
	Isotope	WS/7 (MHz)	WS8-2 (MHz)	Weighted Mean (MHz)
	176 Yb	751 526 021(60) 751 526 031(30)	751 526 029(27)
	173 Yb	751 526 280(60) 751 526 292(30)	751 526 289(27)
	172 Yb & 173 Yb (F =3/2) 751 527 063(60) 751 527 072(30)	751 527 070(27)
	173 Yb (F =7/2)	751 527 122(60) 751 527 131(30)	751 527 129(27)
	171 Yb (F =3/2)	751 527 370(60) 751 527 378(30)	751 527 376(27)
	171 Yb (F =1/2)	751 527 689(60) 751 527 694(30)	751 527 693(27)
	170 Yb	751 527 725(60) 751 527 730(30)	751 527 729(27)

Acknowledgments

The authors would like to thank Emmanuel Bigler for designing and setting up the detection optics, Vincent Giordano for fruitful discussions about the experiment, Jacques Millo for constructive remarks and help with the optical frequency comb operation, and Rodolphe Boudot for his precious comments about the wavemeter calibration and the careful reading of the manuscript.

Funding

This work has been supported by the Agence Nationale de la Recherche (ANR-14-CE26-0031-01 MITICC, ANR-10-LABX-48-01 First-TF, and ANR-11-EQPX-0033 Oscillator-IMP), the Région Bourgogne Franche-Comté, the Centre National d'Études Spatiales and the EIPHI Graduate School (contract "ANR-17-EURE-0002").

Disclosures

The authors declare that there are no conflicts of interest related to this article.


David Rabus

Jean-Michel Friedt

Lilia Arapan

Simon Lamare
email: simon.lamare.pro@hotmail.fr

Marc Baque

Greǵoire Audouin

Fredeŕic Cheŕioux
email: frederic.cherioux@femto-st.fr

Sub-surface H 2 S detection by a Surface Acoustic Wave passive wireless sensor interrogated with a ground penetrating radar

Keywords: Surface Acoustic Waves, Hydrogen Sulfide, Ground Penetrating RADAR, Wireless Measurement, Functional Polymer

Long-term monitoring of organic pollutants in the soil is a major environmental challenge. We propose to meet this issue by the development of a polymer dedicated to selectively react with H2S, coating surface acoustic wave transducers designed as passive cooperative targets with the compound, and probing their response using Ground Penetrating RADAR, thus providing the capability to monitor the presence of H2S in the subsurface environment. The selectivity is brought by including lead(II) cation in a reticulated polymer matrix which can be deposited as a thin layer on a surface acoustic wave sensor. We demonstrate a signal enhancement mechanism in which water absorption magnifies the signal detection, making the sensor most sensitive to H2S in an underground environment saturated with moisture.

Industrial infrastructures or human activities lead in most cases to the pollution of the corresponding soils. Therefore, land restoration after its use is one of the major environmental challenges in many countries for the development of a sustainable society. [1] For instance, the detection of hydrogen sulfide (H2S) in the soil underneath the industrial infrastructures or in some buried industrial infrastructures, like pipelines, is of great importance because H2S is a byproduct of many industrial activities [2] which is lethal at concentrations as low as 500 ppm. [START_REF] Doujaiji | Hydrogen sulfide exposure in an adult male[END_REF] Many works have demonstrated the development of sensors with high selectivity, low detection limit and fast response, in order to monitor the presence of H2S in gas-phase. [4] The pollution of soils is mainly investigated by the classical drilling-treatment-analysis method. However, this method is expensive, slow and marred by errors. Therefore, development of powerful systems, aiming at monitoring the pollution in soils in order to replace the drilling-treatment-analysis, is required. [START_REF] Pamukcu | Underground Sensing[END_REF] Underground sensing systems could emerge as an efficient solution if they were (i) easy to operate, (ii) energy efficient, (iii) accurate and reliable for long-term with referencing capabilities, and (iv) sensitive to preemptive warnings of geo-and environmental hazards. [START_REF] Doujaiji | Hydrogen sulfide exposure in an adult male[END_REF] These milestones require disruptive technological solutions to improve both the sensors and their interrogating systems, in order to unambiguously monitor in-situ the presence of organic pollutants in the subsurface, even over several decades. Within this objective, surface acoustic waves (SAW) sensors have been identified during the last two decades as one of the most promising passive embedded sensor platforms to deliver temporally and spatially ondemand sensing capable for measuring an event signal over a very long period (even years) thanks to cumulative measurement capability with no hindrance from a sensing element local energy source finite lifetime. Indeed, while SAW sensors have been initially designed as radiofrequency analog processing components, [6] they have been widely used for direct detection chemical sensors [7] and as passive cooperative target for short range RADAR sensing systems. [8] Nevertheless, only few investigations use the ability of SAW chemical sensors to be probed through a wireless link. [9] Another challenge then concerns the development of an interrogating system able to monitor the signal backscattered by the acoustic transducer through the soil. As opposed to short range wireless measurement techniques such as RFID, magneto-acoustic or magneto-elastic measurements, [START_REF] Cai | A wireless, remote query magnetoelastic CO2 sensor[END_REF] or even inductively coupled short range interrogation we focus on long range propagative electromagnetic measurement techniques as implemented in RAdiofrequency Detection And Ranging (RADAR). Ground Penetrating RADAR (GPR) is a classical geophysical tool which is widely used to observe sub-surface interfaces, voids, leakages or buried objects (landmines, pipes etc.) associated with electrical conductivity variations. [11] However, up to now, the stability of some of GPR instruments time base prevents the fine measurement of the delay of the signal returned by the target (i. e. echo), mandatory stability of using GPR as interrogation system of a chemical SAW sensor. [12] In this paper, we develop the first complete system allowing the detection of H2S in subsurface environments, e.g. in the context of leak detection, with no claim for quantitative concentration measurement. To achieve this objective, we develop 1) a SAW sensor acting as passive cooperative target and the associated sensitive layer compatible with the transducer to detect the H2S and 2) a Ground Penetrating Radar (GPR), playing the role of the interrogating system of the SAW sensor and operating from the surface controlled by a dedicated software for extracting the sensor signature and hence the detection of H2S in soils, even at a depth of 1 meter.

Results

GPR and SAW sensor design

Surface acoustic wave transducers meet both requirements of sensitive gravimetric sensors as well-known from the field of biosensing [7c] and wireless interrogation capability. [8] Indeed, they do not require any local energy source close to the sensing elements. Unlike the RadioFrequency Identification (RFID) tag communicating with a reader by using impedance modulated digital sentences, and whose measurement capability is hindered by the burden of large analog to digital conversion power consumption and the need to reach the rectifier diode voltage threshold, limiting its measurements range, the SAW transducer is based on the linear inverse piezoelectric conversion of the incoming electromagnetic wave collected by an antenna connected to the interdigitated transducer (IDT) electrode bus in order to generate an acoustic wave (Figure 1). All electrodes are patterned by lift-off using standard clean room facilities including photolithography of a thin layer of 250 nm evaporated aluminum. This acoustic wave propagates on the single crystal piezoelectric substrate with a velocity dependent on the boundary conditions including adlayers of molecules coating the propagation path. Ideally, a reference path with no chemical sensitivity allows for compensating for correlated effects (range from RADAR to sensor, temperature) while a measurement path acoustic velocity is dependent on correlated effects as well as chemical compound concentration (Figure 1). In our implementation, the whole sensor chip is coated with the chemical sensing layer and selective patterning of the reference path is still under investigation. This implementation allows for rejecting RADAR to sensor range nevertheless by analyzing delay differences between two echoes. After reflecting this wave using Bragg mirrors patterned on the piezoelectric substrate, the linear direct piezoelectric effect converts the acoustic wave back to an electromagnetic wave: from a user perspective, the SAW sensor is a radiofrequency electrical dipole returning echoes delayed by a duration dependent on a chemical compound layer thickness adsorbed on the sensor surface. Therefore, the objective of the design of the SAW sensor is mainly to optimize the returned power to extract the signal from noise by the readout system. The optimization of piezoelectric substrate electromechanical coupling coefficient and the lengthening of the delays of the returned echoes beyond clutter are the two design strategies allowing to reach interrogation ranges of multiple meters in low-loss subsurface environments. [START_REF] Friedt | Surface acoustic wave devices as passive buried sensors[END_REF] In this work, as we aim at measuring ppm-range delay variations, we used a RADAR provided by Sensors and Software's SPIDAR control unit powering 100 and 200 MHz unshielded antennas because this apparatus possesses a very good stability of reference timescale clocking. We have indeed observed, on a one-hour long measurement, a standard deviation of 10 ps on a 0.3-microsecond echo delay difference, or a relative stability of 30 ppm at 2-second integration time falling to 10 ppm for a 1-minute integration time. The SAW sensor design includes an IDT, three mirrors and two contact pads. (Figure 1 and supporting information) The contact pads connected to the antenna are used for the electrical excitation of acoustic waves since the IDT induces an electric field in the piezoelectic substrate. Two parameters limit the possible values of the distances D between the three mirrors and IDT. D is given by D = T•c, with c the acoustic wave velocity dependent on the boundary conditions and hence the mass loading of the sensing layer and T the time delay. The lower value of D originates from the clutter suppression achieved by delaying the sensor echoes beyond all subsurface reflections assuming that the GPR receiver sensitivity only allows for recording signals reflected by interfaces located in soil at depths shallower than 100-m. This assumption assesses sensor echo time delay longer than one microsecond, which corresponds to D bigger than 1.5 mm, due to the typical soil permittivity around 2.25 or an electromagnetic velocity around 200 m/µs in soil and typical SAW velocities in the 3000 m/s range. The upper value for D is given by the acoustic propagation losses and sensor compactness which assesses D to be smaller than 4.5 mm (i. e. T below 3 microseconds). From the point of view of signal to noise ratio, the GPR is able to detect SAW sensor echo phase variations of at least few degrees. Indeed, classical time delay analysis using a spectral measurement, e.g. using a network analyzer, measures a phase rather than a time delay. If the system operates at 100 MHz, the SAW sensor acoustic path of 1.5 to 4.5 mm introduces a phase of 3.6•10 4 to 1.10 5 degrees. In a wireless measurement configuration using a pulsed RADAR, the phase measurement is either performed through the Fourier transform of the returned echoes or the cross-correlation of the echoes. [START_REF] Friedt | Surface acoustic wave devices as passive buried sensors[END_REF] Considering the signal to noise ratio of the GPR acquisition allows for 10-ppm phase variation measurements, only phase changes of a few degrees can be accurately measured.

A key-point consists in the suppression of the delay introduced by the path between the reader unit and the subsurface sensor. This is successfully achieved by using differential measurements: since an echo delay is the sum of the time needed for the wave to travel from the GPR to the sensor to which the varying acoustic delay is added, a differential measurement cancels the former to only yield a time delay difference T21=2•(D2-D1)/c. [START_REF] Shmaliy | Limiting phase errors of passive wireless SAW sensing with differential measurement[END_REF] Therefore, instead of only using Bragg Mirrors (M2 and M3) to return the echo from the active zone, another one (M1 in Figure 1) is deposited out of the chemical sensing area. This type of differential experiment requires an instrument with sampling rate or local oscillator stability better than these variations induced by the measurement (i. e. few ppm). In our experiments, D1 = 2.0 mm and D2 = 3.0 mm yielding delays of 1.0 and 1.5 µs, respectively, with a center frequency selected to the GPR operating frequency as shown in Figure 2.

Active sensitive layer

We aim at functionalizing the SAW reflective delay lines with a compound selective and sensitive to hydrogen sulfide. The classical approaches of functionalization of SAW devices are based on the deposition of a thin film or the grafting of a monolayer of molecules onto a piezoelectric substrate. [START_REF] Zhou | Human Immunoglobulin Adsorption Investigated by Means of Quartz Crystal Microbalance Dissipation, Atomic Force Microscopy, Surface Acoustic Wave, and Surface Plasmon Resonance Techniques[END_REF] Then, the exposure of SAW device to the targeted gas molecules lead to an ab-or adsorption (without any covalent link between analyte and the sensing films), on the sensitive layer. This ab-or adsorption promotes an increase of the film mass. This mass loading effect decreases the acoustic velocity and hence increases the time delay for a fixed path length, also observed as a phase decrease at a fixed delay. We aim at developing a tailored organic resist which 1) reacts selectively and irreversibly with H2S, 2) is compatible with the classical processes of microfabrication used in clean-room, 3) yields irreversible reaction processes for cumulative measurements compatible with periodic probing of the sensor response as opposed to continuous monitoring.

The sulfur atom in H2S molecule exhibits a unique reactivity with metal ions, like Pb 2+ or Zn 2+ , which inspired the well-known lead(II) acetate test paper. [START_REF] Cha | Sub-Parts-per-Million Hydrogen Sulfide Colorimetric Sensor: Lead Acetate Anchored Nanofibers toward Halitosis Diagnosis[END_REF] Therefore, we consider introducing lead(II) acetate functional groups in a resist formulation, in order to be consistent with classical cleanroom spin coating for batch sensor fabrication. As many resists used in micro-nanofabrication processes are based on the thermo-or photopolymerization of epoxides with nucleophiles (like alcohol or amine derivatives), we have first synthesized two a,w-dihydroxyalkyl lead(II)carboxylate complexes which are ended by two hydroxyl functions (see Scheme 1 and supporting information for experimental details).

Scheme 1. Synthesis of a,w-dihydroxyalkyl lead(II) carboxylate complexes with n = 1 or 2. This reaction is versatile and efficient. The ring-opening of a lactone in presence of a strong base and lead(II) nitrate gives quantitatively the expected lead(II) carboxylate ended by the OH functions. The length of the alkyl chains in final complexes is fixed by the choice of the starting lactone.

Then, in order to achieve the formation of the resist, we have mixed the tailored a,w-dihydroxyalkyl lead(II)carboxylate complexes (as monomer, 7 mmol) with 1,3-butadiene diepoxide as co-monomer (14 mmol), 1,3,5-trihydroxy benzene as reticulating agent (5 mmol) in a DMF solution (74 mmol). Scheme 2. Synthesis of the reticulated copolymer including lead(II) carboxylate sites as active centers for the sensing of H2S.

As the final co-polymer is expected to be highly reticulated, and consequently quite insoluble, we have directly deposited the mixture on the SAW devices. To achieve this functionalization, a drop of this solution was spincoated (acceleration: 1000 rpm/s, speed: 2000 rpm, duration: 10s) onto a wafer containing the SAW reflective delay lines. Finally, the film is annealed at 180°C for 10 min. This procedure leads to the formation of a solid thin film (thickness: 250 ±20 nm, measured by mechanical profiler) which covers entirely the wafer (see Figure S1 in supporting information). This simple and reproducible procedure is fully compatible with a collective fabrication process at the wafer scale, although selective patterning of the sensing area and leaving the reference area insensitive to chemical reaction is still under investigation, with a lift-off process having been demonstrated but requiring manual removal of the sacrificial resist layer.

Demonstration of the sensing-capability of the SAW-functionalized components

In a first experiment, we aim at demonstrating the capability of the tailored resist to detect efficiently H2S in a gas phase. A typical experiment is based on the exposure of a SAW sensor fully coated with the sensing layer to a flux of H2S and/or H2O, each one being in gas phase. The measurement is performed initially in a wired configuration, with a Rohde & Schwarz ZVC8 vector network analyzer, with the time domain measurement option collecting the frequency domain reflection coefficient and performing the inverse Fourier transform to provide the returned power and phase at the time delay of the various echoes. All following plots exhibit the phase difference between echoes as a function of measurement time.

Three experiments have been planned to investigate the sensing capability. The first one consists in the exposure of the SAW sensor (operating at 380 MHz) to a dry H2S atmosphere at room temperature.

As soon as the SAW sensor is exposed to the dry atmosphere of H2S (20% in dry N2), we observed a strong decrease of the phase (Figure 3). The phase variation, dPhi, is 50° after 800s, corresponding to a sensitivity of 700 ppm (Figure 3). This phase variation of 50° can be assigned to a mass loading during the exposition of the SAW sensor. On the one hand, the gravimetric sensitivity (S=(dPhi/Phi)•(A/Dm) with dPhi the phase variation around the phase rotation along the path Phi induced by a mass Dm variation over area A) of a Rayleigh wave is assumed to be in the S=100 cm 2 /g range. [START_REF] Du | An experimental study of Love-wave acoustic sensors operating in liquids[END_REF] This definition of sensitivity handling surface mass density and relative phase (or frequency) variation because it cancels the actual path length and is valid both for transmission and reflective delay lines. Considering a delay of t = 1.55 µs and an operating frequency of f= 380 MHz, the phase rotation along the path is Phi=360•t•f = 212040 degrees. The corresponding mass variation Dm per unit area is Dm/A = dPhi/(S•Phi) where dPhi is the phase variation (i. e. 50°, Figure 3) and the area A is 0.3•0.3 cm 2 (Figure 1). The numerical application hints a mass variation of about 2.4 µg/cm 2 , since the resist thickness is small enough to only act as a perturbative layer, or a hydrogen sulfide mass increment of 0.21 µg in the layer. On the other hand, the expected absorbed mass of hydrogen sulfide is estimated as 0.16 µg. Indeed, the molecular weight of the polymer is 692.6 g/mol (Scheme 2). Assuming a density of 1.4 g/cm 3 for this polymer, [START_REF] Wypych | Handbook of polymers[END_REF] the volume of sensitive layer of the SAW sensor is 250•10 - 7 •0.3•0.3 cm 3 , which correspond to 4.55 nmol. By considering the molecular weight of 34 g/mol for H2S, the layer loads an added mass of 0.16 µg. This value is very close to the expected 0.21 µg deduced from the acoustic signal analysis assuming a purely gravimetric effect. We conclude that in a dry saturated atmosphere, assuming all lead sensing sites react with hydrogen sulfide, that mass loading is sufficient to interpret the data.

Another set of experiment is made by using an air atmosphere saturated with water moisture at ambient pressure at 20°C and, then exposed to an atmosphere containing 20% of H2S and water saturated pressure at room temperature (Figure 4).

When the SAW sensor (operating at 380 MHz) is exposed to an atmosphere saturated with H2O, the phase variation is 12° (inset in Figure 4). Then, the exposure to H2S to an impressive phase variation of over 650° (or 9000 ppm), which is 13 times stronger than those observed in a dry atmosphere (Figure 3). Exposure of a SAW sensor to atmosphere saturated with H2O, then followed by H2S (20%) in a dry atmosphere. The phase variation due to the H2O is of 12° (inset) while the exposure to H2S leads to a phase variation of 650°.

Finally, another experiment is performed by only inverting the exposure of a SAW sensor to H2S then H2O (Figure 5) by comparison with the previous experiment.

Figure 5. Exposure of a SAW sensor to atmosphere saturated with H2S (20% in a dry atmosphere), then followed by atmosphere saturated with H2O. The phase variation due to the H2S is 40° while the exposure to H2O leads to a phase variation of 100°.

After exposure with dry H2S, the phase variation is of 40°, which is close to those observed in the first experiment (Figure 3). However, as soon as the SAW sensor is exposed to a water-saturated atmosphere, the phase variation is 100° (Figure 5). The analysis upon water loading is much more complex since the mass loading is accompanied with polymer structure change which might affect its Young modulus and hence slow down the wave upon softening. [START_REF] Johannsmann | The quartz crystal microbalance in soft matter research[END_REF]

Subsurface sensor measurement

The previous experiments demonstrated the capability of the SAW sensor to detect H2S in a gas phase with wired configuration and in a controlled environment. We aim at investigating the efficiency of this SAW sensor to detect H2S in a wireless configuration by using a GPR as interrogation system. To achieve this experiment, a wireless reflective delay line SAW sensor is inserted in polyvinyl chloride (PVC) pipes which are buried at different depths ranging from 10 to 100 cm in a sandbox (Figure S2). [12] The sandbox is located outside the building and no specific condition (pression, temperature, weather etc.) was defined in order to simulate a "real condition of use". The sensor is probed from the surface by recording once every 2-seconds a GPR A-scan collected by the Sensors and Software SPIDAR controller connected to PulseEkko emitter and transmitted unshielded 100 or 200 MHz antennas.

The sensor is initially inserted in the PVC pipe, and a baseline is reached after a few minutes as measured continuously by the GPR custom control software. Once the stable baseline is reached, H2S is generated in situ by the reaction of an aqueous solution of HCl acid with a FeS powder, leading to the formation of H2S and FeCl3. The H2S produced gas is pushed by a fan in the PVC tube in order to expose the SAW sensor to this gas. The exposure of the SAW sensor leads to a phase variation of 25° which is easily detected by the GPR through one-meter thick sand (Figure 6). Additional exposure does not lead to a phase variation with these experimental conditions. The repeatability of wireless measurements is checked by using two different 200-MHz sensors located at different depths (0.3 and 0.4 m, see Figures S3 andS4, respectively). The observed phase variations are 130 and 220°, respectively. These three experiments demonstrate that the combination developed in this work based on a tailored sensitive resist, a SAW sensor and a GPR, is efficient to detect H2S in a gas phase in a wireless configuration in soil and in real conditions (weather, temperature, humidity, etc.).

Discussion

On the basis of all experimental data, we demonstrate that the complete system described in this article is efficient to detect in sub-surface the presence of H2S in a gas phase. All phase variations observed are significant (above 25° phase shift or 2000 ppm) with respect to the GPR noise level (6° phase fluctuation) thanks to the mass loading enhancement induced by the presence of moisture, with phase variations consistent with those found in the literature using other mechanisms. Indeed, 3 kHz frequency variation of a 60 MHz SAW device with a gap between IDTs of 2.5 mm, leading to an equivalent phase shift of 900 degrees was reported. [START_REF] Vetelino | Hydrogen sulfide surface acoustic wave gas detector[END_REF] The origin of this phenomenon is explained as follows: the exposure of a lead(II) acetate derivative to H2S leads to the formation of PbS and carboxylic derivative (Figure 7). The formation of PbS explains the color variation in the well-known lead(II) acetate test paper [START_REF] Cha | Sub-Parts-per-Million Hydrogen Sulfide Colorimetric Sensor: Lead Acetate Anchored Nanofibers toward Halitosis Diagnosis[END_REF]21]. In addition, as water has no role in this polymer matrix breaking process, the presence or the absence of water (or moisture) alone is hardly detected by the sensor as observed in Figure 4 before injection of H2S. However, when a sensor is exposed to H2S followed by H2O or H2O followed by H2S, we observe a strong enhancement of the resulting phase variation, respectively 600° and 100°. The enhancement mechanism is explained by the depolymerization of the resist when exposed to H2S. Indeed, the starting resist is initially compact (due to the reticulation) and hydrophobic, which prevents interaction with water (moisture). However, upon reaction with H2S, the lead(II) cations which are involved in the polymer chains are cleaved to produce PbS clusters and short oligomers ended by several carboxylic moieties (Figure 8). These carboxylic groups are hydrophilic due to the formation of multiple hydrogen bonds with water. As the skeleton of polymers is strongly altered and now turns to be strongly hydrophilic, water molecules are absorbed in the organic layer, which acts as a sponge, leading to a significant mass loading enhancement and hence slower acoustic wave, inducing the observed phase variations. This important phase enhancement is useful for the wireless detection in which signal to noise ratio is degraded by radiofrequency signal attenuation when propagating through soil (here sand) leading to poorer detection limits than in the wired configuration. Eventually, as expected from the irreversible reaction of H2S with the Pb 2+ cations, if the initial exposure of the SAW sensor induces complete transformation of all active sites, further exposure does not lead to additional phase variation (Figure 5, date 1000 s).

Conclusion

The combination of a tailored sensitive resist, a SAW sensor and a GPR, is used as an efficient system to detect H2S in a gas phase in a wireless configuration in soil and in real conditions (weather, temperature, humidity, etc.). The strong phase variation, required to compensate for radiofrequency signal attenuation when propagating through soil, is based on an enhancement of the mass loading due to the presence of water or moisture. This remarkable mass loading is justified by the initial design of the resist which was selected in order to induce its depolymerization during the exposure to H2S. As this reaction is irreversible and the sensor is battery-less, it is an ideal candidate for long-term cumulative detection of the presence of H2S in subsurface: phase shifts accumulate as pollutants flow over the sensor and even if sensor monitoring is periodic and not continuous, a punctual pollutant spike is recorded and will be measured during the next sensor probing session. This new strategy paves the way of the developments of sensors, and their interrogation by GPR, for land preservation of industrial sites.

Figure 1 .

 1 Figure 1. SAW sensor layout. The interdigitated transducer (IDT, with w = 3000 µm, a = 2p and p = 20 µm or 10 µm for a sensor working at 100 or 200 MHz, respectively) is located between three Bragg mirrors (M1, M2 and M3) positioned at different distances (D1, D2, D3) from the transducer in order to introduce different delays T. The contact pads are used for the electrical activation of IDT. M1 acts for the reference measurement. The active area is used for the deposition of the sensitive layer, M2 and M3 being used for the detection measurements.

Figure 2 .

 2 Figure 2. Frequency domain (top) and time domain (bottom) characterization of the reflective acoustic delay line used in the experiment depicted in this article. The red vertical markers on the time domain characteristics guide the eye to the echo delays analyzed, at 1.0 and 1.55 µs.

Figure 3 .

 3 Figure 3. Exposure of a SAW sensor to a dry atmosphere of H2S.

Figure 4 .

 4 Figure 4. Exposure of a SAW sensor to atmosphere saturated with H2O, then followed by H2S (20%) in a dry atmosphere. The phase variation due to the H2O is of 12° (inset) while the exposure to H2S leads to a phase variation of 650°.

Figure 6 .

 6 Figure 6. 100-MHz sensor located 1-m deep in sand: H2S is generated at times 570 s and 1000 s. The first exposure to H2S leads to a phase variation of 25°. The second exposure does not induce a phase variation indicating that the reactive sites are saturated upon the initial exposure.

Figure 7 .

 7 Figure 7. Reaction of H2S with a molecular lead(II) acetate derivative.

Figure 8 .

 8 Figure 8. Depolymerisation of copolymer leading to the creation of hydrophilic functions (carboxylic acid) and lead sulfide cluster due to the exposure to H2S.

This work is supported by the French National Agency through the project UNDERGROUND (ANR-17-CE24-0037). The authors also acknowledge partial support of the French RENATECH network and its FEMTO-ST technological facility. The authors thank S. Benchabane, L. Morgenthaler and P. Tanguy (FEMTO-ST), F. Geǵot (SENSeOR), J. R. Ordonez-Varela, A. Credoz, A. Le Beulze and M.-F. Bennassy (TOTAL) for fruitful discussions.

AUTHOR INFORMATION

Synthesis of complexes 53 mmol of g-butyrolactone (or e-caprolactone) and 53 mmol of KOH are dissolved in 150ml of ethanol and 30 ml of H2O. The mixture is heated at reflux for 1h. Then, 26 mmol of Pb(NO3)2 in 40 ml of H2O are added. Then, the mixture is heated at reflux for 20 min. The solvents were removed under reduced pressure. The crude solid is washed in 50 ml of hot ethanol, to give the expected complexes.

Lead(II) 4-hydroxy-butanoate 1 H NMR (300 MHz, D2O) δ = 3.50 (t, J=7.4Hz, 4H), 2.14 (t, J=7.4Hz, 4H), 1.70 (q, J=7.4HZ, 4H). 13 C NMR (75 MHz, D2O) δ = 183.1 (C=O), 61.6 (C-OH), 34.1, 28.4.

Lead(II) 6-hydroxy-hexanoate

1 H NMR (300 MHz, D2O) δ = 3.57 (t, J=6.5Hz, 4H), 2.19 (t, J=7.40Hz, 4H), 1.66-1.55 (m, 8H), 1.66-1.55-1.50 (m, 4H). 13 C NMR (75 MHz, D2O) δ = 184.1 (C=O), 61.5 (C-OH), 38. 3, 31.1, 25.0, 24.9.

SAW device fabrication

Ground Penetrating RADAR operating frequency is a tradeoff between antenna dimensions and interrogation range: reaching the 1-3 m depth targeted to reach the water table or the sub-surface infrastructure requires operating in the 100 to 400 MHz range. Thus, YXl/128° lithium niobate acoustic delay lines were manufactured. With a wavelength of 40 µm and a metal thickness of 2.5% wavelength, the L = 2 to 3 mm long acoustic path introduces a two way trip phase rotation of 2xLx360xf/c =36180 to 54270° phase rotation or a differential phase accumulation of 18090° considering the Rayleigh wave velocity c = 3980 m/s and an operating frequency of f = 100 MHz, or a differential phase rotation of 72360° for the same geometry of a sensor operating at 400 MHz. In order to avoid damping the Rayleigh wave, the organic sensing layer will be designed as thin with respect to the wavelength, aimed at varying the acoustic velocity without excessive acoustic losses: a thickness of about 1% wavelength is considered, yielding 100 to 500 nm thick adlayers for devices operating in the 100 to 400 MHz range.


Design of piezoelectric actuators by optimizing the electrodes topology

Keywords: piezoelectric actuators, topology optimization, polarity optimization

Piezoelectric materials based actuators are highly recognized for the development of microrobotic systems thanks to their high bandwidth, high resolution and high force density. However, one of their main drawback is the low relative stroke (0.1% of actuator's size) that limits the actuator motion range. Overcoming this limitation is challenging but would increase the achievable working space of microrobots. In this paper, topology optimization method is used to maximize the actuator stroke. Instead of optimizing only the material density, we also consider the optimization of the electrodes polarity. This approach allows to combine both material expansion and compression in order to increase the actuator output displacement. To demonstrate this approach, two actuators were designed starting from a full domain considered as a basic reference piezoelectric actuator. The first design considered only the density optimization while the second one took into account the optimization of the topology of electrodes. Both simulation and experiment showed a good agreement between the obtained designs and the fabricated prototypes. The result revealed that the optimized design with polarization has an improved factors of 2.8 and 2.02 compared to full plate and actuator without electrodes optimization, respectively.

I. INTRODUCTION

Accurate and precise positioning is a key element to enhance the capability of microrobotic systems to manipulate and characterize objects in the microworld, i.e. the world of submillimetric objects [START_REF] Régnier | Microrobotics for Micromanipulation[END_REF]. The literature tells us how smart or active materials are widely recognized to the development of actuators for these microrobots and improve their accuracy [START_REF] Kortschack | Smart materials for actuation in microrobots[END_REF]. These active materials are capable to change their physical inherent properties in response to an external stimuli like electrical field, magnetic field, thermal gradient or variation. Such feature allowed them to be largely employed as easily integrated and miniaturized actuators in the above microrobotic systems. Among these materials, piezoelectric materials are of particular interest since they can be used as an elementary actuator with high displacement resolution, large output force, high dynamics response and significant scaling-down possibilities [START_REF] Kortschack | Smart materials for actuation in microrobots[END_REF] [3] [START_REF] Mcclintock | The millidelta: A high-bandwidth, high-precision, millimeter-scale delta robot[END_REF]. However, due to their crystalline arrangement, they provide a low relative deformation (0.1% of actuator's size) that restricts their motion range [START_REF] Grossard | Mechanical and control-oriented design of a monolithic piezoelectric microgripper using a new topological optimization method[END_REF]. Overcoming this limitation would increase the achievable working spcaes of the microrobotic systems in which they are integrated. As reviewed in [START_REF] Wang | A survey of piezoelectric actuators with long working stroke in recent years: Classifications, principles, connections and distinctions[END_REF], empirical and systematic approaches were investigated to the development of piezoelectric actuators. In particular several solutions have been proposed to increase their displacements. These solutions can be classified in four major categories: (i) staking piezoelectric layers, (ii) hybridization, (iii) amplification mechanisms, (iv) structural optimization.

The first class represents the trivial solution that operates on the shape, the polarity and the geometric arrangement of the actuator. It leads to piezostacks or multimorph configurations [START_REF] Devoe | Modeling and optimal design of piezoelectric cantilever microactuators[END_REF]- [START_REF] Huajun Fang | Modeling and design optimization of large-deflection piezoelectric folded cantilever microactuators[END_REF]. Piezostacks are the result of stacking several layers in serial in order to increase the output displacement. However to produce reasonable output displacement, this configuration leads to long and non-compact actuators. Multimorphs consist of several piezoelectric cantilevers layers. They provide larger displacement than piezostacks but at the expense of the generated force. In the hybridization class, the actuator is based on the combination of piezoelectric effect and an additional principle. Thermopiezoelectric actuators were mainly devoloped to this aim [START_REF] Rakotondrabe | Development and dynamic modeling of a new hybrid thermo-piezoelectric micro-actuator[END_REF], [START_REF] Rakotondrabe | Development and force/position control of a new hybrid thermo-piezoelectric microgripper dedicated to micromanipulation tasks[END_REF]. The third class operates on external mechanism to amplify the actuator stroke. Known as amplifying mechanism, several structures have been investigated. The widespread ones are flexural lever structure [START_REF] Uchino | Electrostrictive actuators: Materials and applications[END_REF] and flextensional mechanisms based on empirical or optimized structures such as Moonie [START_REF] Sugawara | Metal-ceramic composite actuators[END_REF] [START_REF] Dogan | High displacement ceramic metal composite actuators (moonies)[END_REF], Cymbal [START_REF] Dogan | Composite piezoelectric transducer with truncated conical endcaps "cymbal[END_REF], nested rhombus structure [START_REF] Ueda | Large effective-strain piezoelectric actuators using nested cellular architecture with exponential strain amplification mechanisms[END_REF] and optimal amplification mechanisms based on passive material [START_REF] Canfield | Topology optimization of compliant mechanical amplifiers for piezoelectric actuators[END_REF], [START_REF] Schlinquer | Displacement amplifier mechanism for piezoelectric actuators design using simp topology optimization approach[END_REF]. Despite their potential, these approaches only deal with geometrical or topological considerations of the amplification mechanism surrounding a predetermined actuator. In general, the amplifying mechanisms which are augmented to the actuators make the whole design bulky and heavy and inappropriate for small scale or micro dimension applications. Finally, the fourth class (structural optimization) is a systematic approach and includes parametric and topology optimization that consider piezoelectric materials within the optimization process [START_REF] Khadraoui | Optimal design of piezoelectric cantilevered actuators with guaranteed performances by using interval techniques[END_REF]- [START_REF] Homayouni-Amlashi | Topology optimization of 2dof piezoelectric plate energy harvester under external in-plane force[END_REF]. The approach consists to find the actuator's optimal dimensions or shapes on the basis of intervals techniques [START_REF] Rakotondrabe | Performances inclusion for stable interval systemstechniques[END_REF] or with topology optimization methods [START_REF] Bendsoe | Topology Optimization: Theory, Methods and Applications[END_REF] respectively. The advantage of the topology optimization is that there is no restriction on the shape of the actuator. Therefore, the whole structure of the actuator can be optimized. Although topology optimization provides a great degree of freedom in terms of structural optimization, the obtained results might bring complexity in terms of fabrication of the designed structures. However this challenge can be tackled by defining proper optimization parameters as it is explained in [START_REF] Bendsoe | Topology Optimization: Theory, Methods and Applications[END_REF]. One of the used topology optimization methods to integrate piezoelectric materials is the SIMP (Solid Isotropic Material with Penalization) method [START_REF] Kögl | Topology optimization of smart structures: design of piezoelectric plate and shell actuators[END_REF]. For instance it was successfully applied to design piezoelectric structures such as gripper [START_REF] Ruiz | Optimal design of a microgripper-type actuator based on AlN/Si heterogeneous bimorph[END_REF]. Designed as a monolithic structure, the achievement has shown a quite good performance by maximizing the in-plane displacement at the output of the gripper, whilst the out-of-plane bending of the gripper's jaw was suppressed. In spite of these interesting and promising achievements, this approach considers only the piezoelectric material distribution while the polarity of electrodes is kept fix.

In this paper, we propose to design piezoelectric actuators with the SIMP method. In addition, not only the density is optimized but we propose a modification of the algorithm such that the structure polarity is also considered. Such consideration is useful for two reasons: (i) it provides several degrees of freedom in terms of optimization where two variables (density and polarity) can be considered, and (ii) it allows to combine material expansion and compression in order to increase the output displacement. Therefore, the proposed approach will optimize the structure of the actuator without augmenting any additional structure. In this case, the optimized designs can be miniaturized for microrobotic applications. To validate the proposed approach, two actuators will be designed where their output displacements are maximized. The design of each actuator starts from a full domain plate considered as a basic reference actuator. Within this domain, the first design is based on density optimization with a fix polarization of the electrodes, whilst the second design is based on both density and polarity optimization. The result will reveal that when polarity is considered, the outcome actuator provides 2.8 and 2.02 as amplification factors in term of displacement compared to full plate and actuator without electrode optimization, respectively. Simulations and experiments are performed to validate the obtained designs and the fabricated prototypes.

II. TOPOLOGY OPTIMIZATION O PIEZOELECTRIC STRUCTURES

This section briefly reminds (i) the principle of topology optimization (ii) the main major steps allowing to derive finite element model of a thin piezoelectric plate and (iii) the integration of piezoeletric model within the topology optimization method.

A. Topology optimization

Topology optimization especially SIMP method is a design methodology aiming to find an optimal structure within a limited design domain [START_REF] Bendsoe | Topology Optimization: Theory, Methods and Applications[END_REF]. Based on material distribution, it consists to minimize systematically an objective function while subjected to one or several constraints. To this end, the method combines finite element formulation and penalization power law to make material intermediate density unattractive and therefore avoid the 0-1 problem of classical topology optimization processes. Presented in [START_REF] Bendsoe | Topology Optimization: Theory, Methods and Applications[END_REF], SIMP method is easy to implement and suitable for the design of passive structures. Since becoming a conceptual tool for structural design, it has been successfully applied for topological design of active structures in particular piezoelectric structures [START_REF] Homayouni-Amlashi | Multi directional piezoelectric plate energy harvesters designed by topology optimization algorithm[END_REF], [START_REF] Homayouni-Amlashi | Topology optimization of 2dof piezoelectric plate energy harvester under external in-plane force[END_REF]. This methodology made possible the consideration of the physics of the material within the optimization problem. To do so, the finite element model of the piezoelectric structure is required which we provide in the next section.

B. Piezoelectric modeling of 2D thin plate

As we target the design of 2D actuators, a plate model is considered in this paper. It consists of a thin piezoelectric layer sandwiched between two electrodes as illustrated in Fig. 1. This configuration is widely used as a framework for the design of planar active mechanism. To derive its model, we first assumed -→ the normal axis of the thin plate which is considered to be the axis of polarization as well.

Then, by applying the plane-stress theory [START_REF] Chawla | Mechanical Behavior of Materials[END_REF], the nominal stresses perpendicular to the xy plane are zero. We chose a 4mm class tetragonal crystal of piezoelectric materials [START_REF] Piefort | Finite element modelling of piezoelectric active structures[END_REF] to derive the corresponding model. This class includes most of the piezoelectric material in particular the well-known PZT material [START_REF] Chung | Microstructural Evolution in Lead Zirconate Titanate (PZT) Piezoelectric Ceramics[END_REF]. According to these hypothesis, the "e" form constitutive equations from [START_REF] Standard | Ieee standard on piezoelectricity[END_REF] are simplified to obtain the following new set of relations:

           1 2 6 3            =         1 2 0 - 2 1 0 - 0 0 6 0 0                    1 2 6 3            (1)
where:

1 = 11 - ( 13 ) 2 33 , 2 = 12 - ( 13 ) 2 33 , 6 = 66 (2) = 31 -33 13 33 , = 33 + ( 33 ) 2 33 (3)
The mechanical coefficients of the tensor given in equation (1) can be rewritten by identifying them with the plane stress isotropic tensor [START_REF] Chawla | Mechanical Behavior of Materials[END_REF]:

[ ] = 1 -2       1 0 1 0 0 0 1- 2       (4)
Therefore, we express the mechanical stiffness tensor of the plan-stress form of constitutive equation in (1) by Young's modulus and Poisson's ratio. Considering the 4mm piezoelectric class property where 6 = 1 2 ( 1 -2 ) [START_REF] Standard | Ieee standard on piezoelectricity[END_REF], the Young's modulus and Poisson's ratio according to 1 and 2 can be expressed as follows:

= 2 1 -2 2 1 , = 2 1 (5)
Taking advantage of planar isotropy of 4mm crystal, the planar Young's modulus and Poisson's ratio of this class of piezoelectric material can be calculated. This is an essential step for the normalization process which will be discussed in the next sections. In the rest of the paper, we consider the commercial piezoelectric material PSI-5H4E from Piezo Systems Inc. By substituting its numerical values, we obtain = 66.2 GPa and = 0.298. For the piezoelectric and dielectric coefficients, we obtain = -28.30 C.m -1 and / 0 = 1932 with 0 the vacuum's permittivity.

C. Finite element model of thin piezo plate

The finite element model of the piezoelectric plate can be derived by following the general approach given in [START_REF] Lerch | Simulation of piezoelectric devices by two-and threedimensional finite elements[END_REF]. It consists to compute the stiffness matrices using the thin plate model described previously and results in the equilibrium equation of the piezoelectric material. For this purpose, we extend the bilinear 4 nodes element to take into account the piezoelectric behavior. In this case, four classical nodes model the mechanical behavior while two additional nodes model the electrical degree of freedom along the -→ that correspond to the electrodes. Each element is considered as a flat cuboid with side length and thickness ℎ. Based on this configuration, interpolation functions for both displacement and electrical potential are defined as follows:

= ( ± 2 ) ( ± 2 ) 4 2 , = ℎ ± 2 2ℎ (6)
where and index each node, either mechanical or electrical. Based on these functions, we compute the mechanical and piezoelectric stiffness matrices following the procedure described in [START_REF] Lerch | Simulation of piezoelectric devices by two-and threedimensional finite elements[END_REF]:

[ ] = 0 , 0 = ℎ 1 -2 (7) = 2 -1 -1 1 -1 1 1 -1 1 1 1 -1 1 -1 -1 1 -1 (8)
with is well-known dimensionless eight by eight matrix (see [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF] for more details) and 0 is the normalization factor which was calculated with the help of the derived Young modules and Poisson's ratio. Using these two matrices, the global finite element system for actuation is expressed as:

{ } = [ ] { } + { } (9)
where and are respectively the nodal force and displacement vectors. Unlike passive material, here is null since the force is generated by the material itself, i.e no external load force. The above equation can be simplified knowing that { } = 1 2 and the potential = 1 -2 :

( { }) = 2 -1 -1 1 -1 1 1 -1 1 = 0 { } ℎ 0 = 2 ( 10
)
where 0 is the normalization factor of the piezoelectric matrix. With this simplification, equation ( 9) can be written as:

{0} = 0 { } + 0 { } (11)
Here 0 stands for a mechanical stiffness while 0 represents the internal force generated by the material when subjected to voltage. For simplification, we denote in the rest of the paper as [ ] and { } as { }, where subscript stands for normalization. This last term can be interpreted as the acting of piezoelectric forces on the element. This leads to a simple and compact finite element equation:

0 [ ] { } = -0 { } (12)
By analogy to pure mechanical behavior, this finite element model defines the equilibrium equation of a piezoelectric material which is mandatory to solve the topology optimization problem. Finally, the displacement vector is normalized for numerical considerations:

0 [ ] { } = -0 { } ⇔ [ ] (- 0 0 ){ } = { } (13) ⇔ [ ] (- 1 0 ){ } = { } (14) ⇔ [ ] { } = { } (15)
Such normalization has two important consequences: (i) the equation is independent from the material coefficients, domain dimensions and input voltage, which greatly improves the element versatility. The physical displacement can still be obtained using the elemental displacement 0 .

It is worth to notice that we are still dependent on the Poisson's ratio through the stiffness matrix, (ii) each member of equation ( 15) is of an order of magnitude between -2 and 0 that helps to reduce rounding error when resolving the linear system described by this equilibrium equation.

D. Integration of piezoelectric material within topology optimization method

Based on the above-derived finite element model, two modifications are brought in order to implement this model within SIMP method. The first modification concerns the extension of the classical penalization scheme that considers both density and polarity variables, whilst the second relies on the reformulation of the objective function and its derivation in respect to these two variables.

1) SIMP penalization scheme extension: Following the original penalization scheme presented in [START_REF] Bendsoe | Topology Optimization: Theory, Methods and Applications[END_REF], several extensions have been proposed to include piezoelectric material [START_REF] Kögl | Topology optimization of smart structures: design of piezoelectric plate and shell actuators[END_REF], [START_REF] Toledo | Design and characterization of in-plane piezoelectric microactuators[END_REF]- [START_REF] Mello | Design of quasi-static piezoelectric plate based transducers by using topology optimization[END_REF]. Most of them were applied on the material properties tensor

, [ ] and . However, these schemes penalize each coefficient of the tensors at the expense of computational time. To tackle this limitation, we propose to apply the penalization directly as a coefficient in front of the normalized stiffness matrix and force vector. This leads to a simple penalization scheme which can be written according to the material density and polarization direction value as follows:

[ ] ( ) = ( ) [ ] (16) { }( ) = ( , ){ } (17) with ( ) = ( - ) + (18) ( , ) = ( - ) (2 -1) + ( 19
)
where and are the stiffness and force penalization coefficients and and are the minimum stiffness and force factors. These minimum values are defined to avoid the singularity problems during the optimization [START_REF] Bendsoe | Topology Optimization: Theory, Methods and Applications[END_REF] [40] [START_REF] Bendsøe | Material interpolation schemes in topology optimization[END_REF]. The consideration of polarization direction as a variable is introduced in a methodology known as "piezoelectric material with penalization and polarization" (PEMAP-P) [START_REF] Kögl | Topology optimization of smart structures: design of piezoelectric plate and shell actuators[END_REF].

To apply this new penalization scheme, all parameters must be known. For , , and , we use the recommended numerical values ( = = 1 and = = 10 -9 ) given in [START_REF] Bendsoe | Topology Optimization: Theory, Methods and Applications[END_REF]. Material Density ( ) and polarization direction ( ) can have values between zero and one. By introduction of penalization coefficients ( , ), the material density ( ) will steer to zero (void) or one (material) during the optimization iterations. However, the penalization coefficients are difficult to estimate since there is no analogy with pure mechanical behavior. Beside, the fact that the penalization coefficients need to be higher than one [START_REF] Bendsoe | Topology Optimization: Theory, Methods and Applications[END_REF], there is no computation of those coefficients for plane stress heterogeneous piezoelectric material to the knowledge of the authors. This challenging issue remains an active research domain as shown in [START_REF] Bisegna | Variational bounds for the overall properties of piezoelectric composites[END_REF]- [START_REF] Ouafi | Analytical and semianalytical modeling of effective moduli bounds: Application to transversely isotropic piezoelectric materials[END_REF]. In this work we set both penalization factors constant to a value of three and no penalization coefficient is defined for the polarization direction since the best value for it is one [START_REF] Kögl | Topology optimization of smart structures: design of piezoelectric plate and shell actuators[END_REF]. Giving rise to consistent results, we did not investigate further this issue.

2) Objective function reformulation: Optimizing piezoelectric actuators can be done following several objectives like electric-mechanical conversion [START_REF] Carbonari | Topology optimization design of functionally graded bimorph-type piezoelectric actuators[END_REF] or optimal displacement [START_REF] Kögl | Topology optimization of smart structures: design of piezoelectric plate and shell actuators[END_REF]. For simplicity, only the latter is considered in this paper. Following the classical mechanical approach reported in [START_REF] Bendsoe | Topology Optimization: Theory, Methods and Applications[END_REF], we rewrote the objective as a minimization of an objective function = -=where is a vector with a value of one that corresponds to the output displacement node and zero otherwise. However, while the variation of this objective function reflects the evolution of the optimization process, its value holds no meaning due to the normalization of the finite element system. The elementary displacement 0 could be used, but this would nullify the versatility of the approach. Instead, we propose to normalize the objective function using its value when the domain is not yet optimized. This initial objective function value is now named 0 . This has an immediate consequence that allows to understand the objective as how much more displacement can be obtained by the optimized design comparing to full domain. This new objective function can be written as = 0 =where the vector includes the normalization factor 0 .

To perform a gradient based optimization, the sensitivities of this objective function with respect to material density ( ) and polarization ( ) are calculated using the well-known adjoint method. In addition, we wrote the system of equation ( 15) as =and the objective function as =in order to simplify the notations.

= (- + Λ ( + )) = ((-+ Λ ) + Λ ) = Λ + (-+ Λ ) + Λ = Λ ( + ) (20)
We obtain the last simplification by specifying -+ Λ = 0. This system can be rewritten as Λ = using the property =

, hence Λ can be seen as a displacement field resulting from a dummy force which is applied at the output displacement node [START_REF] Bendsoe | Topology Optimization: Theory, Methods and Applications[END_REF]. Following the same procedure, the sensitivity of objective function with respect to polarization can be obtained:

= Λ ( ) (21)
The main difference with pure mechanics is the dependency on the sensitivities of the piezoelectric force and . The stiffness matrix sensibility can be directly computed as follows:

= ( [ ]) = [ ] = ( - ) -1 [ ] (22)
The same procedure is applied for the piezoelectric force sensibility:

= ( - ) (2 -1) -1 { } (23) = 2( - ) { } (24)
This last step complete the integration process of the piezoelectric material within the SIMP method. In terms of programming, we re-use and fit the 88 lines code [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF] in order to include the 2D finite element presented above. Then, the OC (Optimal Criteria) is substituted with the MMA (Method of Moving Asymptotes) in order to improve the stability and performances of the optimization algorithm.

III. O

This section focuses on the design of one degree of freedom piezoelectric actuators using the approach described previously. Two actuators are designed in order to highlight the importance of optimizing the electrodes topology. The first actuator was based on density optimization with uniform electrodes while the second is based on the combination of density and polarity optimization. Thus the objective is to maximize the displacement of this point when optimizing the design. According to section II-D.2, we proceed by modeling this output as a force in the opposite direction and an artificial spring with low stiffness in comparison to piezoelectric material itself, as shown in red in Fig. 2.b. This is a well known procedure in designing the compliant mechanisms and further details can be found in [START_REF] Bendsoe | Topology Optimization: Theory, Methods and Applications[END_REF]. We consider the stiffness of the spring to be 0.01 of the derived young modulus of the piezoelectric material in section II-B. Furthermore, we model the clamping of the plate by adding supports at the bottom edge of the domain (blue boundary). At this stage, these specification are the same for both designs. For the first design a fixed uniform voltage over the upper electrode is considered and the lower electrode is grounded. In contrast, the second design considers electrodes polarity as an optimization variable. In this case, the electrodes can be set to ground, positive or negative potential depending on the optimization process. Finally, we take advantage of the symmetry to have the entire working domain in order to reduce the computation time and memory usage.

To solve this problem, we start by discretizing the design domain using the finite element approach presented in Sec. II. The size of each element is set to 200 µm, resulting in a mesh of 125×250 elements. Finally, the volume fraction is set to 0.3, meaning that only 30% of the full domain will be retained for the optimized design.

B. Post-processing

The obtained designs after optimization can be seen in Fig. 3.a and Fig. 4.a. The efficiencies (output mechanical energy/input electrical energy [START_REF] Zheng | Topology optimization of energy harvesting devices using piezoelectric materials[END_REF]) of the full piezoelectric plate, optimized design without polarity optimization and with polarity optimization are 0.482, 0.474 and 0.445 respectively. A negligible decrease in the efficiency of the optimized design can be seen. However, the details of force and displacement should be investigated after post processing. Although the raw designs are already well defined, we apply a classic post-processing approach to ease the exportation of the design as CAD model. First, the resolution is increased by a factor of three using bilinear interpolation. This operation allows to perform an anti-aliasing effect by smoothing the boundaries for cleaner exportation results. Next, a volume preserving threshold is performed to convert the intermediate binary densities, while maintaining the global volume of the domain. This is done using a simple dichotomy approach to find the optimal threshold. Doing so, it is easy to export the design to CAD software following the afterwards steps: (i) perform a contour detection of the post processed design under Matlab, (ii) export the resulted contour under Space Claim software from ANSYS (ANSYS, Inc, Canonsburg, USA) as a cloud of point, (iii) interpolate the sketch using splines and extrude the thickness of the piezoelectric plate. Then, with the help of the obtained CAD files, the performances of the optimized designs are analyzed with the commercial finite element software COMSOL Multiphysics (COMSOL, Inc.,Burlington, USA)

C. Simulation

Figure 3.a presents the result of the topology optimization algorithm for the case of no polarization optimization. Different interesting features can be observed from the result. The overall structure can be described as a three members lattice beams. When powered, the beams straighten, increasing their length. As a result, two amplification effects are created. First, the beams on the left side play a role of mechanical lever that converts the symmetrical extension to horizontal movement by rotation. Second, the right beam acts as an extension rod, further pushing forward. This behavior is well illustrated in Fig. 3.b presenting the displacement field norm superposed with the resulting deformed domain. We conclude on the raw design by evaluating its performance as presented in section II-D.2 resulting in an amplification ratio of 1.38, i.e. we have obtained +38% displacement than a full plate of the same dimensions.

The result of the topology optimization with polarization optimization is shown in Fig 4 .a. In this figure, the working principle of the grey areas which have the same polarization direction is similar to the one of previous result without polarization optimization. However, the major difference is the part with inverse polarization direction. With the same applied voltage, this part will be retracted and the combination of retraction and extension in the design will result to considerable improvement of displacement amplification. For example, the simulation in Fig 4 .b shows that displacement of the optimized design with optimization of the polarity is almost twice the displacement of the one without the optimization of the polarity in Fig. 3.

The numerical data of simulation is reported in Table I. In this table the displacement per volt, the blocking force, energy density and the related gains for an applied voltage of 5 volt are reported. Energy density is the mechanical energy divided by piezoelectric mass. The mechanical energy is calculated as presented in [START_REF] Wood | Optimal energy density piezoelectric bending actuators[END_REF]. As can be seen in Table I, the displacement gains with respect to full plate for optimized design without polarization and with polarization optimization are 1.42 and 2.82 respectively. On the other hand, the gains of blocking force are decreasing by factors of 0.66 and 0.33 respectively. This is expected since increasing the displacement is always accompanied by reduction of output force. However, it should be noted here that optimized designs have only 30 percent of the volume of the full plate. Therefore, a great improvement in energy density can be seen for optimized designs as it is reported in Table I.

IV. E

A. Prototype fabrication

The fabrication process of the prototypes start by converting the optimzed designs to an exploitable CAD model. Based on the CAD model, a Starlight laser engraving machine from SiroLasertec (Siro Lasertec GmbH, Pforzheim, Germany) is used to achieve cutting operations. Starting from a piezoelectric full plate with a dimensions of 70x70x0.127 mm (from Piezo Systems), the three prototypes shown in Fig. 5(a,b,c) were fabricated: (a) full plate, (b) prototype without optimized polarity and (c) prototype with optimized polarity. For the last one, two cutting operations were necessary. The first one consisted to cut the main structure of the actuator while the second consisted to separate the electrodes following the optimized electrode profile. It is worth to note that a rectangle domain of 5 mm was added on the resulting design. This extension allows to fix the design on a custom made PCB to ensure mechanical clamping. Furthermore, silver-based conductive glue was used to ease the electrical connection with the lower and upper electrodes.

B. Experimental Bench

Figure 5.d shows the experimental set-up used to characterize the prototype. It is composed of:

-an interferometer (SIOS, GmbH company) able to measure displacements (deformations) with nanometric resolution, -the prototypes to be characterized and fixed on PCB supports, -an acquisition system (computer with Matlab + dSPACE board) to generate excitation signals, and to record, process and analyze measured experimental data. The prototypes and the interferometer probe were placed on an anti-vibration table to avoid vibrations that may come from the ground.

C. Characterization of Results and Discussion

The objective of experimental investigation is to measure the displacements in order to demonstrate the efficiency of optimized designs in comparison to full plate and to validate the results obtained in the simulation part. We use sine voltage signals with low amplitude in order to respect linearity assumption of the modelling. A series of measurements is performed under an excitation that goes from 0V to 5V with a step of 0.5V and a frequency of 2Hz. This frequency is low enough in comparison with the plate and prototypes natural frequency such that the static model of the system which was presented in section II-C is respected. The domain of applied voltage is chosen to keep the material behavior only in its linear domain, as the model used for the optimization is based on linear assumptions. Higher voltage and higher frequencies can result in nonlinear behavior as hysteresis [START_REF] Rakotondrabe | Multivariable classical prandtl-ishlinskii hysteresis modeling and compensation and sensorless control of a nonlinear 2-dof piezoactuator[END_REF]. Each series is applied alternatively on the full plate, prototype without polarization and prototype with polarization. The resulting average displacements are reported in Fig. 6 and First, based on Fig. 6, it is obvious that there is a satisfying agreement between the experimental and the simulation data. However by increasing the voltage, the difference between the experimental and the simulation data increases. This is due to the inherit nonlinear behavior of piezoelectric materials. Though, within the domain of the applied voltage, linear assumption in the modelling produces negligible error.

Then, Fig. 6 demonstrates better performances of the optimized designs with respect to full plate design. The improvement factor from the optimized design without polarization optimization is 1.38. On the other hand, the improvement factor from the optimized design with polarization optimization is close to 2.8. Hence, the optimized polarization design proposes 2 times more displacement than the nonoptimized polarization design. This improvement primarily proves that reducing the material volume and optimizing the topology of the structures within the structures domain can provide more flexible design in terms of producing the desired displacement. Furthermore, it is illustrated that polarization optimization allows both expansion and contraction within the design domain and thus amplifies the displacement amplitude significantly.

The results in Fig. 6 and Table I demonstrate two accom-plishments in term of precision and displacement. However, it is true that the displacement amplification is acquired by decreasing the output force of the actuator. Although the experimental measurements of the piezoelectric actuator's output force is highly challenging due to the low thickness to length ratio of the piezoelectric plate, the decrease in the output force is calculated by FEM simulation and it is reported in Table I. In fact, by increasing the stiffness of the modeled spring at the tip of the actuator, the optimization will produce a simple bar beam to maximize the force. On the other hand, the in-plane piezoelectric force is much higher than the out of plane force. Therefore, the output force of the optimized designs are still enough for the applications like atomic microscopy or in applications where the actuator stroke is essential. In summary, the method allows to reduce drastically the material amount. Indeed, only 30% of the material was optimally distributed in order to provide a displacement greater than the displacement of an actuator with a uniform polarization. This leads to a compact and economical design. This is particularly interesting in the context of miniaturization, since the non-occupied space can be utilized to implement additional functionalities like sensors or electronic circuits.

V. C

In this paper, the optimal design of a piezoelectric actuators was presented. Based on topology optimization approach, the design improved the actuator maximum stroke comparing to conventional actuators, owing to the combination of material density and electrodes polarity optimization. To do so, first a finite element model of piezoelectric plate with the help of plane-stress assumption was established and a normalization was introduced to increase the generality of the proposed method. In addition, the polarization was considered as a variable of optimization to give further degree of freedom to the topology optimization algorithm for finding the optimal layout of the piezoelectric actuator. The objective function was defined to maximize the stroke of the actuator based on the defined stiffness of the movable object. To better investigate the effect of polarization, two optimized design with and without polarization optimization were obtained. The performance of these two optimized designs were investigated by numerical simulation and experimental validation. Satisfying agreement between the numerical and experimental results were observed. It was proved that optimization of polarization brings both expansion and contraction within the design domain which will magnify the amplitude of the displacement significantly. The stroke obtained from the optimized design with the polarization optimization is twice larger than that from the design without polarization optimization.

Future work extend the approach to take into account the out-of-plane displacement. Indeed, a large class of piezoelectric actuators use multilayered structure which benefit from large bending displacement in [START_REF]Smart materials-based actuators at the micro/nano-scale. Characterization, Control and Applications[END_REF]. This extend will pave the way for the design of multi-degrees of freedom actuators.

Fig. 1 .

 1 Fig. 1. Piezoelectric material sandwiched between electrodes.

Fig. 2 .

 2 Fig. 2. Pusher geometric and mechanical specification. Gray area represents the design domain

Fig. 3 .

 3 Fig. 3. Optimized design without polarization optimization. a) The obtained layout by optimization b) COMSOL simulated deformation after applying 5(V) actuation voltage.

Fig. 4 .

 4 Fig. 4. Optimized design with polarization optimization. a) The obtained layout by optimization b) COMSOL simulated deformation after applying 5(V) actuation voltage.

Fig. 5 .

 5 Fig. 5. Fabricated prototypes and experimental bench. a) full plate actuator, b) actuator prototype without optimized polarity, c) actuator prototype with optimized polarity following the profile of Fig. 4 (a), d) general view of the experimental set-up.

Fig. 6 .

 6 Fig. 6. Experimental characterization of prototypes.

TABLE I S

 I Simulation ( Input voltage = 5V ) Full plate Opt without pol Opt with pol

	Displacement (nm/V)	57	81	161
	Displacement gain w.r.t.f.p	-	1.42	2.82
	Blocking force (N)	2.56	0.21	0.18
	Blocking force gain w.r.t.f.p	-	0.08	0.07
	Energy density ( / 3 )	4.55	1.81	3.10
	Energy density gain w.r.t.f.p	-	0.39	0.68
		Experiment ( Input voltage = 5V )
		Full plate Opt without pol Opt with pol
	Displacement (nm/V)	62	86	174
	Displacement gain w.r.t.f.p	-	1.38	2.8
	* w.r.t.f.p : with respect to full plate

Table I

 I summarizes the obtained performance gains for both simulations and experimental validation.

This work was supported by the national CODE-TRACK project (ANR-17-CE05-0014-01, Control theory tools for optimal design of piezoelectric energy harvesters devoted to birds tracking devices) and was partially supported by the Bourgogne Franche-Comté region COMPACT project.