Céline Décuq
email: celine.decuq@inrae.fr

Marjolaine Bourdat-Deschamps

Pierre Benoit

Colette Bertrand

Rachid Benabdallah

Baptiste Esnault

Brigitte Durand

Benjamin Loubet

Clémentine Fritsch

Céline Pelosi

Sabrina Gaba

Vincent Bretagnolle Df

Carole Bedos

A multiresidue analytical method on air and rainwater for assessing pesticide atmospheric contamination in untreated areas

Keywords: landscape, concentration, impact, wet deposition, SBSE Graphical Methodology, Conceptualization, Validation Conceptualization, Resources, Supervision, Review Colette Investigation, Conceptualization Project administration, Review, Resources, Funding acquisition Curation, Investigation, Brigitte Durand a : Resources, Data Curation, Investigation, Conceptualization, Investigation, Review, Supervision

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

Pesticides are commonly used in agriculture to protect crops from pests and diseases leading to contamination of the environment. The atmospheric compartment contributes to their dispersal at different distances from the application areas [START_REF] Van Dijk | Atmospheric Dispersion of Current-Use Pesticides: A Review of the Evidence from Monitoring Studies[END_REF]FOCUS 2008).

Consequently, organisms in untreated areas (e.g. [START_REF] Pelosi | Residues of currently used pesticides in soils and earthworms: A silent threat?[END_REF] and also human populations (e.g. Dereumeaux et al., 2020[START_REF] Teysseire | Assessment of residential exposures to agricultural pesticides: A scoping review[END_REF]) can be exposed. Pesticides are emitted into the atmosphere during application by dispersion of spray droplets and after application mainly by volatilization. Once in the atmosphere, compounds can partition into their gaseous and particulate phases, be degraded, transported in the wind, and finally deposited by dry deposition (gaseous or particulate phases deposited by turbulence) or wet deposition (by rain) [START_REF] Woodrow | Pesticides and Related Toxicants in the Atmosphere[END_REF].

Both dry and wet deposition should be considered in risk assessments for non-target ecosystems. Indeed, FOCUS (2008) recommends that during registration procedures dry deposition should be taken into account in risks assessments for non-target ecosystems located close to the treated plots, in addition to the drift already considered. Through modeling studies coupling volatilization and local atmospheric dispersion models, Asman et al. (2001) andBedos et al. (2013) found that the contribution of dry deposition to water surface contamination may be higher than the contribution of drift for volatile pesticides. Regarding wet deposition, [START_REF] Potter | Assessing pesticide wet deposition risk within a small agricultural watershed in the Southeastern Coastal Plain (USA)[END_REF] analyzed three years of rainfall monitoring in a small watershed in the USA and concluded that wet deposition should be included in pesticide risk assessments. Wet deposition rates can reach high levels, depending on the compound. For instance, studies in Europe [START_REF] Van Dijk | Atmospheric Dispersion of Current-Use Pesticides: A Review of the Evidence from Monitoring Studies[END_REF]Asman et al., 2005) found wet deposition rates ranging from a few mg.ha -1 .y -1 to several g.ha -1 .y -1 . More recently in Sweden, Kreuger et al. (2017) carried out long-term monitoring of pesticide concentrations in rainwater and calculated wet deposition rates ranging from 9 to 800 mg.ha -1 .y -1 , which represents up to 0.1 % of the applied dose for some compounds. All these deposition studies relied on pesticide monitoring in the air and rainwater.

Regarding airborne pesticide concentrations, most monitoring studies, which are generally carried out to assess population exposure, focused on background concentrations. Datasets on background concentrations have been available in various regions of France since the early 2000s, thanks to campaigns conducted by the Approved Air Quality Monitoring Associations (AASQAs) (data available in the Phytatmo database 2019; see also Désert et al., 2018 and[START_REF] Villiot | Temporal and seasonal variation of atmospheric concentrations of currently used pesticides in Champagne in the centre of Reims from 2012 to 2015[END_REF]. There are however, no continuous monitoring programs of rainwater contamination organized at the French national level. However, data is available in reviews such as Dubus et al. (2000), Asman et al. (2001) or Bedos et al. (2002), which concluded that the range of concentrations in rainwater was large, from non-detectable to values higher than 100 ng.L -1 . More recently, [START_REF] Scheyer | Pesticides analysed in rainwater in Alsace region (Eastern France): Comparison between urban and rural sites[END_REF] observed high pesticide concentrations in Eastern France, particularly during periods of intensive applications on crops. The quantitative assessment of the contribution of atmospheric contamination and deposition to untreated areas, including exposure of the human population, remains a challenge as agricultural practices are constantly evolving (e.g., banning, introducing new molecules on the market, using biocontrol products). Monitoring programs are therefore still needed to assess the state of environmental contamination and its temporal trends, as well as to evaluate the effectiveness of environmental regulation policies, as stressed by Asman et al., (2005).

Various measurement systems are available to quantify pesticide contamination of the atmosphere. This is a challenging measurement because the concentrations to be detected can be relatively low, less than 2 ng.m -3 (Phyatmo 2019), and complex mixtures of compounds are present. Regarding the gaseous and particulate phase, active and passive samplers are used to trap pesticides from the air [START_REF] Yusà | Sampling and analysis of pesticides in ambient air[END_REF][START_REF] Coscollà | Pesticides and Agricultural Air Quality[END_REF]. Passive samplers involve exposing a sorbing phase to ambient air to collect the gas phase. The sampled pesticides migrate by molecular diffusion to the adsorbent where they are trapped .

The sampling time varies from a few weeks to several weeks. Because an electricity supply is not required, the use of passive samplers over long periods in remote areas remains a convenient technique. The main drawback of this sampling technique is that the volume of air passing through the adsorbent during exposure must also be accurately assessed [START_REF] Wania | Development and Calibration of a Resin-Based Passive Sampling System for Monitoring Persistent Organic Pollutants in the Atmosphere[END_REF]Lévy et al., 2020). According to [START_REF] Coscollà | Pesticides and Agricultural Air Quality[END_REF], although improvements are possible, passive sampling will always remain more uncertain than active sampling. Active samplers actively pump air through an adsorbent to collect the pesticides present in particulate and gaseous phases. The sampling time is shorter than for passive samplers, typically a few hours to a few days, and it can be reduced depending on measurement objectives and the sensitivity of chemical analyses. In contrast to passive samplers, active samplers allow accurate measurement of the volume of air sampled using a flow meter, but require a power supply, which can be difficult in remote areas and this is the main drawback of this method.

Once trapped in the adsorbent, extraction of the sorbing phase is required for both techniques.

