E Prost

V Loriot

E Constant

I Compagnon

L Bergé

F Lépine

S Skupin

Air-photonics terahertz platform with versatile micro-controller based interface and data acquisition

We present a recently developed terahertz platform. An air plasma produced by an ultrashort two-color laser pulse serves as a broadband terahertz source, with electric field that has peak amplitude in the MV/cm range. Air biased coherent detection of the terahertz field is employed where a peak detector associated with a micro-controller board acquires the signal coming from an avalanche photodiode. The temporal and spectral profiles of the produced terahertz electric field are presented, in excellent agreement with numerical simulations of the whole setup. We illustrate the capabilities of this platform by performing spectroscopy on water vapor and of a polystyrene reference sample.

I. INTRODUCTION

The terahertz (THz) part of the electromagnetic spectrum, usually defined between 0.1 THz and 30 THz, has attracted great interest in various research fields in the last decades. This recent progress in THz science has led to the advent of numerous applications [START_REF] Tonouchi | Cutting-edge terahertz technology[END_REF][START_REF] Redo-Sanchez | Review of terahertz technology readiness assessment and applications[END_REF] . Terahertz time domain spectroscopy [START_REF] Jepsen | Terahertz spectroscopy and imaging -modern techniques and applications[END_REF][START_REF] Xie | The application of terahertz spectroscopy to protein detection: A review[END_REF][START_REF] Neu | Tutorial: An introduction to terahertz time domain spectroscopy (thz-tds)[END_REF] (THz-TDS) is one of those and it is extensively used in different domains of research ranging from material science 6 to chemistry [START_REF] Fischer | Chemical recognition in terahertz time-domain spectroscopy and imaging[END_REF][START_REF] Bawuah | Advances in terahertz time-domain spectroscopy of pharmaceutical solids: A review[END_REF] . Nowadays, commercial THz-TDS systems are available, making this technique a widespread tool for numerous applications (see [START_REF] Neu | Tutorial: An introduction to terahertz time domain spectroscopy (thz-tds)[END_REF] and references therein). One of the usual limits of THz-TDS systems is the bandwidth of the THz emitters and detectors as well as the amplitude of the produced fields. Usually, emitters and detectors rely on photoconductive antennas or optical rectification in nonlinear crystals with bandwidths limited to 10 THz and suffering from bandwidth gaps due to absorption by optical phonons in the material.

The improvement of THz-TDS systems goes hand in hand with the development of appropriate THz sources [START_REF] Hafez | Intense terahertz radiation and their applications[END_REF][START_REF] Zhang | Intense terahertz radiation: generation and application[END_REF] . An interesting alternative source relies on the use of an air plasma generated by a two color laser pulse [START_REF] Cook | Intense terahertz pulses by four-wave rectification in air[END_REF][START_REF] Kress | Terahertz-pulse generation by photoionization of air with laser pulses composed of both fundamental and second-harmonic waves[END_REF][START_REF] Kim | Terahertz emission from ultrafast ionizing air in symmetry-broken laser fields[END_REF] . Such THz source can be conveniently coupled with air biased coherent detection [START_REF] Karpowicz | Coherent heterodyne time-domain spectrometry covering the entire "terahertz gap[END_REF][START_REF] Ho | Design and performance of reflective terahertz air-biased-coherent-detection for time-domain spectroscopy[END_REF] (ABCD). Using air as the generation and detection medium has several advantages, most notably the ability to produce broadband high amplitude THz fields without any risk of damaging the source or detector. Moreover, the propagation and detection of these broadband THz spectra was demonstrated, with bandwidths going up to 75 THz [START_REF] Kim | Coherent control of terahertz supercontinuum generation in ultrafast laser-gas interactions[END_REF] and beyond [START_REF] Matsubara | Ultrabroadband coherent electric field from far infrared to 200 thz using air plasma induced by 10 fs pulses[END_REF] . Not relying on solid material also means avoiding the appearance of gaps in the THz spectra due to optical phonons.

Here we report on the development of a THz platform based on the above mentioned air plasma source a) Electronic mail: emilien.prost@univ-lyon1.fr and ABCD setup. In contrast to previous works using lock-in amplifiers in the detection process, our data acquisition scheme relies on the use of a peak detector to collect the data and a micro-controller board as an analog to digital converter. The first part of this paper is dedicated to the presentation of the layout of the platform. After describing the optical setup and the signal detection, we introduce the micro-controller based data acquisition and interface of the platform. In the second part of the paper we present experimental results measured with the platform. We begin by showing the temporal and spectral envelope of the produced THz pulses. Numerical simulations of the platform including the THz source and the detection support our measurements and reveal important differences between generated and measured THz spectra. Then we demonstrate the spectroscopic capabilities of the platform. First by measuring the fine structure imprinted by water vapor on the spectrum. Secondly, by using a polystyrene reference sample and comparing the results to measurements obtained with a Fourier transform infrared (FTIR) spectrometer.

II. THZ PLATFORM SETUP

A. Optical setup

Our optical setup relies on THz generation by a twocolor ultra-short laser pulse [START_REF] Cook | Intense terahertz pulses by four-wave rectification in air[END_REF][START_REF] Kress | Terahertz-pulse generation by photoionization of air with laser pulses composed of both fundamental and second-harmonic waves[END_REF][START_REF] Kim | Terahertz emission from ultrafast ionizing air in symmetry-broken laser fields[END_REF] composed of the fundamental harmonic (FH) and its second harmonic (SH).