This can be carried out either using organic solvents or thermal desorption (TD). Using solvents increases the uncertainties due to the numerous preparation steps required before analysis (extraction, concentration) and it is also time consuming. Thermal desorption does not require a concentration step. Furthermore, when coupled directly with gas chromatography (GC) and mass spectrometry (MS) systems, it has a lower limit of quantification than the solvent extraction technique (Bedos et al., 2006). Thus, active sampling coupled with TD-GC-MS is an interesting analytical technique for quantifying airborne pesticides. (Bedos et al., 2006 ;[START_REF] Clément | Adsorption/thermal desorption-GC/MS for the analysis of pesticides in the atmosphere[END_REF][START_REF] Briand | Assessing the transfer of pesticides to the atmosphere during and after application. Development of a multiresidue method using adsorption on Tenax and thermal desorption-GC/MS[END_REF]. Two drawbacks of this technique are that it can only be used for non-thermosensitive compounds, and the sample is destroyed during the analysis, preventing further analysis.

To quantify pesticide concentrations in rainwater, Stir Bar Sorptive Extraction (SBSE)

techniques (Baltussen et al., 1999) can be used, especially for small sample volumes. A larger amount of sorptive extraction phase is available in comparison to other extraction techniques such as SPME (Solid Phase Microextraction). As a result, more pesticides can be extracted using SBSE, which increases the sensitivity and decreases the limits of quantification of this technique [START_REF] David | Stir bar sorptive extraction for trace analysis[END_REF]. The Stir Bars (Twister®) are directly extracted by TD [START_REF] Sandra | Multi-residue screening of pesticides in vegetables, fruits and baby food by stir bar sorptive extraction thermal desorption capillary gas chromatography mass spectrometry[END_REF] and analyzed by GC-MS. The same analytical pipeline can therefore be used for both gas and rainwater matrices. This is an advantage that facilitates monitoring concentrations in both air and water and reduces measurement uncertainties. The multiresidue method has the advantage of decreasing the number of analyses. However, it is generally a complex process that requires compromises due to the diversity of the chemical properties of the compounds.

This study was part of a project being conducted in the Long Term Socio-Ecological Research (LTSER) (Bretagnolle et al. 2018a), aiming to investigate the contamination of several environmental matrices (air, rainwater and soils) and soil macro-fauna (earthworms, carabids and small mammals) by pesticides in an intensive agricultural landscape [START_REF] Pelosi | Residues of currently used pesticides in soils and earthworms: A silent threat?[END_REF]. Here, we focused on assessing the contribution of atmospheric pollution to this exposure. Among the pesticides frequently used in the sampling area and those showing potential emission to the atmosphere, 27 molecules were targeted. Our objective was to develop a multiresidue analysis method using a single analytical pipeline to quantify both atmospheric and rainwater pesticide concentrations in sample collected in natura.

Our objective was to develop a multiresidue analysis method using a single analytical pipeline to quantify both atmospheric and rainwater pesticide concentrations in sample collected in natura. A comparison of the data thus obtained on rainwater contamination with existing data on air contamination in the same region is presented together with an analysis based on compound physico-chemical properties. Wet deposition (µg.m -2 ) was calculated given the rainfall amount and pesticide quantified concentration to assess its contribution to the contamination of untreated areas.

2 Materials and methods

Study area, rainwater and air sampling

This study was carried out in the Long Term Socio-Ecological Research (LTSER) area, known as the --en-Bois,

France (Bretagnolle et al., 2018a). Biodiversity, ecological functions as well as farmer practices and land use are annually monitored in this intensive cereal plain of 435km². In 2018, agricultural practices were recorded (especially pesticides) during farmer interviews in at least 100 fields per year (see Bretagnolle et al. 2018a, b). Farmer interviews allow very accurate data on pesticide use to be recorded (products, active molecules, timing of application etc.).

Rainwater was collected at the research station « Centre d'Etudes Biologiques de Chizé » on a weakly basis from March to November 2018 with a refrigerated collector at 4°C (Eigenbrodt NSA 181/KS, DE). The samples were stored in polyethylene terephthalate bottles at -20°C before analysis. The amount of rainwater per week was measured with an ARG100 Tipping Bucket rain gauge and was recorded with a CR10 data logger (Campbell Sci, USA). Rainwater sampling was in a forested area and assumed to represent the background level of wet deposition in the whole area.

For air sampling, two 1 km 2 landscapes windows with contrasting agricultural management practices were selected: the first contains a high proportion of conventionally farmed arable crops (98 % hereafter CF landscape ) and the second a high proportion of organically farmed arable crops (69 % hereafter OF landscape ).

Air sampling was conducted daily from the 22 nd May to 1 st June 2018. Ambient air was sampled at 1.6 m above the ground surface at a flow rate ranging from 0.52 to 1.13 L.min -1 over 24 hour sampling periods. Pesticides were trapped on cartridges filled with 60 mg Tenax TA adsorbent (Gerstel®, DE). The deployed sampling system operated on a 12V supply and therefore did not require a nearby power source. Volumetric controllers (Gallus 2000, Schlumberger, FR) were used to measure the volume of air sampled by each cartridge. Before the sampling campaign, cartridges were pre-conditioned by heating at 280°C under a helium stream of 60 L min -1 for 6 hours. After exposure, cartridges were closed and stored in the dark at 4°C until analysis.

Selected pesticides

Herbicides, fungicides, and insecticides were studied. In total, 27 pesticides belonging to 14 chemical substance families were analysed (Table 1). These pesticides were chosen as they were frequently used in the sampling area (see [START_REF] Pelosi | Residues of currently used pesticides in soils and earthworms: A silent threat?[END_REF], were compatible with multiresidue analysis, and have significant potential for atmospheric emission. Some have isomers as specified in Table 1 and a total of 33 molecules was thus considered.

Analytical method development

Rainwater sample extraction

SBSE was used to extract pesticides from rainwater. SBSE is based on the partitioning of pesticides between the water sample studied and a specific polymeric phase supported by a magnetic bar (called a "Twister , Gerstel®10mm×0.5mm with a 24 µL Polydimethylsiloxan (PDMS) coating). Several factors can modify the equilibrium and affect the extraction yield. to be determined (Bourdat-Deschamps et al., 2007) and followed the equation given in Eq.S1.

Before being used, Twister® were preconditioned using Gerstel TCD2 C200. They were heated at 280°C under a helium stream of 60 mL.min -1 for four hours. Up to two Twister® bars could be extracted and analyzed at the same time (in the same analytical tube called a TDU TUBE)

by the analytical pipeline used in the laboratory. It was thus possible to choose up to two different extraction conditions. This procedure helped determine the best compromise for the extraction of all molecules, since one condition was more efficient for some compounds and another for the other compounds. The optimised extraction method applied to rainwater samples was the following: 30% (m/v) NaCl was added to 10 mL of rainwater sample and a first extraction was carried out for 135 min (1st Twister®). 10% (v/v) acetone was added to a second aliquot of 10 mL rainwater sample and a second extraction was carried out for 135 min with a second stirring Twister®. Then both Twister® bars were placed together in a (TDU tube described in the following section) for analysis.

Analysis by TD-GC-MS

TD-GC-MS was used to analyze the pesticides in rainwater and air samples. Tenax tubes and Twister® were desorbed using a Thermo Desorption Unit (TDU, Gerstel, DE). The TDU was programmed to desorb pesticides from tubes and Twister® from 50 to 280 °C for 8 min at a rate of 60 °C.min -1 . Pesticides were then cryo-focused in the programmable temperature vaporization injector at -20 °C using a baffled glass liner. Compound separation was carried out using an Agilent 7890B gas chromatograph equipped with a capillary column (30 m length, 0.25 mm internal diameter, 0.25 µm df, HP-5MSUI column, Agilent). The oven temperature was initially set at 70 ºC hold 2 min, heated at a rate of 25ºC.min -1 to 150 ºC, then heated at a rate of 3°C.min -1 to 200°C and finally at 8°C.min -1 to 300°C hold 1 min. Helium was used as carrier gas, with a flow rate of 1.5 mL.min -1 . As mentioned in Table 1, several of the studied pesticides have isomeric forms: cypermethrin (four isomers), cyproconazole (two isomers), propiconazole (two isomers) and tau-fluvalinate (two isomers). With the chromatographic conditions used, a total of 32 compounds out of the 33 possible compounds were separated: the 27 selected pesticides given in Table 1, three isomeric forms of cypermethrin separated in addition to the main one, one isomeric form of propiconazole and one isomeric form of taufluvalinate. We also found two isomers of cyproconazole but they could not be separated.

Therefore, only 32 molecules could be quantified. Pesticide detection was performed with an Agilent 5977A mass spectrometer. The Electronic Impact (EI) mode was at +70 eV and a SIM/SCAN mode was used. The SCAN monitoring range was set from 30 to 510 m/z. Two ions per molecule were selected in the SIM mode (Table S. 2 and their isomers into three batches. Dilutions in acetone (Carlo Erba) of these commercial solutions were made to obtain standard concentrations from 0.01 to 2 ng.µL -1 (Table S.1). To quantify pesticides in rainwater, a calibration curve was built using Twister® extraction of 10 mL of ultrapure water (MilliQ, Merck) placed in a light-resistant flask and spiked with 10 µL of each of the three batches of standard solutions at increasing concentrations (0.01; 0.02; 0.04; 0.05; 0.10; 0.20; 0.40 ng.µL -1 ). These calibration samples were extracted at the same time as the rainwater samples, to obtain a calibration curve. For trueness recoveries (see section 2.4) below 80 %, the re To quantify pesticides in the air, clean Tenax® tubes were spiked with 1µL of standard solution (0.04; 0.1; 0.2; 0.4; 0.6 and 2 ng.µL -1 ) taken up with a glass syringe and deposited on the grid of the Tenax tube. The tubes were then dried to evaporate the acetone solvent with a helium stream at 0.5 ± 0.05 L.min -1 flow rate for 15 min. This was carried out three times to successively spike the same calibration Tenax tube with each standard batch. To quantify the mass of compound trapped, field tubes were analyzed on the same day as the calibration ones to ensure that the analytical conditions were the same. To evaluate the efficiency of pesticide extraction from rainwater, the trueness experimental recovery was determined as follows: some ultrapure water (MilliQ, Merck) (free from pesticides) and a rainwater sample taken in the ZAPVS (previously analyzed and found to be free of pesticides) were spiked with each of the studied molecules from a standard solution at a concentration of 0.4 ng.L -1 . The SBSE protocol established above was then carried out on two replicates of each water sample. Experimental trueness of the recovery of each pesticide studied was evaluated by comparing the quantities of molecules extracted from the rainwater (Prain) and the ultrapure water (Ppure). It was calculated as follows:

Eq. 2 and compared to the theoretical one calculated with Eq1.

To evaluate the efficiency of the chosen extraction conditions, following the experimental design results, we compared them to the modeled optimum for each compound by calculating a ratio (given as a percentage of the optimum, PO) based on the following equation: The limits of detection (LOD) and quantification (LOQ) were established for both air and water matrices. For air, we used a Tenax tube spiked with 1 µL of standard solution at 0.1 ng.µL -1 of all compounds. For water samples, we used 10 mL MilliQ water spiked with 10 µL of standard solution at 0.01 ng.µL -1 for all compounds. LOD and LOQ were calculated as signal to noise ratios of 3 and 10, respectively. Repeatability was also assessed for air and water samples (n=3) including the sample preparation steps. Results are presented in Table S.3.

Wet deposition calculation

Given the pesticide concentrations in rainwater ( in ng.L -1 ) and rain intensity for each sample ( in mm), wet deposition ( in µg.m -2 ) was calculated for each compound as cumulated over all samples in the nine months according to the following equation:

Eq. 4

2.6 Monitoring data from ATMO Nouvelle Aquitaine: sampling location and selection.

The AASQAs monitor and assess ambient air quality in every French administrative region (metropolitan and overseas). Even though pesticides are not subject to regulatory monitoring, some AASQAs carry out ad hoc and local campaigns to assess pesticide contamination of the atmosphere. Here we used the data measured weekly by Nouvelle Aquitaine (AASQA ATMONA) from February to November 2018. Active sampling was carried out using a PUF trapping media for the gas phase and a quartz filter (47mm) to trap the particles.

We selected four locations among the seven locations that the AASQA ATMONA monitors (Figure 2). These were Poitiers and Limoges (urban sites) as they are located within a similar agricultural profile as the ZAPVS studied (arable crops), as well as Medoc (a rural site vineyard)

and Cognaçais (a rural site, mixed environment of arable crops and vineyard) as they were located almost upwind of the ZAPVS and thus pesticide applications in this area could result in contamination of air masses near our site.

3 Results and discussion interactions had negative effects. To find optimal conditions, the of stat graphics was used.

Statistical models explained more than 90% of the variability in the quantity of the extracted molecules for each studied pesticide. The estimation of optimum extraction conditions was thus reliable. Thiamethoxam was the only molecule that could not be extracted by SBSE due to its very low octanol/water coefficient (LogKow = -0.13, Table 1): 31 molecules out of 32, which could be separated chromatographically (see 2.

Rainwater sample extraction efficiency

As explained previously, it was difficult to choose optimum conditions for each studied compound. The quantity of molecules extracted with the conditions fixed above represented more than 80% (PO) of the response that would have been obtained under the optimum conditions simulated by the models. Thus, the extraction conditions chosen above were a good compromise to efficiently extract all the molecules under study. Experimental recoveries (Eq 2) were higher than 80% for 22 out of the 31 molecules, and ranged from 40% to 59% for only eight molecules (Table S. 2: no available data for cypermethrin1). There was a significant increase between theoretical recoveries (Eq 1) and experimental recoveries (Eq 2), of two to ten fold for a majority of the molecules. However, five molecules and their isomers (cypermethrin, deltamethrin, lambda-cyhalothrin, bifenthrin, tau-fluvalinate) showed lower experimental recovery, two-fold lower than the theoretical yields, even though the experimental design clearly showed that the solvent conditions were the most favorable conditions for extracting them. These compounds all belong to the pyrethroid family. [START_REF] Serôdio | Development of a stir-bar-sorptive extraction-liquid desorption-large-volume injection capillary gas chromatographic-mass spectrometric method for pyrethroid pesticides in water samples[END_REF] observed a decreased signal with higher stirring speeds, one of the parameters that was not tested in the Box-Behnken design. The experimental extraction yield was taken into account by applying a correction factor for these five molecules in final calculations of their concentration in rainwater. SBSE was used here for 14 different chemical families covering a wide range of physico-chemical properties (log Kow ranging from 1.36 to 7.02). The choice of contrasted extraction conditions described above allowed the best compromise for the majority of the studied molecules.

All methodological tests and evaluations (LOD, LOQ, repeatability etc.) are summarized and presented in Table S. 3. These criteria were taken into account in the analysis and interpretation of the results.

Quantification of pesticides in rainwater samples

The cumulated rainfall amount over the experimental period was 621 mm, which is typical for the Deux-Sèvres department at this time of the year (ranging from 417.5 to maximum 1071 mm, in 2009 to 2019, as indicated in the INRAE CLIMATIK database) (Table 2). No rainwater was collected for nine of the weeks so a total of 22 samples was analyzed. The concentrations of 12 pesticides belonging to nine different chemical families were quantified in rainwater, ranging from 0.5 ng.L -1 to 173.9 ng.L -1 . All pesticide types were detected: six fungicides, five herbicides and one insecticide (Table 2). The highest concentrations were found for herbicides, which also showed the largest variability in concentration. Bifenthrin was the only insecticide detected with concentrations similar to those of fungicides. Over the entire period, rainwater contamination was mostly observed in spring (end of March-beginning of June) for both fungicides and herbicides (Figure S.3). The herbicides pendimethalin and diflufenican were also found in early March, in October and November. Note that winter was less well represented than the other seasons as no samples were taken before the 7 th of March 2018 and none after the 28 th of November 2018.

Pendimethalin and S-metolachlor were the most frequently detected pesticides in this study and are also among the most reported in rainwater with concentrations ranging from 160 to 350 ng.L -1 and from 100 to 510 ng.L -1 , respectively, in Europe (Dubus et al., 2000;Asman et al., 2005). Potter and Coffin (2017) observed even higher concentrations between 1.5 and 1700 ng.L -1 for pendimethalin and 15 to 3800 ng.L -1 for S-metolachlor. In another recent study, Kreuger et al. ( 2017) detected boscalid, epoxiconazole, diflufenican, pendimethalin, and Smetolachlor in 30% of the rainfall events between 2012 and 2015 in Sweden with concentrations ranging from 1 (most frequently detected pesticides) to 100 ng.L -1 (more occasionally detected pesticides), which is closer to the range of concentrations found in this study.

Quantification of pesticides in air samples

The developed analytical method (section 2.3) was carried out on the air samples collected in May 2018.

Six pesticides were quantified in air samples: three fungicides, two herbicides and one insecticide (Table 3). Concentrations ranged from 0.03 ng.m -3 to 1 ng.m -3 for the majority of the compounds except for pendimethalin for which concentrations were higher, between 1.5 ng.m -3 to 22 ng.m -3 . As for rainwater, we observed that herbicide concentrations were higher than those of fungicides. Tau-fluvalinate, the only insecticide detected, showed the same range of concentrations as fungicides. No significant variations in concentrations were observed during this short sampling period except for S-metolachlor in the organic farming window on the 30th and 31st of May, where there was a higher concentration for one sampling period over the four sampling periods (by a factor of three). This could be due to a S-metolachlor application close to the sampling location. Observed pesticide concentrations were relatively similar in the two studied windows, except for pendimethalin for which a higher concentration was observed in the CF (conventional farming) window (four fold higher). Six and four pesticides were detected in the OF (organic farming) and CF windows, respectively. The OF window was close to treated areas (Figure 1). The conventional farming fields may thus have had an impact on pesticides in the air above this area. However, the sampling was short (24 hours) and the number of samples collected was low, so these results should be considered with caution. It would be necessary to interpret them with regards to the actual treatments carried out within the studied 1 km 2 landscapes and at their edges (information not available at this scale).

Canadian or American studies published before 2013, for a large number of pesticides, including epoxyconazole, fenpropidine, pendimethalin and S-metolachlor, The concentrations quantified in our study are within the ranges reported by Coscolla et al. (2016). To go further on, we compared our results to the measurements performed by AASQA ATMONA at selected representative sites in 2018. Among the 27 active substances (not taking into account isomers), ten were not investigated in the AASQA ATMONA campaign: cloquintocet-mexyl, cycloxydim, dimethachlor, metconazole, metrafenone, napropamide, pirimicarb, pyraclostrobin, tau-fluvalnalate, and thiamethoxam. Among the 17 others, seven pesticides were quantified by AASQA ATMONA: aclonifen, boscalid, diflufenican, fenpropidine, metazachlor, S-metolachlor and pendimethalin. Two pesticides detected here (Table 3) were also quantified by AASQA ATMONA in May 2018, namely S-metolachlor and pendimethalin.

While S-metolachlor concentrations were consistent with the AASQA ATMONA data (0.34 to 1.32 ng.m -3 ), the pendimethalin concentrations measured in this study were higher (AASQA ATMONA data: 0.17 to 1.32 ng.m -3 ). Regarding fenpropidine, at the end of April, AASQA ATMONA reported a concentration of 0.22 ng.m -3 in the air, which is very close to that we measured. AASQA ATMONA did not detect epoxiconazole (LOQ < 0.15). The overall agreement between this study and AASQA ATMONA shows that the method developed in the present work is relevant for pesticide quantification in the atmosphere. These findings correspond to farmer pesticide usage recorded through interviews in the study site. For instance, pendimethalin and epoxiconazole were used in 23% and 15% of the fields for which data on farming practices were collected in 2018 (see Bretagnolle et al. 2018a,b for details on inquiries).

We found only one surprising discrepancy for S-metolachlor. From the interviews, few farmers used it (< 5% of fields for which data on agricultural practices are known) while it was frequently found in the atmosphere in this study. This pesticide is mostly used in maize for weeds. If maize is not a particularly frequent crop at ZAPVS scale (c.10% of the area, representing 8% of agricultural practices datas collected from farmer interviews), it was particularly abundant in this conventional farming landscape studied and sampled here (>50%), which could explain our observations.

Seasonality of rainwater and airborne concentrations

Some pesticides (aclonifen, diflufenican, pendimethalin, S-metolachlor) were detected in rainwater and reported by AASQA ATMONA in the atmosphere over the same periods. Both datasets can be used to analyze the seasonal trend of these pesticides in the region as illustrated in Figure 4 for pendimethalin and S-metolachlor, the two most frequent compounds quantified in air by AASQA ATMONA and in rainwater in this study.

Overall, there was a good correlation between the seasonal variations in pesticide concentrations measured in the rainwater and atmosphere. This is expected because of the equilibrium between gaseous compounds and their solubilized counterparts in water. The ratio between rainwater and atmospheric concentrations (Expratio) was calculated and evaluated theoretically as proposed by [START_REF] Ligocki | Trace organic compounds in rain III. Particle scavenging of neutral organic compounds[END_REF] in order to parameterize the scavenging of gaseous organic compounds by rain. Two estimates of this ratio were calculated either based (T1ratio), or on the ratio between vapour pressure and water solubility (T2ratio). This approach assumes that all the pesticide is in the gas phase. Theoretical and experimental ratios were in the same order of magnitude except for those for pendimethalin, which showed a much larger experimental ratio (Figure S.4). These differences may be due to the uncertainty in the physico-chemical properties of these compounds. Indeed, Van Pul et al.

(1999) evaluated that the uncertainty in wet deposition of gaseous pesticides is a factor of two to three, mainly because of uncertainties in dependency. The observed differences may also be due to the assumption that the compounds are entirely in the gas phase, while a fraction may be in the particulate phase, thus affecting the scavenging of the compounds by rain. For further comparisons of the experimental ratios to the theoretical ratios, the gas-particle partitioning of pesticides should be taken into account, which is still a challenge.

3.5 Wet deposition and consequences for soil contamination Kreuger et al. found similar orders of magnitude (2017) for propiconazole (somewhat lower than 1 µg.m -2 ) and epoxiconazole (close to 1 µg.m -2 ). Higher deposition was found here for pendimethalin and S-metolachlor comp this is consistent with the fact that the use of pendimethalin and S-metolachlor were forbidden in Sweden when the samples were taken while they are still allowed and used in France. Our estimates were however lower than the annual mean deposition of 92 and 80 µg.m -2 found by [START_REF] Potter | Assessing pesticide wet deposition risk within a small agricultural watershed in the Southeastern Coastal Plain (USA)[END_REF] for pendimethalin and S-metolachlor, respectively. In terms of application dose (based on the recommended application dose), for pendimethalin this deposition rate corresponded to 0.03% of the application dose being deposited over the nine months. This result is in agreement with these assumptions, the yearly contribution of wet deposition to soil concentrations was evaluated to range between 1% and 2.4% of the pesticides found in untreated soil by [START_REF] Pelosi | Residues of currently used pesticides in soils and earthworms: A silent threat?[END_REF]. The rainwater was sampled in an area of the ZAPVS that is not located exactly in the same fields as those sampled by [START_REF] Pelosi | Residues of currently used pesticides in soils and earthworms: A silent threat?[END_REF] and not in the direct vicinity of treated fields, thus we may have underestimated the contribution of rainwater to the level of soil contamination observed by [START_REF] Pelosi | Residues of currently used pesticides in soils and earthworms: A silent threat?[END_REF]. Spatial and temporal synchronization of sampling is necessary to go further with this analysis, which should also be extended to a wider range of compounds. In addition, untreated areas may also be contaminated by gaseous deposition of pesticides in the air or by droplets drifting from treatments into adjacent treated areas. Modeling studies taking into account local plant protection practices and landscape features are required to assess the contribution of each process.

Conclusion

This study developed a reliable method for quantifying mixtures of several pesticides in rainwater and air. The proposed methodology based on using a single analytical method for both compartments represents a clear advantage: sources of uncertainties associated with the methods involving several laboratories and techniques can be reduced. The use of Box-Behnken factorial experiments made it possible to identify the best compromise in conditions to extract the 27 active substances in rainwater. The descriptive quality of the statistical models obtained from this experimental design made it possible to model the influence of analytical choices on extraction yields without conducting new laboratory experiments that would be time-, material-and money-consuming. Using this statistical tool, we could study the influence of several factors in a minimum number of experiments in the laboratory. With regards to the gas phase, a sampling time of 24h may be possible to quantify atmospheric concentrations, in view of the limits of quantification obtained, thus allowing the relatively fine monitoring of the temporal evolution of concentrations, which may be rapid. TD-GC-MS appears to be a good analytical solution to detect and quantify pesticides in both air and rainwater. The limits of quantification obtained were low and could be further improved by avoiding molecule splitting towards the analysis system. Nevertheless, the proposed approach also has some limits. For example : TD is not suitable for thermolabile compounds and the sample is destroyed so additional analyses of the same sample are not possible , GC is intended for volatile or semi-volatile molecules.

The feasibility of air monitoring using an approach based on a mobile active sampling system that does not require a nearby power supply was confirmed by comparing the concentration levels with AASQA ATMONA regional measurements in the region. To assess the involvement of dry deposition in the contamination of untreated areas, a broader-scale deployment of the measurements within agricultural landscapes would be required. Developed air sampling systems that are self-contained, particularly in terms of energy supply, would meet this need.

For wet deposition, a seasonal effect was clearly observed showing a correlation between molecules detected in rainwater and air. The differences between theoretical and measured airto-rainwater concentration ratios raise the question of uncertainties in rain scavenging processes and points towards the need for further monitoring in field conditions. Using the same analytical approach for both compartments is a real advantage for addressing this issue.

25.4

Table 2 : Concentrations of pesticides detected in rainwater in ng.L -1 , rainfall amount in mm, concentration minimum and maximum quantified during sampling period in ng.L -1 and quantification frequency in % (ratio number of times detected /number of rain events) . ND : not detected;

<LOQ : detected but below LOQ; "-" : no rainfall or sample volume collected was too small for analysis.

  Some factors were fixed according to Prieto et al. (2010), and based on our instrumentation: the stirring speed (800 rpm), the extraction temperature (25°C), the sample volume (10mL) and the acceptor phase volume (24 µL PDMS). A Box-Behnken experimental design (Box and Behnken, 1960) was conducted to study the influence of the following factors : the addition of an inert salt (NaCl, Factor A), the addition of an organic modifier (Acetone, Factor B) and the extraction time (Factor C). Box-Behnken designs are based on three-level incomplete factorial designs. A series of experiments was set up to study the relationships between factors and response. The experimental domain and the Box-Behnken design used are detailed in Table S.1. Statistical analysis was performed to model the response Y (peak areas corresponding to the quantity of extracted molecules) as a function of the chosen factors A, B, C (extraction conditions). The obtained Pareto diagram was used to identify which factors significantly affected the quantity of extracted molecules. The statistical model allowed the factor levels that maximized the response (simulated optimum conditions)

  ). The temperature of the ion source and Quadruple analyzer were set at 230 and 150°C, respectively. TD-GC-MS data were processed with MassHunter software (version B.07.04.5560, Agilent Technologies Inc.) ensure an accurate identification and quantification of the studied molecules, we worked on the EIC (Extracted Ion Chromatogram) SIM chromatograms. The extracted ion chromatogram is generated by separating the ions of interest from the SIM data file (Figure S.1). All the chromatographic peaks were checked one by one by analyzing (a) the retention time of the molecules (b) the presence for each molecule of its specific ions (c) the intensity ratios of the 2.3.3 Quantification step Standard solutions at 100 ng.µL -1 in acetone of each studied molecule were purchased from Restek. To avoid inter-molecule reactivity, the supplier suggested separating the 27 pesticides

2. 4

 4 Method characterization SBSE theoretical recovery of a compound from a water sample was calculated with the following equation (Baltussen et al.,1999): Eq. 1 where is the mass of analyte in the SBSE, is the amount of analyte present in the water sample, is the octanol-water partitioning coefficient, and are the water sample volume and the volume of polydimethylsiloxane (PDMS) phase on the Twister ® , respectively.

  the peak area representing the quantity of molecules extracted in the first extraction, the peak area representing the quantity of molecules extracted in the second extraction, the peak area simulated representing the quantity of molecules extracted under simulated optimum conditions. To validate the effectiveness of standard tube preparation for air extractions, sequential doping (three different spikes on the Tenax tube) as described in section 2.3.3 was compared to the routinely used method for doping (only one spike on Tenax tube) (Figure S.2). The thermal desorption efficiency step was also tested by analyzing the same sample twice and calculating the ratio of the amount of compound measured after the second desorption to the amount measured with the first desorption.

3. 1

 1 Method development and characterization3.1.1 Experimental design results: choice of conditions for rainwater sample extractionTableS.4 summarizes the results obtained for all compounds and Figure3shows the experimental design results for one compound, bifenthrin, as an illustration. The Pareto chart (Figure3a) represents the standardized effect (means of the effect divided by its standard error) of the factors, and ranks the factors in order of importance. For bifenthrin, three factors were statistically significant. The extraction time (C) and the addition of acetone (B) had a positive effect, i.e. they improved its extraction, while the addition of NaCl (A) had a significantly negative effect. One of the advantages of the experimental design approach was also to study factor interactions. For bifenthrin, NaCl/Extraction time (AC) and NaCl/Solvent (AB)

  3.2 subsection) were thus considered. According to the optimum results of the Box-Behnken experiment, the first factor to enhance extraction of a majority of molecules was the addition of salt (30%) (13 molecules out of 31) followed by the addition of acetone (10 %) (nine out of 31 molecules). These two extraction conditions seemed to be therefore unavoidable. The results were less contrasted for the remaining molecules. For clomazone, pirimicarb, and propiconazole 1, adding salt had a significant positive effect whereas acetone had a negative effect. Adding acetone did not have an effect on cloquintocet-mexyl and diflufenican extraction. Salt addition had no effect on the extraction of pyraclostrobin and metrafenone. Each molecule would therefore be extracted efficiently in both cases. Overall, the addition of salt and acetone appeared to be a good compromise in the choice of extraction conditions.The results of Pareto diagrams also indicated that NaCl/acetone interactions had a significant negative effect on twelve molecules. Therefore, the two extraction conditions were not carried out on the same sample at the same time. However, since it is possible to analyze two Twister® bars together in the same TDU tube, the two extraction conditions were carried out on two separate aliquots of the rainwater samples. Concerning the extraction time, the simulations predicted that more than 200 minutes would lead to optimum extraction. However to carry out all the extractions in a reasonable time with good repeatability, a compromise was sought to reduce the extraction time. Based on the simulations, 135 min was found as a good compromise that only reduced the quantity of extracted molecules by less than 20% compared to the optimum.The following extraction conditions were hence chosen: 30% (m/v) NaCl was added to 10 mL of rainwater sample and a first extraction was carried out for 135 min (1st Twister®). 10% (v/v) acetone was added to a second aliquot of 10 mL rainwater sample and a second extraction was carried out for 135 min with a second stirring Twister®. Both Twister® bars were extracted and analyzed simultaneously by TD-GC-MS.

  Kreuger et al. (2017) who found that deposited amounts corresponded to 0.1-0.0001% of the applied dose in the field.[START_REF] Potter | Assessing pesticide wet deposition risk within a small agricultural watershed in the Southeastern Coastal Plain (USA)[END_REF] found 1.2% and 1% for pendimethalin and S-metolachlor, over the agricultural watershed over 3 years.[START_REF] Pelosi | Residues of currently used pesticides in soils and earthworms: A silent threat?[END_REF] detected diflufenican, epoxiconazole and boscalid in the top 5 cm of soils (sampled in 2016) in untreated areas such as grasslands (temporary and permanent) and hedgerows. The lowest concentrations (between 0.3 and 0.5 ng.g -1 dry soil) were found in permanent grasslands (untreated). Combining our results and those from[START_REF] Pelosi | Residues of currently used pesticides in soils and earthworms: A silent threat?[END_REF], bearing in mind that the sampling years were different in the two studies, we evaluated the contribution of wet deposition to soil contamination based on the following assumptions. First, since the sales of boscalid, diflufenican and epoxiconazole were relatively constant in the Deux-Sèvres to 2018 (BNVD 2021), we assumed that the results of both experimental campaigns are representative of the overall contamination of rainwater and soils. Second, we considered that wet deposited pesticide molecules remained in the top 5cm of soil (depth sampled by Pelosi et al. (2021)), thus ignoring further infiltration and degradation given the relatively high DT50lab (time required for the concentration to decline to half of the initial value), of each compound in the soil (DT50lab,: 94.5, 353.5 and 484.4 days for diflufenican, epoxiconazole and boscalid respectively, from PPDB). Based on

Figure 1 :

 1 Figure 1: Spatially nested maps showing the location of A) the ZAPVS study site, B) the

Figure 2 :Figure 3 :

 23 Figure 2: AASQA ATMONA map, location of the ZAPVS and the different monitoring sites

Figure 5 :

 5 Figure 5: Cumulated wet deposition (grey bars) over all samples (µg.m -2 ) from March to

Table 1 .

 1 Céline Décuq a , Marjolaine Bourdat-Deschamps a , Pierre Benoit a , Colette Bertrand a , Rachid Benabdallah a , Baptiste Esnault a , Brigitte Durand a , Benjamin Loubet a , Clémentine Fritsch b , List of the studied pesticides and their physico-chemical properties, type and chemical substance family. (PPDB Pesticide Properties DataBase University of Hertfordshire 2020)

	5 Acknowledgments :
	metaprogram INRAE SMaCH Call 2017. The study also benefited from the samples collected
	d Food and the

https://doi.org/10.1016/j.apr.2020.03.014 assessing pesticide atmospheric contamination in untreated areas

Ministry for an Ecological and Solidary Transition, with the financial support of the French , using the fees for diffuse pollution from the Ecophyto Plan through the national agency ONEMA. We thank the ZAPVS for helping with access to field plots, especially Jean-Luc Gautier, CEBC, in collecting rainwater samples and A. Hulin from Atmo Nouvelle Aquitaine for providing information on AASQA ATMONA Data. We thank also Leigh Gebbie for her english reviewing and Pascal weekly by AASQA ATMONA at the same sampling period, on the site selected in section 2.5.

ND: not detected, NR: not researched, LOD: limit of detection.

Declaration of interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:


Guilhem Madiot

Sylvain Barbay

Rémy Braive
email: remy.braive@c2n.upsaclay.fr

Vibrational resonance amplification in a thermo-optic optomechanical nanocavity

Vibrational resonance amplifies a weak low-frequency signal by use of an additional non-resonant high-frequency modulation. The realization of weak signal

enhancement in integrated nonlinear optical nanocavities is of great interest for nanophotonic applications where optical signals may be of low power. Here, we report experimental observation of vibrational resonance in a thermo-optically bistable photonic crystal optomechanical resonator with an amplification up to +16 dB. The characterization of the bistability can interestingly be done using a mechanical resonance of the membrane, which is submitted to a strong thermoelastic coupling with the cavity.

Phase transitions and double-well potentials have been extensively exploited to amplify or detect weak signals. [START_REF] Rajasekar | Nonlinear resonances[END_REF] Such general physical concepts, observed in various fields of science are at the heart of the phenomenon of vibrational resonance 2 (VR). Introduced as an analogy with the well known stochastic resonance, [START_REF] Dykman | Stochastic resonance in perspective[END_REF] VR uses a high-frequency (HF) periodic signal to amplify a low-frequency (LF) input signal. It has been theoretically studied in different types of nonlinear systems, e.g. in neural network, [START_REF] Deng | Vibrational resonance in neuron populations[END_REF] in excitable systems [START_REF] Zaikin | Vibrational resonance in a noise-induced structure[END_REF] or in biological networks. [START_REF] Rajasekar | Vibrational resonance in biological nonlinear maps[END_REF] Several experimental demonstrations have also been conducted in electronic circuits, [START_REF] Ullner | Vibrational resonance and vibrational propagation in excitable systems[END_REF] bistable VCSELs [START_REF] Chizhevsky | Experimental Evidence of "Vibrational Resonance" in an Optical System[END_REF][START_REF] Chizhevsky | Experimental evidence of vibrational resonance in a multistable system[END_REF] or electromechanical Duffing resonator, [START_REF] Chowdhury | Weak signal enhancement by nonlinear resonance control in a forced nano-electromechanical resonator[END_REF] for example. Despite the growing interest for photonic nanodevices and their use for signal processing and sensing, VR has not yet been demonstrated using nonlinear optical nanocavities.

Among optical nonlinearities, the thermooptic effect overwhelmingly manifests itself in photonic micro and nano-devices. Resulting from the temperature dependence of a material refractive index, it is commonly utilized to introduce tunability and to realize elementary computational all-optical components, such as photonic switches, [START_REF] Watts | Adiabatic thermo-optic Mach-Zehnder switch[END_REF] phase-shifter [START_REF] Harris | Efficient, compact and low loss thermo-optic phase shifter in silicon[END_REF] or integrated spectrometers. [START_REF] Souza | Fourier transform spectrometer on silicon with thermo-optic nonlinearity and dispersion correction[END_REF] Due to the ultimate electromagnetic field confinement achievable in optical nanocavities, strong thermo-optic effects can lead to multistability, [START_REF] Uesugi | Investigation of optical nonlinearities in an ultra-high-Q Si nanocavity in a two-dimensional photonic crystal slab[END_REF][START_REF] Gao | Cavity-enhanced thermo-optic bistability and hysteresis in a grapheneon-Si3N4 ring resonator[END_REF][START_REF] Gray | Thermo-Optic Multistability and Relaxation in Silicon Microring Resonators with Lateral Diodes[END_REF] which has been widely used to study nonlinear dynamical behaviors such as, e.g., self-stability, [START_REF] Carmon | Dynamical thermal behavior and thermal self-stability of microcavities[END_REF] selfpulsing, 18 excitability 19 or soliton modelocking. [START_REF] Joshi | Thermally controlled comb generation and soliton modelocking in microresonators[END_REF] The usefulness of optical nanocavities has also largely been demonstrated to sense or manipulate mechanical vibrations. [START_REF] Aspelmeyer | Cavity optomechanics[END_REF] Thus complex dynamics can emerge due to a combination of thermo-optic and optomechanical nonlinearities enabling e.g. utilization in sensing applications, [START_REF] Deng | Thermooptomechanical oscillator for sensing applications[END_REF] chaos generation, [START_REF] Navarro-Urrios | Nonlinear dynamics and chaos in an optomechanical beam[END_REF] or electrooptomechanical self-oscillation [START_REF] Allain | Electro-Optomechanical Modulation Instability in a Semiconductor Resonator[END_REF] .

While stochastic resonance has been observed in optomechanical systems, [START_REF] Mueller | Optomechanical stochastic resonance in a macroscopic torsion oscillator[END_REF][START_REF] Monifi | Optomechanically induced stochastic resonance and chaos transfer between optical fields[END_REF] VR remains unexplored both theoretically and experimentally. Within this framework, VR in opti- cal nanocavities supporting simultaneously mechanical modes and nonlinearities is a key development for potential nanophotonic and optomechanical applications.

In this letter, we demonstrate VR amplification of a weak LF signal in a suspended thermo-optic photonic crystal (PhC) nanoresonator in the bistable regime. Thanks to a strong thermo-elastic coupling between one mechanical resonance of the PhC membrane and the optical cavity load, the state of the thermooptic cavity can be readout through the mechanical mode spectral position. It further enables to characterize the threshold modulation amplitude below which the modulation signal does not trigger switching between the states of the bistability. Then by adding a non-resonant HF signal, we observe up to +16 dB amplification of the weak LF signal.

The system is a 12 × 21 µm 2 rectangular InP membrane with thickness 265 nm suspended over a 250 nm airgap. The membrane is pierced with a 2D PhC in which four L3 defectcavities are diagonally arranged (fig. 1a) such that the electromagnetic field confined in each defect can evanescently couple to the neighboring ones. [START_REF] Chalcraft | Mode structure of coupled L3 photonic crystal cavities[END_REF] A SOI waveguide is placed underneath the membrane for integrated light injection. [START_REF] Tsvirkun | Integrated III-V Photonic Crystal-Si waveguide platform with tailored optomechanical coupling[END_REF] The electromagnetic modes confined in the PhC molecule are characterized by injecting a broadband superluminescent diode through the waveguide. The IR spectral analysis of the waveguide transmission, shown in fig. 1b, evidences four resonance dips. Here, the use of a PhC molecule enables higher optomechanical coupling to the MHz drum modes of the membrane -i.e. whose associated displacement field is extended on the whole membrane -and on which we focus here.

The PhC eigenmodes can alternatively be injected with a preliminary amplified tunable laser (see fig. 1c) enabling coherent and resonant excitation of the photonic modes. Doing so, the brownian motion of the suspended membrane is optomechanically transduced into the optical field. The analysis of the output field noise spectrum is done via an electrical spectrum analyzer (ESA). The membrane carries several mechanical modes in the MHz domain. Focusing on the 6 MHz fundamental mechanical mode, we plot the noise spectrum of the output optical field as a function of the input laser wavelength in fig. 1d. The laser is scanned over the 2 nd and 3 rd photonic modes and the estimated power sent to the photonic molecule is P in = 2.4 mW. The use of sufficiently high circulating power reveals a thermo-optic nonlinearity which leads both photonic modes to shift and to become bistable. Note that the 2 nd mode experiences a stronger red-shift than the 3 rd mode. The thermal origin of this nonlinear effect is confirmed in the following by the temporal response of the optical cavity. Interestingly, the intra-cavity field intensity is partially absorbed which also provokes the photothermal shift of the mechanical frequencies through the thermo-elastic effect. Indeed, as the laser wavelength is upwardly swept from off-resonance to on-resonance the mechanical frequency experiences a linear photothermal shift. The abrupt intra-cavity intensity drop occurring when passing the edge of the first bistability (2 nd photonic resonance) translates into an abrupt red-shift of the mechanical frequency (black arrow). The same phenomenon happens when scanning the 3 rd photonic bistable resonance (fig. 1d). The downward scan can only be performed manually and gives access to the hysteresis turning point. The backward jumps of the mechanical frequency are measured and shown with the gray arrows. The corresponding trajectory superposes with the upward trajectory except in the bistable regions. Therefore it is possible to evaluate the intra-cavity field intensity by measuring the frequency shift of a mechanical mode. In the following we set the laser wavelength at the center of the second photonic mode bistability (λ = 1566.75 nm) to symmetrize the double-well potential. [START_REF] Chowdhury | Phase Stochastic Resonance in a Forced Nanoelectromechanical Membrane[END_REF] Thus we name hot (blue) and cold (green) mechanical states the corresponding mechanical frequencies Ω hot m /2π and Ω cold m /2π (see inset in fig. 1d). A lorentzian fit of the mechanical lineshape (white lines) returns a mechanical linewidth Γ m /2π = 4.3 kHz, hence a mechanical quality factor of about 1400.

When the optical resonator state flips, it requires a certain time to reach its new stable state. Under strong modulation of the input field intensity, the optical state is asked to flip twice per modulation period, which is possible only if the switching time τ s is short enough. Therefore, the maximum modulation frequency allowing the system to accurately flip is given by 2τ -1 s . We estimate this cut-off modulation frequency by measuring the typical switching time of the system. To do so, we modulate the input laser field slowly enough such that the transitory regime -in which the system leaves a stable state to reach the other one -can be measured. The laser wavelength is modulated in an intensity electro-optic modulator (EOM, with half-wave voltage V π = 7.5V ) on which we apply a square electrical signal carrying amplitude V mod and frequency Ω mod . A femtowatt sensitivity photodetector converts the transmitted laser field into a voltage V out we input in an oscilloscope (see fig. 1c).

Modulating the input field with V mod = 1 V and Ω mod /2π = 10 kHz, the system response is recorded over several hundred of modulation cycles. The data are averaged cycle by cycle and plotted in fig. 2a. During one cycle, the resonator, initially set in the hot state, transits towards the cold state and then switches back to the hot state at half a modulation cycle. Note that since the cavity load is accessed through the waveguide transmission, the cavity hot state (resp. cold state) corresponds to the low (high) transmission level. This transient regime manifests as an exponential decay from one state to the other. Fitting the averaged data provides the switching time τ s ≈ 4 µs. This value is similar for both transitions hot→cold or cold→hot. The corresponding cutoff frequency, as discussed above, is of the order of 125 kHz, which is consistent with previous estimation made in a similar InP photonic crystal nanocavity. [START_REF] Brunstein | Thermo-optical dynamics in an optically pumped Photonic Crystal nano-cavity[END_REF] In order to clearly resolve the mechanical frequency, we use in the following a LF weak-signal with frequency of Ω mod /2π = 10 Hz. In principle, the use of faster signal should not prevent from the realization of amplification as demonstrated below. However, the modulation tones present in the driving field couple with the optomechanical cavity such that they are transduced into the mechanical frequency domain. [START_REF] Pelka | Floquet control of optomechanical bistability in multimode systems[END_REF] As a result, the mechanical response at resonance would be weakened as sidebands would be created. Keeping Ω mod Γ m guarantees here the demonstration of VR.

Once the laser wavelength and modulation frequency are set, the switching between optical stable states can be achieved only if the modulation amplitude passes a certain threshold. The demonstration of VR amplification relies on the use of a weak signal, i.e. with amplitude lying below this threshold, whose knowledge is consequently required. To evaluate the latter, we increase the modulation amplitude and check the optical response. When the amplitude is sufficiently high, the optical state flips periodically at the modulation frequency.

In our case, the optical state can be accessed both through the transmission signal or through the spectral position of the mechanical resonance. In parallel of the direct output measurement (V out ), the output signal is analyzed in the ESA via a fast photodetector. From the spectrum, we determine the two positions of the mechanical peak Ω cold m /2π = 6.0572 MHz and Ω hot m /2π = 6.0164 MHz for the cold and hot optical states, respectively. The RF signal is therefore amplified and demodulated using a lock-in amplifier. Two demodulation channels are used at frequencies Ω cold m and Ω hot m and with 1 kHz wide passband filters. The two corresponding demodulated signal amplitudes V cold m and V hot m are recorded in the oscilloscope. Due to the demodulation filter whose bandwidth is smaller than the mechanical linewidth, noise in the demodulated signal is induced by frequency fluctuations of the mechanical resonance.

The resulting time traces are shown for three values of the modulation amplitudes in fig. 2b. For each situation, we show the reference mod-ulation signal (black), the transmission signal (red) and the demodulation amplitudes (blue and green). Around a threshold value V t = 980 mV, the optical state starts to flip but tends to remain in the cold state. The dynamics of the mechanical frequency is in perfect correlation with these jumps. For V mod > V t , the optical state jumps are perfectly synchronous with the modulation reference signal. Note that the signal-to-noise ratio in the mechanical spectrum is higher in the hot state, due to higher optomechanical coupling which is proportional to the energy stored in the optical cavity. This difference is also visible in the amplitude of the mechanical response (inset of fig. 1d).

For a given time trace, one can calculate the residence probability of the optical state. It consists in evaluating the amount of time spent by the system in the cold (or hot) optical state. To do so, a threshold line is arbitrarily chosen in between the two corresponding amplitude levels. The probability for the system to set in the initial state, is 100% for low modulation amplitudes (see fig. 2c). It quickly converges to 50% around the threshold amplitude V t = 980 mV. Such data treatment performed on V out , V hot m or V cold m all provide the same residence probability shown in fig. 2c.

With the full knowledge on the cut-off frequency and the threshold modulation amplitude, we investigate the amplification of a weak signal via an additional HF intensity modulation with frequency Ω hf and amplitude V hf . Thus the field now carries both tones Ω mod and Ω hf :

V ref = V mod sign(sin(Ω mod t) + V hf cos(Ω hf t)
When V mod < V t and V hf = 0, we know that the square signal cannot flip the optical state. However, the increase of V hf comes with a significant distortion of the bistability until the switching process can be re-activated. Ensuring that Ω mod Ω hf , 10 we arbitrarily set the frequency of the HF signal to Ω hf /2π = 80 kHz and V mod = 250 mV, and record the output optical response for increasing value of V hf . The results are shown through three examples. In each case, we plot the time trace (fig. 3a) and

1 0 -3 1 0 -2 1 0 -1 1 0 0 1 0 -3 1 0 -2 1 0 -1 1 0 0 0 . 9 1 . 0 1 . 1 1 0 -3 1 0 -2 1 0 -1
V h f = 4 0 0 m V b ) a ) 0 . 0 0 . 1 0 . 2 0 . 3 V o u t ( V ) 0 2 4 6 8 1 0 0 . 0 0 . 1 0 . 2 0 . 3 V h f = 7 5 0 m V ( Ω m o d / 2 π) ×t r e f V h f = 6 4 0 m V F F T a m p l i t u d e ω/ Ω m o d c ) S N R ( d B ) V h f ( m V ) V m o d 1 5 0 m V 2 5 0 m V + 1 6 d B
Fig. 3: a) Time-domain and b) frequencydomain transmitted signal V out for three values of the high frequency signal amplitude V hf = 400, 640 and 750 mV. c) SNR as the function of the HF signal below and above the amplification threshold, resp. with V mod = 150 mV (gray) and V mod = 150 mV (black). We use Ω hf /2π = 80 kHz.

the corresponding FFT spectrum (fig. 3b). At frequency Ω mod , a peak testifies the presence of the LF square modulation. Its amplitude increases when the switching process becomes more efficient. While at V hf = 400 mV the system remains in the hot state, switching is obtained at 640 mV, with 100% fidelity to the LF signal. In both cases, the noise remains constant in the FFT spectrum. At V hf = 750 mV, the fidelity drops which results from a clear increase of the noise floor. For a given LF signal amplitude V mod , the signal-to-noise ratio (SNR) is evaluated as a function of the HF signal amplitude V hf . [START_REF] Chizhevsky | Experimental Evidence of "Vibrational Resonance" in an Optical System[END_REF] The SNR is taken as the amplitude of the LF modulation peak in the FFT spectrum, to the noise level around this peak. This measurement is reproduced for several values of V mod . The resulting amplification curves are shown in fig. 3c. Amplification occurs when the SNR is higher than its value for V hf = 0. These reference levels are indicated with dashed lines. We observe a significant amplification (V mod = 250 mV, black curve) up to +16 dB. The range of V hf -in which amplification is observed -is delimited by transitory regions where the switching is imperfectly achieved, resulting in a high noise level, as discussed above. Weakening further the LF signal amplitude modifies amplification features as observed with the gray curve in fig. 3c, where we use V mod = 150 mV. An amplification of 5 dB is observed around V hf = 900 mV. It is likely that higher amplification could be obtained by increasing V hf .