The key process driving the THz generation has been identified as an asymmetric current caused by photoionization and subsequent acceleration of free electrons [START_REF] Kim | Coherent control of terahertz supercontinuum generation in ultrafast laser-gas interactions[END_REF][START_REF] Babushkin | Ultrafast spatiotemporal dynamics of terahertz generation by ionizing two-color femtosecond pulses in gases[END_REF] . The produced THz field is measured with an ABCD setup [START_REF] Karpowicz | Coherent heterodyne time-domain spectrometry covering the entire "terahertz gap[END_REF][START_REF] Dai | Detection of broadband terahertz waves with a laser-induced plasma in gases[END_REF][START_REF] Lu | Broadband terahertz detection with selected gases[END_REF] where the THz pulse and a synchronized FH laser pulse are focused and superposed in the presence of a DC bias field. Four wave mixing due to the χ [START_REF] Jepsen | Terahertz spectroscopy and imaging -modern techniques and applications[END_REF] non-linearity of air produces a SH radiation, a process also called terahertz field induced second harmonic generation 21 (TFISH). In order to be able to retrieve the THz electric field from the SH pulse, an additional high voltage DC bias field is applied during the four wave mixing process. The produced SH signal and its treatment are discussed hereafter. The optical setup is presented in Fig. 1. It uses 20% of the output of a Ti:sapphire chirp-pulse amplifier (Coherent Legend Elite Duo HP USX) producing 30 fs pulses at a repetition rate of 5 kHz with 2 mJ per pulse at a central wavelength of 800 nm. The 2 W incoming beam is split equally between two arms, the first one being dedicated to THz generation and the second one to THz detection. The full THz beam-path is placed in a dry air enclosure (see gray area in Fig. 1), in order to prevent the strong absorption bands of water in the THz range.

The SH component of our two-color laser pulses is produced in a 100 µm thick β-Barium Borate (BBO) crystal adjusted for type I phase matching and a subsequent ultra-thin half-waveplate at 800 nm (WPD03-H800-F400-SP from Newlight Photonics). The second harmonic generation (SHG) process converts part of the fundamental pulse at the frequency ω into its SH at frequency 2ω. Optimized type I SHG produces a SH field with its polarization axis orthogonal to the FH field. The THz yield is much larger when both components have parallel polarization [START_REF] Kosareva | Polarization control of terahertz radiation from two-color femtosecond gas breakdown plasma[END_REF] , therefore we use a half-waveplate to align both the FH and SH polarization axis. The use of an ultra-thin waveplate is mandatory to keep some temporal overlap between the two colors. The BBO crystal and the waveplate are placed 70 mm before the geometrical focus of a 300 mm lens. Another half-waveplate is placed before the lens to adjust the polarization axis of the FH field such that the polarization axis of the produced linearly polarized THz field and the FH probe pulse are parallel in the detection area.

A few millimeters before the geometrical focus of the lens a plasma is generated by the two color laser pulse and serves as an elongated THz point source. The emitted THz field and the two color laser pulse are collected and collimated by an off axis parabolic mirror. After this first mirror the two color optical pulse is removed by a silicon wafer that transmits about 50% of the THz power [START_REF] Kaltenecker | Attenuation of thz beams: A "how to" tutorial[END_REF] . The transmitted THz field is focused using a second off axis parabolic mirror, before being collimated by a third one. The samples are placed at that intermediate THz focus. A motorized sample holder allows us to perform measurements of several samples without opening the dry air enclosure. The probe part of the FH laser field propagating through a delay line (DL325 from Newport) is dedicated to THz detection. The THz pulse is focused again using the last off axis parabolic mirror. A hole is this mirror allows the FH probe pulse, at frequency ω, to be focused by a 300 mm lens at the same position as the THz pulse. At focus, a static 2 kV bias voltage is applied along the probe and THz polarization direction by a pair of 2 mm long electrodes separated by 2 mm. This results in TFISH through a four wave mixing process between the probe FH field E ω and the low frequency electric field E Ω = E THz + E bias where E THz and E bias are the THz and bias fields, respectively. We use a bandpass filter (FB400-40 from Thorlabs) centered at 400 nm to isolate the produced SH signal, which is then detected by an avalanche photodiode (APD440A2 from Thorlabs). The probe and low frequency fields are both time dependent and the resulting SH signal depends on their crosscorrelation. The SH signal S 2ω received by the photodiode, for a given time delay τ between the probe and THz field, can be expressed as 14

S 2ω (τ ) ∝ χ (3) I ω 2 E 2 THz (τ ) + E 2 bias + 2E bias E THz (τ ) ,
where χ [START_REF] Jepsen | Terahertz spectroscopy and imaging -modern techniques and applications[END_REF] denotes the effective third-order susceptibility of air and I ω the intensity of the probe pulse. Usually, THz-TDS relies on the use of an external low-frequency modulation of E bias and lock-in amplification to extract the THz electric field from the SH signal. Here, we choose a different approach that requires to perform two successive measurements, S + 2ω and S - 2ω , with bias field of opposite sign, E bias = ±10 kV/cm, as shown in Fig. 2(a). By subtracting the two signals S + 2ω and S - 2ω we obtain a signal S THz directly proportional to the THz electric field,

S THz (τ ) = S + 2ω (τ )-S - 2ω (τ ) ∝ 4 χ (3) I ω 2 |E bias | E THz (τ ).
Note that by comparing the magnitude of SH signals with and without bias field, the THz field amplitude can be estimated when E bias is known [START_REF] Iwaszczuk | Terahertz field enhancement to the mv/cm regime in a tapered parallel plate waveguide[END_REF] .

B. Acquisition and interface electronics

A schematic diagram of the interface and data acquisitions used in our THz platform is presented in Fig. 3. The core of the electronics is a micro-controller board (Arduino UNO Rev3) that serves as an interface between various elements of the platform and the computer. We use the LabVIEW Hobbyist Toolkit library 25 to control and interact with the micro-controller. The only other element directly controlled by the computer is the delay line, all other communications are performed via the micro-controller. The desired sign of the applied bias field is sent to the high voltage supply. Motorized beamblocks, composed of a piece of anodized aluminum fixed on a servomotor (Parallax Standard Servo #900-00005), are placed on both pump and probe arms. Their positions are controlled by pulse-width modulation and can thus be controlled by the micro-controller. The Arduino board does not only control all other elements but also serves as an analog to digital converter. We use it to record the experiment signal as well as the humidity inside the enclosure with a sensor (Honeywell HIH-4000-001) placed near the sample holder. The measured signal does not come directly from the photodiode but from a homemade peak detector [START_REF] Achtenberg | Review of peak signal detection methods in nanosecond pulses monitoring[END_REF] , which is used to convert the photodiode signal into a DC voltage readable by the micro-controller board. At its output, the avalanche photodiode delivers a voltage pulse, whose amplitude is proportional to the power of the received light pulse, with a maximum value of 2V. The peak detector takes the highest value of this input electronic signal and stores it as a constant output voltage. Before being processed, this output voltage is amplified by a factor 2.5 to use the whole dynamic range of the analog to digital converter of the micro-controller. A reset of the peak detector is performed before each measurement.

-0.3 -0.2 -0.1 0 0.1 0.2 0.3 0 0.5 1 (a) -0.3 -0.2 -0.1 0 0.1 0.2 0.3 -0.2

III. THZ PLATFORM RESULTS

A. THz pulses produced by the platform

Figure 4 shows experimentally measured temporal and spectral profiles of the THz pulses produced by the platform. For those measurements, no sample was placed in the path of the THz field and the relative humidity inside the enclosure was kept around 9% by injecting dry air. The measured temporal profile is short, with an envelope around 50 fs, and seems to be an almost single cycle pulse. Using measurements taken with a known applied bias voltage of 10 kV/cm and without any bias [START_REF] Iwaszczuk | Terahertz field enhancement to the mv/cm regime in a tapered parallel plate waveguide[END_REF] , the THz electric field strength has been evaluated in the MV/cm range. The THz spectrum is obtained by taking the Fourier transform of the measured temporal electric field. Thus, the chosen temporal step and span for the acquisition also defines the spectral ones. For the measurement presented in Fig. 4, the temporal step is set to 2 fs and the data are taken over a 2 ps span giving a 0.5 THz spectral step and a 250 THz span. These settings are well suited to measure the spectral envelope of the THz field and the shape of the temporal THz field. The spectral intensity profile shown in Fig. 4(b) has a smooth profile extending from 1 THz to 35 THz. This measurement was performed by averaging over 10 laser shots for each data point leading to a total acquisition time of a couple of minutes.

B. Numerical simulations

To understand the different stages of the THz platform and the results presented in the previous section, we also performed numerical simulations in the exact same experimental conditions. For the description of the air-plasma based THz source, we use direct fully space and time resolved computations based on a vectorial version of the unidirectional pulse propagation equation (UPPE) [START_REF] Kolesik | Nonlinear optical pulse propagation simulation: From Maxwell's to unidirectional equations[END_REF] . This numerical model governs the forward- propagating transverse electric field components E x , E y of elliptically-polarized pulses, which are subject to linear dispersion and diffraction; together with third-order nonlinear polarization, photo-ionization and related losses [START_REF] Bergé | Optical ultrashort filaments in weakly-ionized, optically-transparent media[END_REF] . More detail about this vectorial model can be found in [START_REF] Tailliez | Terahertz pulse generation by twocolor laser fields with circular polarization[END_REF] while the physical parameters for air are taken from [START_REF] Nguyen | Spectral dynamics of thz pulses generated by two-color laser filaments in air: the role of kerr nonlinearities and pump wavelength[END_REF] .

In the simulations presented thereafter, we used our experimental laser parameters, namely an incident linearly polarized 115 µJ laser pulse at 800 nm with 30 fs FWHM duration and 6 mm FWHM beam width on the focusing lens L. This fundamental harmonic (FH) pulse is then frequency doubled in the BBO crystal to create a co-propagating second harmonic (SH) pulse (2.7% energy conversion). The polarizations of the FH and SH pulses are aligned by means of a half-waveplate before the pump intensity reaches values above the threshold for ionization and a two-color air-plasma is generated about 2 mm before the geometrical focus, here assumed at z = 0. The simulated peak electron density as a function of the propagation distance z is shown in Fig. 5(a), together with the energy build-up of the generated spectral components below 90 THz. Figure 5(b) reveals that the generated THz spectrum broadens significantly during the generation process, with higher frequencies generated later, that is, closer to the end of the plasma.

It turns out to be crucial to take into account the dispersion properties of the BBO, the quartz half-waveplate as well as the ambient air in order to get the correct beam and pulse profile of the SH component as well as its (slightly elliptic) polarization state at focus. Figure 5(c)-(f) shows the FH and SH field components at the beginning (z = -2 mm) and at the end (z = 0) of the interaction region. In particular the temporal walk-off between the two components is clearly visible, as the SH pulse is delayed by about 25 fs compared to the FH pulse. While at z = -2 mm some overlap near zero delay still enables the photocurrent mechanism for THz generation, plasma defocusing in the trailing part of the FH pulse has completely separated the two colors in the spatio-temporal domain. This is why THz generation stops roughly in the middle of the plasma at z = -1 mm, and the expected THz energy yield remains in the order of few nJ.

Finally, we want to shed some light on the THz detection process, that is, the part of the generated THz field that we can expect to detect in our measurements. One could incorrectly expect to measure a pattern close to the final spectrum at the end of the plasma (brown curve z=0 mm) in Fig. 5(b). This is not the case because those spectra are integrated over the whole transverse plane. In the experimental detection setup only the part of the THz field overlapping with the focused IR probe pulse and the high voltage field is detected. Therefore, it is instructive to compute the THz field in the detection area. The four parabolic mirrors are essentially forming two consecutive 4f systems and thus an image of the THz source is produced. Figure 6(a) shows the evolution of the THz fluence as a function of the propagation distance in the detection area (accounting for 50% transmittance of the silicon wafer). For convenience, the position of the image of the geometrical focus of the laser pump is set to z = 0. The highest THz fluency (and intensity) is obtained at z = -1.5 mm, which corresponds to the position of maximum THz generation of the source. It is reasonable to focus the IR probe pulse at this position in order to maximize the ABCD signal. Consistently with our measurements the IR focal FWHM spot size can be estimated to less than 20 µm, that is, much smaller than the THz spot-size, and thus essentially the on-axis THz intensity is detected.

The simulated on-axis THz spectral intensity at z = -1.5 mm in the detection area is shown in Fig. 6(b). Compared to the spectrum emitted by the source shown in Fig. 5(b) lower THz frequencies are strongly suppressed. The reason for the suppression of lower THz frequencies in the ABCD scheme is that the THz spot size is strongly frequency dependent, as it is limited by the wavelength. Thus, for those THz frequencies the on-axis intensity decreases. Moreover, beam divergence increases with decreasing frequency, and parts of the low frequency field may not be collected by the parabolic mirror.

It is important to note that the air-plasma based THz source is very sensitive to the two-color pump waveform. This is demonstrated by adding small chirp to the incoming FH pulse. This chirp is imprinted to the SH pulse in the generation process in the BBO, and as a result both frequency and phase relation of FH and SH in the region of temporal overlap in the plasma is affected. Even a small chirp corresponding to a group delay dispersion of -100 fs 2 will alter the recorded THz spectrum significantly, as shown by the dotted curve in Fig. 6(b). Both amplitude and width of the spectrum are strongly decreased, rendering the simulated spectrum much closer to the experimental one presented in Fig. 4 of this order of magnitude is likely to be present in the upcoming experiments. Despite this, our measurements are reproducible, thus showing the stability of the platform.

(b). A chirp

C. Spectroscopic features of the THz pulses

Let us now use the platform to investigate the spectroscopic features imprinted in the THz pulses through its beam path. For these purpose, we increase the spectral resolution for a better identification of sharp spectral lines that may appear in the spectra. To this end, we keep the temporal step of 2 fs in our acquisition (giving a spectral span of 250 THz), but we extend the temporal span to 10 ps, thus reducing the spectral resolution down to 0.1 THz. The acquisition time using these parameters and averaging over 10 shots for each data point is below 10 minutes. Resulting temporal and spectral profiles (still without any sample at the THz focus) are presented in Fig. 7. Figure 7 (a) and (b) show different parts of the same signal with adjusted vertical scales. It can be seen that after the main pulse (a) [identical to the pulse shown in Fig. 4(a)], a weak signal is present in the trailing part of the pulse (b). In our setup, this signal is mainly attributed to the presence of ambient water vapor [START_REF] Babilotte | Observation of the field-free orientation of a symmetric-top molecule by terahertz laser pulses at high temperature[END_REF] in the dry air enclosure that is kept at 9% relative humidity. In the spectral domain, it corresponds to absorption bands in the spectrum, as can be seen in Fig. 7(c). Even if the amplitude of this signal is significantly smaller than the main peak, it carries the essential spectral information on the water absorption bands. This amplitude ratio highlights the need for a high signal to noise ratio, which has been evaluated to about 100 for our acquisition system. agreement, giving a good standard to validate our measurements. Another weak absorption band in the measured spectrum, close to 18.1 THz, can be attributed to a phonon resonance in the silicon wafer [START_REF] Kaltenecker | Attenuation of thz beams: A "how to" tutorial[END_REF] .

D. THz time-domain spectroscopy of polystyrene

We finally perform time-domain spectroscopy of a polystyrene sheet and compare it with data extracted from a Fourier Transform Infra-Red (FTIR) spectroscopy database. Such sheets are frequently used as a reference to calibrate FTIR spectrometers in both line positions and amplitudes. The acquisitions parameters are the same as in the previous section, that is, the spectral resolution is set to 0.1 THz. Figure 8(a) shows the measured temporal profiles. The signal with sample is delayed due to global change of the group velocity inside the polystyrene sheet. The polystyrene dispersion stretches the signal significantly, leading to some oscillations appearing between the main pulse and about 0.6 ps.

Absorption leads to the appearance of spectral bands between 10 THz and 30 THz, as well as the water lines discussed previously.

To obtain a better visualization of the polystyrene spectral bands we calculated the transmission spectrum of the sample by dividing the sample spectrum by the reference spectrum. The obtained transmission spectrum is shown in Fig. 8(c), and compared with a reference spectrum from the OPUS software Bruker ATR-Polymer-Library 33 . The main line positions between 15 THz and 30 THz are in good agreement. For lower frequencies there are no absorption lines to compare, and for higher frequencies the signal-to-noise ratio becomes too low to perform any comparison. The signal near the 21 THz line drops to the noise level [see Fig. 8(b)], thus the transmission value calculated around this frequency might not be relevant. It may explain why the measured line shape does not match the one observed in the FTIR spectrum.

IV. CONCLUSION

We have presented an air photonics based terahertz platform with micro-controller interface for basic operations and data acquisition. The platform is thus easy to use and customizable for future experiments. Our detection setup does not require a lock-in amplifier, but only relies on a peak detector and the application of a high voltage bias of switchable sign. Using this simple setup we report a signal to noise ratio above 100. The platform produces short THz pulses with a broadband spectrum reaching 30 THz and a peak electric field strength in the MV/cm range. This makes for a versatile platform usable for THz-TDS but also for other type of applications employing THz pulses as pump or probe such as THz streaking [START_REF] Ardana-Lamas | Temporal characterization of individual harmonics of an attosecond pulse train by thz streaking[END_REF] , molecular orientation [START_REF] Xu | Long-lasting molecular orientation induced by a single terahertz pulse[END_REF] or nonlinear spectroscopy [START_REF] Tcypkin | Giant third-order nonlinear response of liquids at terahertz frequencies[END_REF] .

Comparison with comprehensive numerical simulations of both the THz source and detection reveals the origin of the characteristic spectral shape obtained with ABCD measurements. We also demonstrated the importance of taking into account pump propagation effects before the plasma for THz generation is formed, in particular in the BBO crystal and waveplate. Such detailed confrontation of theoretical results and measurements performed on a given setup is unprecedented. We thus expect strong synergistic effects between experiments and numerical simulations in the future improvement and exploitation of our THz platform.

FIG. 1 .

 1 FIG. 1. Scheme of the experimental setup. BS : Beam splitter (50%); L : focusing lens (f = 300 mm); DL : Delay line; HWP : Half waveplate; BH : composed of a BBO crystal and an half-waveplate; BD : Beam dump; SiW : Silicon wafer; S : Sample; HV : High voltage electrodes; F : 400 nm band-pass filter; APD : Avalanche photodiode. The gray area highlights the part of the setup within an air controlled enclosure.

FIG. 2 .

 2 FIG. 2. (a) Second harmonic signal measured by applying a bias field E bias of 10 kV/cm with positive (red) or negative (blue) sign. (b) Reconstructed terahertz field amplitude from the difference between the two previous measurements. All signals are normalized to the maximum measured signal.

FIG. 3 .

 3 FIG. 3. Scheme of the interface and data acquisition. The arrows represent the direction of the data flow.

FIG. 4 .

 4 FIG. 4. Measured (a) temporal profile and (b) reconstructed spectral intensity of the THz electric field.

  FIG. 5. THz (a) Peak plasma density and THz energy (ν < 90 THz) versus propagation distance relative to the geometrical focus at z = 0. (b) Snapshots of the transverse integrated THz spectrum at several positions z along the plasma. Snapshots of the main FH (c) and SH (d) polarization component of the pump field near the beginning of the plasma at z = -2 mm. (e) and (f) show the same fields near the end of the plasma at the geometrical focus z = 0.

- 1 )FIG. 7 .

 17 FIG. 7. (a) Temporal profile of the main THz pulse [identical to the pulse shown in Fig. 4(a)]. (b) Trailing part of the THz pulse containing the spectroscopic features, with vertical scale adjusted for clarity. (c) Reconstructed spectral intensity (blue curve) of the temporal profile of (a) and (b). The red curve represents the absorption coefficient α extracted from the HI-TRAN database for H2O convoluted by a 0.1 THz FWHM Gaussian. The dashed lines serve as a guide to the eye by indicating most intense line positions. The dotted line indicates the position of a Si absorption band.