The observation of VR reported here can be consistently modeled with a single waveguidecoupled optical mode including a photothermal nonlinearity. Numerical simulations, available in the Supporting information, provide a qualitative agreement with the experimental findings.

VR amplification could not be performed by exploiting the demodulated signals, V hot m and V cold m at the mechanical resonances because of the appearance of sidebands through the HF tone. Indeed, the presence of the latter (with Ω m > Ω hf Γ m ) can produce strong sidebands reducing the effective mechanical response at resonance, and whose intensity depend on the modulation depth and frequency. [START_REF] Pelka | Floquet control of optomechanical bistability in multimode systems[END_REF] Therefore, optimal adjustment of these parameters is expected to render possible the characterization of VR in the mechanical domain. Other solutions might also be investigated. Self-sustained oscillation regimes enabled in high optical quality factor and low phase-noise optomechanical cavities 32 would constitute an ideal support to explore manipulation of mechanical oscillations through the presently described phenomenon of vibrational resonance.

In conclusion, we have demonstrated amplification of a weak LF signal using vibrational resonance in a waveguide-coupled thermo-optic nano-resonator. The characterization of the bistable system is made both through the waveguide transmission optical field and through the thermo-mechanical effect induced by photo-thermal absorption in the optomechanical cavity. This strong thermo-elastic coupling could be exploited in sensing applications [START_REF] Deng | Thermooptomechanical oscillator for sensing applications[END_REF] or in the realization of FM modulation of nanomechanical oscillators, as demonstrated here. Our system made of coupled PhC cavity constitutes an ideal test-bed for such investigations. Additionally, the amplification of weak signals could be enhanced and enriched by the presence of multistability, which can occur when two or more bistable resonances overlap. [START_REF] Gray | Thermo-Optic Multistability and Relaxation in Silicon Microring Resonators with Lateral Diodes[END_REF] approach:

ȧ(t) = j(ω L -ω c ) - κ t 2 a(t) + √ κ e × s i (t) (1)
where κ t = κ i + κ e is the optical mode total resonance linewidth. We introduce the thermo-optic nonlinearity through the cavity temperature dynamics, with ∆T the temperature elevation in the optical cavity:

∆T = K t κ th |a| 2 -∆T (2)
with K t the thermalization rate of the cavity and κ th the thermal resistance (in K/J). The optical mode resonance wavelength is thermally shifted by an amount:

∆λ = λ 0 n 0 dn dT ∆T (3)
with λ 0 = 2πc/ω 0 the optical mode wavelength -ω 0 and n 0 being respectively the cavity natural frequency and refractive index, both taken at room temperature -and dn dT the thermo-optic coefficient of the material. In the numerical simulation, we integrate eqs. (1) to (3) using the real and imaginary parts of a, and considering a temperature dependent resonance frequency for the optical mode:

ω c = ω 0 (1 - 1 n 0 dn dT ∆T ) (4)
The ODEs are integrated after timenormalization t -→ K t t, as K t constitutes the cut-off frequency of the thermo-optic nonlinearity. The input field is modelled with a square modulation at Ω mod plus a high frequency sinusoidal signal at Ω hf . Both components have respective amplitudes given by the modulation depths, respectively α mod and α hf : s i (t) = P in 1 + e jπ -1 2 +α mod f (Ω mod t)+α hf cos(Ω hf t) (5) where P in is the laser input power and f (t) is the square signal function,