Figure 7 (FIG. 8 .

 78 FIG. 8. (a) Temporal profile of the THz electric field with (blue) and without (red) polystyrene sample both normalized to the reference spectrum peak value. (b) Retrieved intensity spectrum with the polystyrene sample. (c) Reconstructed transmission spectrum with polystyrene sample (blue) and polystyrene absorption spectrum extracted from Bruker ATR-Polymer Library (red). The dashed lines serve as a guide to the eye and indicate the positions of the main absorption lines.

ACKNOWLEDGMENTS

We thank B. Moge, I. Aguili, F. Khalid Balyos and C. Clavier for technical support. We thank P. U. Jepsen and B. Zhou from the Technical University of Danemark for fruitful discussions. We acknowledge financial support from CNRS, the ANR ALTESSE2 (ANR-19-ASMA-0007) project, and the Qatar National Research Fund (NPRP 12S-0205-190047). Numerical simulations were performed using resources at Grand Équipement National De Calcul Intensif (GENCI) (A0100507594).

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.


Simon D Carrière
email: simon.carriere@upmc.fr

Bertille Loiseau

Cédric Champollion

Chloé Ollivier

Nicolas K Martin-Stpaul

Nolwenn Lesparre

Albert Olioso

Jacques Hinderer

Damien Jougnot

First evidence of correlation between Evapotranspiration and Gravity at a daily time scale from two vertically spaced superconducting gravimeters

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Evapotranspiration (ET) is a key variable to understanding eco-hydrosystems because it links water, carbon, and energy cycles [START_REF] Fisher | Land-Atmosphere Interactions[END_REF]. ET is a major component of the water cycle, accounting for approximately 70% of precipitated water at the continental scale [START_REF] Oki | Global hydrological cycles and world water resources[END_REF][START_REF] Fisher | The future of evapotranspiration: Global requirements for ecosystem functioning, carbon and climate feedbacks, agricultural management, and water resources[END_REF]. Several approaches exist for estimating ET, but large uncertainties persist and estimates remain difficult to validate, particularly at large scale [START_REF] Verstraeten | Assessment of evapotranspiration and soil moisture content across different scales of observation[END_REF][START_REF] Müller | The use of lysimeters in forest hydrology research in north-east Germany[END_REF]. New approaches for estimating ET are crucially needed to improve our understanding of related processes. Gravimetry is a geophysical method that measures the gravitational acceleration field (g). During the last decades, progress has led to the development of devices that can now achieve precision of a few nm/s² and obtain highly accurate temporal monitoring [START_REF] Van Camp | Hydrogeological investigations at the Membach station, Belgium, and application to correct long periodic gravity variations[END_REF]]. These advances now allow to study the dynamics of ET. Data from the GRACE (Gravity Recovery and Climate Experiment) satellite have been used to study ET [START_REF] Rodell | Basin scale estimates of evapotranspiration using GRACE and other observations[END_REF][START_REF] Long | Uncertainty in evapotranspiration from land surface modeling, remote sensing, and GRACE satellites[END_REF][START_REF] Pan | Detection of human-induced evapotranspiration using GRACE satellite observations in the Haihe River basin of China[END_REF][START_REF] Pascolini-Campbell | A 10 per cent increase in global land evapotranspiration from 2003 to 2019[END_REF]. The last studies have shown a correlation between GRACE data and regional-scale ET estimated by global models [e.g. [START_REF] Zhang | The use of stable isotopes to partition evapotranspiration fluxes into evaporation and transpiration[END_REF][START_REF] Mu | Improvements to a MODIS global terrestrial evapotranspiration algorithm[END_REF]. Superconducting gravimeter (SG) data are commonly used in hydrogravimetric studies that approach ET indirectly via water budgets [e.g. [START_REF] Hector | Hydrogravimetry in West-Africa: first results from the Djougou (Benin) superconducting gravimeter[END_REF][START_REF] Fores | Monitoring saturation changes with ambient seismic noise and gravimetry in a karst environment[END_REF][START_REF] Chaffaut | New insights on water storage dynamics in a mountainous catchment from superconducting gravimetry[END_REF]. Only [START_REF] Van Camp | Direct measurement of evapotranspiration from a forest using a superconducting gravimeter: evapotranspiration and gravity[END_REF] have directly interpreted the gravity signal as soil-toatmosphere water transfers at the daily time step. They found similar dynamics between the decreased signal of an underground SG (i.e., below the root levels) and decreased soil water content measured by TDR (time domain reflectometry) probes at 30, 45, and 60 centimeters. The relationship was observed in an oak and beech forest in Belgium on rain-free days. Direct comparison of ET to the gravity signal was not possible in [START_REF] Van Camp | Direct measurement of evapotranspiration from a forest using a superconducting gravimeter: evapotranspiration and gravity[END_REF]'s study for two reasons: i) the daily water content change measured by TDR probes was not necessarily representative of ET at the plot scale, and ii) the limited signal to noise ratio constrained the authors to stack the signal over 10 days so as to limit perturbations. Indeed, [START_REF] Van Camp | Direct measurement of evapotranspiration from a forest using a superconducting gravimeter: evapotranspiration and gravity[END_REF] discussed gravity signal processing difficulties in removing various perturbations: tidal effects, regional hydrological effects, atmospheric load, instrumental drift or even polar motion. In this study, we broaden the approach: i) by working on the surface-to-depth difference between two vertically spaced superconducting gravimeters, one at the surface and the other at a depth of 512m. Such configuration allows to bypass signal processing issues by eliminating all common large-scale perturbations. ii) We consider local hydrological effects by using measured discharges at the outlet of the karst hydrosystem (Fontaine-de-Vaucluse spring). iii) We explore the direct link between the gravity signal and the daily ET dynamic calculated with the SimpKcET model validated on three (3) flux towers near the site [START_REF] Ollivier | An evapotranspiration model driven by remote sensing data for assessing groundwater resource in karst watershed[END_REF]. For the first time, these conditions make it possible to relate gravity directly to ET signals at the daily time step.

Materials and Methods

Site

The experimental site is located in southern France (43°56'28.7''N; 5°29'01.1''E) within the Fontaine-de-Vaucluse hydrosystem which covers 1200 km². The site is part of the French critical zone observation network (OZCAR -http://www.ozcar-ri.org/). The karst hydrosystem is developed in Cretaceous limestones. The upper formations (≈ 300 m thick; Fig. 1) are lower Aptian / upper Barremian sub-reef limestones [Masse and Fenercci-Masse 2011] with high porosities (5-20% in the study area [START_REF] Léonide | Diagenetic patterns and pore space distribution along a platform to outer-shelf transect (Urgonian limestone, Barremian-Aptian, SE France)[END_REF]). Seasonal variations in water content, measured by MRS (magnetic resonance sounding), revealed the strong capacitive hydrogeological role of this limestone [START_REF] Carrière | The role of porous matrix in water flow regulation within a karst unsaturated zone: an integrated hydrogeophysical approach[END_REF][START_REF] Mazzilli | Surface Nuclear Magnetic Resonance Monitoring Reveals Karst Unsaturated Zone Recharge Dynamics during a Rain Event[END_REF]. The underlying formations are lower Barremian limestones with significantly lower porosities (< 5% [START_REF] Masse | Contribution a l'étude de l'Urgonien (Barrémien -Bédoulien) des Monts de vaucluse et du Luberon[END_REF][START_REF] Guglielmi | Carbonate reservoirs -The Rustrel-Fontaine de Vaucluse LSBB: a unique facility for deep carbonate reservoir[END_REF]). This last formation has mainly a transmissive hydrogeological role. Water discharge at Fontaine-de-Vaucluse spring is representative of the deep flow under the LSBB (Fig. 1C) for several reasons: i) soil and rock located between the two SGs (upper Barremian / lower Aptian limestones) have a strong capacity effect regarding water transfer due to high porosity [START_REF] Carrière | The role of porous matrix in water flow regulation within a karst unsaturated zone: an integrated hydrogeophysical approach[END_REF], ii) rocks located under the deep SG (lower Barremian limestones) have only a transmissive role due to low porosity and high permeability through karst conduits, iii) the Fontaine-de-Vaucluse hydrosystem is a highly developed karst where water transfers within the saturated zone are very rapid (piston effect in a few hours) as demonstrated by [START_REF] Garry | Etude des processus d'écoulements de la zone non saturée pour la modelisation des acquifères karstiques[END_REF] using an Mg 2+ concentration time series analysis and [START_REF] Lesparre | New insights on fractures deformation from tiltmeter data measured inside the Fontaine de Vaucluse karst system[END_REF] using tiltmeter measurements, iv) Fontaine-de-Vaucluse discharge dynamics are similar to piezometric dynamics observed in a borehole near the site LSBB. The topography of the study area is a southward slope. The elevation ranges between 500 and 1100 m.a.s.l. whereas the top of the saturated zone within the karst is at approximately 100 m.a.s.l. (Fig. 1). The thick vadose zone made it possible to excavate a tunnel 3.8 km long with a depth ranging between 0 and 512 m mm deep. This structure was a military bunker converted into a scientific research laboratory named "LSBB" in 1997 (http://lsbb.eu/presentation/). This underground facility cuts through the limestone and provides unique access to the karst unsaturated zone (Fig. 1).

Figure 1. Site presentation:

A/ Location map of the experimental site in France. B/ LSBB experimental site located in the Fontaine-de-Vaucluse catchment with a simplified geological context. C/ Simplified water balance of the system included between the two gravimeters where P is precipitation, ET is evapotranspiration, and Q is deep water discharge. Runoff is commonly neglected due to the very high permeability of the karst. D/ Simplified geological cross-section showing gravimeter locations and the influence cone of 90% of gravity signal.

The climate is Mediterranean with hot dry summers; most rainfall occurs during spring and autumn. Between 2003 and 2015, average annual rainfall was 909 mm and ranged between 407 and 1405 mm. Average annual temperature over the same period was 12.9 °C. Runoff is null at the scale of the massif due to the high karst permeability. No streams are present in the study area. Vegetation cover in the study area is primarily forest, dominated by Holm oaks on the southern slope [START_REF] Carrière | Process-Based Vegetation Models Improve Karst Recharge Simulation Under Mediterranean Forest[END_REF]] and white oaks and Montpelier maples on the plateau.

Superconducting gravimeter

The LSBB is a unique site in the world for hydrogravimetric research because it is the only place to date where two superconducting gravimeters (SG) are installed in vertical alignment. This configuration allows us to study hourly mass changes in the area between the two gravimeters using the surface-to-depth approach developed by [START_REF] Jacob | Timelapse surface to depth gravity measurements on a karst system reveal the dominant role of the epikarst as a water storage entity[END_REF] and [START_REF] Champollion | Estimating epikarst water storage by time-lapse surface-to-depth gravity measurements[END_REF]. The SG iOSG#024 [START_REF] Boy | Superconducting Gravimeter Data from Rustrel (LSBB)-Level 1[END_REF] was installed in 2015 at a depth of 512 m (507 m.a.s.l.) and the SG iGrav#031, installed at 1019 m.a.s.l., has been operational only since 2020. Both devices collect continuous data at the minute time step (iGrav#031 and iOSG#024). The instrumental drift of each gravimeter is corrected using repeated measurements of absolute gravity. In this paper, we only explore the surface-to-depth residual signal obtained by subtracting the deep SG signal from the surface SG signal. Residuals from single gravimeters are presented as supplementary material (S1). This simple difference makes it possible to circumvent the usual processing needed to eliminate all large-scale effects (e.g. tide, polar motion) while obtaining a signal with very low residual noise. This noise cancellation has already been demonstrated on laterally separated gravimeters [START_REF] Kennedy | Direct measurement of subsurface mass change using the variable baseline gravity gradient method[END_REF][START_REF] Goto | Simultaneous gravity measurements using two superconducting gravimeters to observe temporal gravity changes below the nm s-2 level: ocean tide loading differences at different distances from the coast[END_REF]. Compared to single gravity measurements, the sensitivity of differential gravity is modified by: i) removing common largescale effects; ii) enhancing the signal to noise ratio; iii) the sensitive volume is smaller than a single gravimeter (focused on the zone between both SG); iv) the sensitive volume extend laterally. Indeed, this vertical configuration offers a greater sensitivity on the sides of the device compared to a single gravimeter [START_REF] Cooke | Characterisation of a new mobile absolute quantum gravimeter : application in groundwater storage monitoring[END_REF]; Supplementary Information S2). The lateral sensitivity combined with the reduced common noise is a great advantage for studying near surface processes such as evapotranspiration. We calculated the Daily Change in Gravity (DCG) by subtracting the signal average between 0h and 5h on day "d" from the same average on the next day "d+1" (Fig. 2B). Standard deviations were calculated each day over the same hourly period and then summed by pair of days to represent error bars in Figure 3b. During rain-free periods, the DCG consists mostly of signals from the lower (Q) and upper (ET) outlets of the karst system (Fig. 1C). Consequently, we subtract the gravitational attraction caused by changes in groundwater storage due to discharge at Fontaine-de-Vaucluse (Q) from DCG and compare the result to ET. Downward flow Q is expressed in mm by dividing the flow at the outlet of the Fontaine-de-Vaucluse spring by the surface of the entire hydrosystem (1200 km²). Surface-to-depth residuals can be expressed in mm by dividing acceleration (initially in nm/s²) by the sum of absolute Bouguer constant for LSBB site [START_REF] Kumar | Groundwater monitoring and characterization by a vertical dipole of superconducting gravimeters in a karst aquifer[END_REF]. They found -0.34 nm/s².mm -1 for the deep gravimeter and 0.39 nm/s².mm -1 for the surface gravimeter. This double application of Bouguer's constant is equivalent to dividing the gradiometer signal by two in order to cancel the double detection of mass variations between the two gravimeters. As soil water is not spatially and temporally homogeneous, one can expect temporal changes of the Bouguer constant due to the precipitation patterns and the heterogeneity of the vegetation. Estimating a variable Bouguer coefficient in such complex medium as karst cannot be done from gravity data alone and numerical modeling including eco-physiological processes are needed.

Evapotranspiration modeling

The SimpKcET model was developed by [START_REF] Ollivier | An evapotranspiration model driven by remote sensing data for assessing groundwater resource in karst watershed[END_REF] to quantify ET over large karst areas. It was coupled to the KaRaMel model, which simulates water flow within karst systems [START_REF] Ollivier | Karst recharge-discharge semi distributed model to assess spatial variability of flows[END_REF], to improve the simulation of karst aquifer processes. SimpKcET was adapted from the FAO56 model [START_REF] Allen | Crop evapotranspiration-Guidelines for computing crop water requirements-FAO Irrigation and drainage paper 56[END_REF]] by considering parameterizations of soil and plant water transfers adapted to the particularities of karst surface layers (strongly heterogeneous with large amounts of coarse elements). SimpKcET uses the vegetation fraction cover derived from remote sensing data to force the evolution of plant transpiration from reference evapotranspiration inputs and a simple description of daily soil water balance including its impact on vegetation water stress. MODIS data was used to derive vegetation fraction cover at 1 km spatial resolution (calculated from the MODIS EVI). Meteorological inputs, such as precipitation and reference evapotranspiration, were obtained from SAFRAN reanalysis at a spatial resolution of 8 km [START_REF] Vidal | A 50-year highresolution atmospheric reanalysis over France with the Safran system[END_REF]. SimpKcET was favorably evaluated over three (3) flux towers on two (2) typical Mediterranean forested fields and one crop field in southeastern France [START_REF] Ollivier | An evapotranspiration model driven by remote sensing data for assessing groundwater resource in karst watershed[END_REF]. In this paper, we used SimpKcET to model the evapotranspiration that was compared to the gravity data. ET was extracted from the 15 pixels covering the area 2 km around both SGs. These 15 pixels were used to evaluate ET (average) together with the ET variability.

Results

The surface-to-depth gravity residuals reveal a low-noise time series (Fig. 2A). Two phases can be clearly identified on the time series: i) gravity signal growth related to rainfall that increases the system mass between the two SGs, ii) signal decrease related to water loss from the system (ET and Q). Throughout the gravity time series, we notice several spikes, which are not related to hydrological processes but mainly to earthquakes and explosions due to the excavation of a new gallery within the LSBB. During winter and spring, the residual time series is jagged because aquifer recharge occurs during this period. In fact, signal increases are synchronous with rainfall and increased flow at Fontainede-Vaucluse. During the summer, we observe a prolonged gravity signal decrease, which is linked to ET and Q processes. We can observe a consistency between, on the one hand, the slope reduction of the gravity time series and, on the other hand, the decrease of ET and Q related to the beginning of drought. At the daily time step, we can clearly see that the signal decreases in a step-wise pattern due to ET, as described by [START_REF] Van Camp | Direct measurement of evapotranspiration from a forest using a superconducting gravimeter: evapotranspiration and gravity[END_REF] on the gravity signal or as other authors describe on piezometry [e.g. [START_REF] Healy | Using groundwater levels to estimate recharge[END_REF][START_REF] Carlson Mazur | Estimating evapotranspiration and groundwater flow from water-table fluctuations for a general wetland scenario[END_REF] and the soil humidity time series [e.g. [START_REF] Nachabe | Evapotranspiration of two vegetation covers in a shallow water table environment[END_REF][START_REF] Jackisch | Estimates of tree root water uptake from soil moisture profile dynamics[END_REF]]. The gravity signal drops during the day as a result of water loss in the system due to ET whereas the signal is more stable during the night (Fig. 2B). In fall, gravity time series becomes jagged again due to the beginning of seasonal rainfall and a new recharge period.

Interpretation and discussion

The dataset allowed us to directly relate the gravity signal to ET. We exploited the unique configuration of two vertically spaced superconducting gravimeters to explore surface-to-depth gravity differences. We studied the June to September 2020 period in detail as it is a period of increased ET and low rainfall. Under these favorable conditions, we were able to relate a) mass changes measured by gravimeters minus water discharge (Q) measured at the karst system outlet to b) ET modeled using the SimpKcET model [START_REF] Ollivier | An evapotranspiration model driven by remote sensing data for assessing groundwater resource in karst watershed[END_REF].

Figure 3B shows strong Pearson (0.73) and Spearman (0.71) correlations between ET and DCG -Q. These relations are broadly significant (p-values < 4.10 -10 and 9.10 -9 , respectively). The small bias (0.55 mm) highlight a slight underestimate of modeled ET. It can be seen that the Q signal is fairly constant over the summer 2020 period (Fig. 2) because it shapes an extremely low gradual decay. As a result, the DCG signal can be primarily attributed to daily ET. Our approach cannot be used for situations such as: i) rainy periods because water inflow into the system and related fast flows create a signal that dominates ET, ii) atmospheric conditions such as storms and periods of high or low pressure, iii) earthquakes, iv) explosions related to tunnel excavation within the LSBB. Because of these disturbances, 40 days indicated by the red areas (Fig. 3A) were removed from the analysis.

On the zoom of surface-to-depth residuals (Fig. 2B) it can be seen that the gravity signal rises during the afternoon. The rise is also seen in the stacked data of [START_REF] Van Camp | Direct measurement of evapotranspiration from a forest using a superconducting gravimeter: evapotranspiration and gravity[END_REF]. We assume this phenomenon to be related to the change in lower atmosphere conditions (pressure, temperature, humidity) during the afternoon. Indeed, the water evapotranspired during the day is removed by the convective exchange between the boundary layer and the free troposphere and trapped in the atmospheric boundary layer at the end of the day when convection vanished [START_REF] Champollion | GPS water vapour tomography: preliminary results from the ESCOMPTE field experiment[END_REF]. Despite this misunderstood phenomenon, our analysis remains largely unaffected because we calculate DCG over two consecutive days at fixed times.