f (t) = 4 π N p=0 sin (2p + 1)t 2p + 1
In the numerics we use N=15. All the constants are given in table 1.

Numerical simulations

In fig. 4a the mode intensity is numerically evaluated as a function of the laser wavelength under forward and backward scans (resp. black and red curves). Tuning the input power P in , the spectral span (∼ 1 nm) of the bistability (red region) is adjusted to qualitatively match the experiment. Note that playing with the power of with the thermo-optic coefficient is equivalent in this model. We find the best agreement for P in = 10 mW.

Using the laser wavelength λ L = 1566.3 nm), we numerically resolve the ODE for increasing modulation depth α mod . The modulation frequency is set to Ω mod /2π = 500 Hz rather than 10 Hz in the experiments, in order to reduce the integration time. The time trace |a| 2 (t) is used to calculate the residence probability of the optical mode (see fig. 4b).

Enabling the high-frequency modulation with Ω hf /2π = 10 kHz, the amplification curve is obtained by evaluating the signal-to-noise ratio (SNR) of the weak low-frequency modulation in the mode response, as a function of the HF signal amplitude, α hf . The results are shown for two different values of α mod , in fig. 5. Significative amplification is obtained for α mod = 0.005 (black curve), but not for α mod = 10 -5 (gray Note that the model here accounts for a single optical mode. Therefore the significative overlap between bistable optical modes that we observe experimentally is not captured by the present simple model.

Fig. 1 :

 1 Fig. 1: a) SEM image of the suspended photonic molecule -made of L3-defect cavities (arrows)and underneath SOI waveguide (red stripe). b) Waveguide transmission spectrum with four dips corresponding to the PhC molecule photonic eigenmodes. c) Experimental setup. d) Noise spectrum of the output optical field as a function of the injected laser wavelength. Scan performed over the 2 nd and 3 rd photonic modes. The fundamental mechanical resonance experiences a photothermal shift proportional to the intracavity energy. The black (resp. gray) arrows indicate the abrupt jumps of the mechanical frequency observed at the bistabilities' edges under upward (backward) scan of the laser wavelength. inset: noise spectra taken when the optical resonator is the hot state (blue) and in the cold state (blue). The mechanical quality factor Q m ∼ 1400 is estimated from the Lorentzian fits (white).