Conclusions

Determining daily ET by gravimetry is a major challenge that approaches the technological and signal processing limits of current gravimetry. Indeed, the objective is to detect daily variations that are sometimes less than 1 nm/s². To reach conclusive results, we have enumerated the conditions that make it possible to highlight the expected signal: i) we selected days without rain during a summer period (June -September) when the ET signal was quite strong (> 1 mm/j).

ii) we used data from two gravimeters to subtract the signals with two consequences: suppression of undesirable large-scale effects (e.g., tide) and maximization of the gravity signal between the two SGs. Considering these last points, we were able to highlight that, in the studied system (Fig. 1B), the daily water outflow by ET (from SimpKcET model) agreed reasonably well with daily mass change (from residuals gradients) subtracted from the discharge at the outlet of the system (from Fontaine-de-Vaucluse discharge). The water stock dynamic is a key variable in hydrology that is generally difficult to quantify and the surface to depth superconducting gravimetry can quantify this key parameter with a precision around 1 nm/s². The high accuracy allows to indirectly estimate the ET at daily time step via the water stock variation. This work is at the interface between geophysics and hydrology. The results open new perspectives for both communities: i) in gravimetry, which can now consider ET to be a process that can significantly create noise in the gravity signal at the subdaily time step, according to the results of this study ii) in hydrology, which can consider gravimetry as a new method for approaching ET, with four main advantages: a large integration volume, a stable influence area over time (unlike Eddy covariance or scintillometry), and minimal required maintenance. However, the installation of a SG is sometimes more cumbersome than a flux tower (depending on the vegetation cover) and requires a main power supply.