(

  Fig.2: a) Output signal (V out ) averaged during hundreds of modulation periods with Ω mod /2π = 10 kHz. The switching decay times is fitted with exponential fits. b) For V mod = 900, 980 and 1200 mV, from top to bottom: reference modulation signal, transmitted signal V out and demodulated noise amplitude at the hot (cold) state mechanical frequency Ω hot m (Ω cold m ). Ω mod /2π = 10 Hz. c) Probability for the optical state to be in the hot (blue) or cold (green) state.

  Fig.4: a) We evaluate the mode energy |a| 2 while the laser wavelength λ L = 2πc/ω L is swept forward (red) and backward (black), using P in = 10 mW, giving rise to a ∼ 1 nm wide bistability (red area). α mod = α hf = 0. b) Residence probability as a function of the modulation depth for λ L = 1566.30 nm and Ω mod /2π = 500 Hz.

Fig. 5 :

 5 Fig.5: Signal-to-noise ratio (SNR) plotted as a function of the high-frequency signal amplitude, α hf , for α mod = 0.005 (black) and α mod = 10 -5 (gray). Ω mod /2π = 500 Hz and Ω hf /2π = 10 kHz.

Table 1 :

 1 List of parameters used in the numerics.

	Physical parameter	value
	λ 0 : resonance wavelength	1566.3 nm
	κ i : internal loss rate	28.2 GHz
	κ e : external loss rate	18.4 GHz
	n 0 : refractive index	3.16
	dn dT : thermo-optic coef.	1.9298•10 -4 K -1
	K t : Thermalization rate	125 kHz
	κ th : thermal resistance	1.62 K.fJ -1

Acknowledgement This work is supported by the French RENATECH network, the European Union's Horizon 2020 research innovation program under grant agreement No 732894 (FET Proactive HOT), the Agence Nationale de la Recherche as part of the "Investissements d'Avenir" program (Labex NanoSaclay, ANR-10-LABX-0035) with the flagship project CON-DOR and the JCJC project ADOR (ANR-19-CE24-0011-01).

Supporting information Theoretical model

Here we derive a model to describe the vibrational resonance amplification of a weak signal in a thermo-optic waveguide-coupled optical mode. We consider a single optical mode with complex amplitude a. We write ω c , κ i , κ e , s i and ω L respectively the mode frequency, its intrinsic loss rate, its external loss rate, the laser amplitude and the laser frequency. The system can be modelled within coupled mode theory