Our results encourage collaborative work among hydrologists and geophysicists to improve signal processing in order to deduce ET from the gravimetric signal at experimental sites equipped with a single SG. To achieve this goal, it may be useful to implement additional measurements to improve signal processing (e.g. full parameter meteorological station, sap flow sensors on trees). We can also expect that in the future decades gravimeters pairs (or even network of gravimeters) will be deployed at multiple experimental sites with one SG at the surface and the other deep inside a borehole. Finally, it will be interesting to make a direct comparison between gravity measurements and ET measurements (e.g. eddy-covariance, scintillometry).

Figure 2 .

 2 Figure 2. A/ Time series during the year 2020 for: precipitation (P); evapotranspiration (ET) calculated with the SimpKcET model [Ollivier et al. 2020] where gray lines are individual pixel time series and the sky blue line is the mean value calculated using 15 pixels within 2 km of the superconducting gravimeters; daily Fontaine-de-Vaucluse (FdV) discharge (Q) calculated in mm/d by dividing discharge (in m 3 /d) by surface of FdV catchment (1200 km²); and temporal dynamics of the surface-to-depth gravimetry. B/ two-day close-up on temporal dynamics of the surface-to-depth residuals to show how daily change in gravimetry (DCG) is determined.

Figure 3 .

 3 Figure 3. Comparison of measured surface-to-depth gravimetry signal minus discharge (Q) measured at the Fontaine-de-Vaucluse spring with the modeled evapotranspiration (ET). A/ Temporal graphs. The red areas are periods excluded from the analysis because they are affected by various phenomena (earthquake, storm, rain and explosion). ET gray lines are individual pixel time series and the sky blue line is the mean value calculated within 15 pixels surrounding the superconducting gravimeters. The cumulated daily flow at Fontaine-de-Vaucluse was divided by the surface of hydrosystem (1200 km) to express discharge in mm/d. P is precipitation. B/ Correlation between daily ET modeled and daily DCG -Q measured (r: Pearson correlation; RMSE: root mean square error; n: number of value).

Acknowledgments, Samples, and Data

 The authors would like to express their gratitude to CIRAME and to all the LSBB team.

Data used in this study were collected by the SNO KARST and H + . SNO KARST and H + are part of the OZCAR Research Infrastructure, which is supported by the French Ministry of Research, French Research Institutions and Universities. SimpKcET was developed within the framework of the PITEAS project and the TRISHNA project, which received financial support from CNES. The superconducting gravimeter iGrav#031 was funded by the Equipex CRITEX project and the superconducting gravimeter iOSG#024 was funded by the Equipex MIGA project.

Open Research

 Superconducting gravimetry data is accessible using the following links: surface SG iGrav#031 (http://igets.u-strasbg.fr/SG/ru031/Level1/) and deep SG iOSG#024 (http://igets.u-strasbg.fr/SG/ru024/Level1/)