Mark Geoghegan
email: mark.geoghegan@newcastle.ac.uk

Marta M Mróz
email: marta.mroz@polimi.it

Chiara Botta
email: chiara.botta@scitec.cnr.it

Laurie Parrenin
email: laurie.parrenin@gmail.com

Cyril Brochon
email: cyril.brochon@enscbp.fr

Eric Cloutet
email: eric.cloutet@enscbp.fr

Eleni Pavlopoulou
email: epavlopoulou@iesl.forth.gr

Georges Hadziioannou

Tersilla Virgili
email: tersilla.virgili@polimi.it

Optical Gain in Semiconducting Polymer Nano and Mesoparticles

Keywords: transient absorption spectroscopy, conjugated polymers, stimulated emission, nanoparticles, mesoparticles

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

Introduction

Semiconducting polymer films may be a viable alternative to inorganic materials when high performance is less important than cost. Polymers have a significant advantage in that they may be solution processed, which permits film preparation routes such as spin coating, doctor blading, and drop casting. However, polymers need not be dissolved in solution but rather can be dispersed in a liquid medium as latex particles before being cast as films. This is appropriate for certain synthetic routes, such as emulsion polymerization or nanoprecipitation, which result in a particle suspension [START_REF] Pecher | Nanoparticles of Conjugated Polymers[END_REF]. These are all compatible with various cross-coupling routes that can be used to synthesize conjugated polymers [START_REF] Geoghegan | Polymer Electronics[END_REF][START_REF] Murad | Conducting Polymers for Optoelectronic Devices and Organic Solar Cells: A Review[END_REF]. This wide choice of synthesis and film preparation routes allows fine tuning of the morphological and structural properties of the film. Semiconducting polymer particles can be optimized for light emission [START_REF] Gao | Site Isolation of Emitters within Cross-Linked Polymer Nanoparticles for White Electroluminescence[END_REF][START_REF] Piok | Organic Light-Emitting Devices Fabricated from Semiconducting Nanospheres[END_REF], energy generation [START_REF] Kietzke | Novel approaches to polymer blends based on polymer nanoparticles[END_REF][START_REF] Ganzer | Ultrafast spectroscopy on water-processable PCBM: rod-coil block copolymer nanoparticles[END_REF], or phototherapies [START_REF] Yang | Convertible Organic Nanoparticles for Near-Infrared Photothermal Ablation of Cancer Cells[END_REF][START_REF] Zangoli | characterization, and biological applications of semiconducting polythiophene-based nanoparticles[END_REF], among other purposes [START_REF] Jana | Photon Harvesting in Conjugated Polymer-Based Functional Nanoparticles[END_REF][START_REF] Xu | Nanoparticles made of π-conjugated compounds targeted for chemical and biological applications[END_REF][START_REF] Macfarlane | Functional nanoparticles through πconjugated polymer self-assembly[END_REF].

The optoelectronic properties of polymer films depend, at least indirectly, on their morphology as well as their electronic structure. Because optical and electronic properties are inextricably linked, electronic and excitonic transport need to be optimized for the required optical properties. For main-chain conjugated polymers, charge transport benefits from both interchain and intrachain phenomena and also from a good degree of crystallization. Optical properties, however, are less affected by crystallization than transport properties but perturbing the conjugation length of even rigid polymers is known to affect them [START_REF] Kim | Optoelectronic and Charge Transport Properties at Organic-Organic Semiconductor Interfaces: Comparison between Polyfluorene-Based Polymer Blend and Copolymer[END_REF]. Indeed, where comparisons between amorphous and crystalline polymers have been made, it is noted that changes in optical properties are due to the concomitant change in conjugation length as opposed to chain ordering [START_REF] Cadby | Film morphology and photophysics of polyfluorene[END_REF][START_REF] Khan | Morphology-dependent energy transfer within polyfluorene thin films[END_REF].

Semiconducting polymer films that display stimulated emission often have excellent electroluminescent properties, and so effort spent optimizing optical gain is often performed with wider applications in mind [START_REF] Jiang | Organic solid-state lasers: a materials view and future development[END_REF]. Although optical gain in semiconducting polymer films has been known for some time [START_REF] Moses | High quantum efficiency luminescence from a conducting polymer in solution: A novel polymer laser dye[END_REF] and is not rare, it is not commonplace either, and some thought is needed to achieve it. Interchain interactions are often a means of quenching excitations, so control of the photophysics of conjugated polymers can be obtained by enhancing single chain properties through dilution in a matrix of an amorphous polymer [START_REF] Virgili | Ultrafast resonant optical switching in isolated polyfluorenes chains[END_REF], encapsulation in nanochannels [START_REF] Martini | Controlling optical gain in semiconducting polymers with nanoscale chain positioning and alignment[END_REF], or through the use of an appropriately structured film [START_REF] Portone | Tailoring optical properties and stimulated emission in nanostructured polythiophene[END_REF]. The separation of semiconducting polymers from each other (dilution) does not allow for good electronic transport, but other work [START_REF] Yap | Simultaneous optimization of charge-carrier mobility and optical gain in semiconducting polymer films[END_REF] showed that, by mixing different polyfluorenes, it was possible to have films of semiconducting polymers with both optical gain and excellent charge transport properties. Further, alternative, approaches include mimicking the effects of dilution by adding a small amount of dimethyl phenylene into polyfluorene chains [START_REF] Zhang | Simultaneously Enhancing Photoluminescence Quantum Efficiency and Optical Gain of Polyfluorene via Backbone Intercalation of 2,5-Dimethyl-1,4-Phenylene[END_REF] or by the supramolecular separation of chains by threading cyclodextrin rings around the conjugated polymer to create polyrotaxanes [START_REF] Mroz | Amplified Spontaneous Emission in Conjugated Polyrotaxanes Under Quasi-cw Pumping[END_REF].

Films are usually produced by spin coating or doctor blading from a polymer solution, which give rise to good uniformity and control of thickness. Rough films may result in weaker optical properties due to the scattering of light by the film surface. Such scattering may be exacerbated in particle films, particularly for particles of optical dimensions. Nevertheless, it could be reasonably argued that nanoparticle films convey some advantages because confinement inherent in the particles allows some element of dilution by enhancing intrachain properties due to the extended conjugation length with respect to interchain interactions. However, transient absorption spectroscopy revealed no evidence of stimulated emission in an aqueous nanoparticle suspension of poly[N-9"-heptadecanyl-2,7-carbazole-alt-5,5-(4,7-di-2-thienyl-2′,1′,3′-benzothiadiazole)] (PCDTBT) [START_REF] Virgili | Size-Dependent Photophysical Behavior of Low Bandgap Semiconducting Polymer Particles[END_REF]. In contrast, a mesoparticle suspension of PCDTBT in propanol did show stimulated emission in the infrared region of the spectrum. Here, transient absorption is used to show that the photophysics of PCDTBT mesoparticle films (of diameter ~450 nm) is very different to those of nanoparticle films (~50 nm), and in particular, the mesoparticles exhibit good optical gain. Nanoparticles on the other hand exhibit photoinduced absorption due to the formation of charged states, rather than stimulated emission. Both polymers used for these films exhibit the same optical properties when spin cast from chloroform solution to form a continuous PCDTBT film. PCDTBT films are largely amorphous when spin coated from chlorobenzene solution [START_REF] Cho | A Thermally Stable Semiconducting Polymer[END_REF], although there is some evidence of limited order when cast from chloroform [START_REF] Aïch | Electrical and Thermoelectric Properties of Poly(2,7-Carbazole) Derivatives[END_REF]. These results show that it is possible to create good optical gain by the preparation route used to synthesize the films, rather than by using elaborate means to dilute semiconducting polymers.

Results

Figure 1 shows the absorption and emission spectra of the spin-coated nanoparticle, mesoparticle, and continuous polymer films taken as reference. The absorption spectrum of the continuous polymer film (red line) exhibits two broad transitions with peaks at 390 (3.19 eV) and at 550 nm (2.26 eV). In agreement with earlier work [START_REF] Virgili | Size-Dependent Photophysical Behavior of Low Bandgap Semiconducting Polymer Particles[END_REF][START_REF] Banerji | Exciton Formation, Relaxation, and Decay in PCDTBT[END_REF][START_REF] Etzold | Ultrafast Exciton Dissociation Followed by Nongeminate Charge Recombination in PCDTBT:PCBM Photovoltaic Blends[END_REF] the two bands are identified as the π-π*-transition of the first and second excited singlet states (S1 and S2). The emission spectrum shows a broad and unstructured emission between 620 (2 eV) and 920 nm (1.35 eV) (red line), peaking at 703 nm (1.77eV). The spectra of the nanoparticles (blue lines) are very similar to that of the continuous polymer film, with a red shift of around 65 meV for both peaks in both absorption and emission. This is likely to indicate a more planar chain conformation and an increased conjugation length [START_REF] Brown | Effect of interchain interactions on the absorption and emission of poly(3-hexylthiophene)[END_REF][START_REF] Eggimann | How β-phase content moderates chain cenergy transfer in polyfluorene films[END_REF]. The absorption spectrum of the mesoparticle film (black line in Figure 1), prepared from a lower molar mass polymer, is quite different from those previously reported for PCDTBT [START_REF] Banerji | Exciton Formation, Relaxation, and Decay in PCDTBT[END_REF][START_REF] Etzold | Ultrafast Exciton Dissociation Followed by Nongeminate Charge Recombination in PCDTBT:PCBM Photovoltaic Blends[END_REF]. In contrast to what would be expected [START_REF] Shao | Probing the molecular weight dependent intramolecular interactions in single molecules of PCDTBT[END_REF], both transition peaks are redshifted by around 150 meV, and also their relative intensity is changed, with the transition at 607 nm being more intense than that at 415 nm. The long wavelength absorption tail is likely to be due to the Rayleigh scattering of light from the 450 nm particles [START_REF] Mcquarrie | Statistical Mechanics[END_REF]. Despite the red-shifted absorption, the photoluminescence (PL) emission (black line in Figure 1) is similar to the one of the polymeric film with the main peak at 703 nm, resulting in an evident reduction of the Stokes shift.

The left panel of Figure 2 presents the pump-probe spectra at different probe delays for the continuous polymer film after excitation at 400 nm (Figure 2a) and at 610 nm (Figure 2c). After excitation at 400 nm the pump-probe spectrum shows the presence of an initial bleaching band centred around 550 nm (low energy peak in the absorption spectrum) and a photoinduced absorption band between 650 nm and 780 nm attributed to charged states [START_REF] Virgili | Size-Dependent Photophysical Behavior of Low Bandgap Semiconducting Polymer Particles[END_REF]. In the first ~10 ps a red shift of the bleaching is present (Figure 2a, red line) indicating an energy migration towards the lower energetic sites of the absorption spectrum. In this spectral region an overlap between the positive bleaching signal and the negative one due to charge formation is also present. To highlight the charge formation, the dynamics at 750 nm are considered (Figure 2b), which shows a formation time of around 1 ps.

The low energetic sites at 610 nm were also excited. The pump-probe spectrum shows an instantaneous bleaching of the high and low energetic excited states (see red and black line in Figure 2d) inducing an initial enlargement of the bleaching band of around 60 nm (0.19 eV). Nevertheless, even with this pump excitation the photoinduced absorption (PIA) band due to charged states is present with a low intensity. It is therefore possible to create charges even for excitation of the tail of the polymeric film of the absorption spectrum. In this way fewer charges are created, but a gain region in the continuous polymer film is nevertheless not observed. The pump-probe spectra are presented at different probe delays for the nanoparticle film after excitation at 400 nm (Figure 3a) and at 610 nm (Figure 3b) and for the mesoparticle film (Figure 3c,d). The difference between the films is evident.

The nanoparticle film is first considered. After excitation at 400 nm, the pump-probe spectrum shows the presence of an initial bleaching band centred at around 570 nm (low energy peak in the absorption spectrum) and a photoinduced absorption band at 700-800 nm attributed elsewhere to charged states [START_REF] Virgili | Size-Dependent Photophysical Behavior of Low Bandgap Semiconducting Polymer Particles[END_REF]. The spectrum is similar to the one of the polymer film even if the dynamics show an instantaneous charge formation (blue line in Figure 4a) indicating a more efficient process. Even in this case there is a small red shift of the bleaching band toward low energetic states.

After excitation at 610 nm, a very different behavior was observed compared to that after 400 nm excitation. The positive band extended up to 750 nm indicating the presence of a long-lived gain region that is not present when exciting at 400 nm. It is therefore possible to have optical gain after excitation on the tail of the absorption spectrum. The direct excitation at 610 nm creates few charges (formation time ~10 ps, see Figure 4c), allowing the nanoparticles to have gain (see Figure 5a for the schematic energy diagram after each excitations). Excitation of the mesoparticle films at 400 nm produces a large PIA band in the visible region and a positive band in the near-infrared region. As already seen in solution [START_REF] Virgili | Size-Dependent Photophysical Behavior of Low Bandgap Semiconducting Polymer Particles[END_REF], the positive band is attributed to stimulated emission, while the instantaneous negative band present in the visible region was not present in the transient transmitted spectra of the solution. Looking at the decay traces (Figure 4b), the signal at 550 nm is instantaneously created. It has a rapid initial decay and then a recovery of the signal after ~200 fs. The positive signal at 750 nm is not instantaneously created and has a similar formation time of ~200 fs. The initial PIA band (PIA1) at 550 nm is attributed to the Sn-Sm singlet transition, it is followed by a fast internal conversion that increases the population of the state S1, which in turn produces the stimulated emission (gain) at 750 nm and a signal (PIA2) from S1 to Sn (see the schematic energy diagram in Figure 5b) in the visible region. Moreover, after ~50 ps, a negative band centred at 650 nm (PIA3), attributed to the formation of charged states, overlaps with the initial photoinduced signal [START_REF] Virgili | Size-Dependent Photophysical Behavior of Low Bandgap Semiconducting Polymer Particles[END_REF]. In the mesoparticle film, after excitation at 400 nm, this represents a different photophysics compared to the continuous and nanoparticle films: an internal transition from Sn to S1, a strong gain signal in the near-infrared region, and slow charge formation are evident.

After excitation of the mesoparticle film at 610 nm, the pump-probe spectra are similar to the 400 nm excitation. A large negative band in the visible region and a positive one (gain) in the near infrared are present. The decay traces show that the signals at 750 and 550 nm are both instantaneously created and behave similarly (Figure 4d). This is then a direct excitation of the S1 state which produces a stimulated emission in the near infrared and PIA2 in the visible region (S1-Sn transition). In this case, charges are not created. In the mesoparticle films, it is possible to have optical gain after excitation at both excitation wavelengths. After excitation at 400 nm, a 200 fs internal conversion brings the excitation to the lower energetic site from which stimulated emission is observed. The direct excitation at 610 nm allows an instantaneous emission with no charge formation.

Discussion

Before comparing the different films, it is important to recall that when these two particles are dissolved in a common solvent, and spin cast to form a continuous film, they exhibit the same photophysical properties [START_REF] Virgili | Size-Dependent Photophysical Behavior of Low Bandgap Semiconducting Polymer Particles[END_REF]. This observation provides strong evidence that differences between the optical properties of the mesoparticle and nanoparticle films are due to size effects, rather than their differing molar masses.

The red shift in the absorption spectrum of the nanoparticle films, with respect to the continuous polymer films, can be associated to an increased conjugation length. The Stokes shift of ~0.49 eV (~150 nm), however, is similar for the two films, with the absorption and emission spectra red-shifted by the same amount. The mesoparticle films, on the other hand, display a much smaller (0.27 eV) Stokes shift, leaving the emission spectrum in the same position as that for the uniform film. The differences in the steady-state absorption-emission spectra of polymeric assemblies are commonly analyzed within an H/Jaggregate model [START_REF] Baghgar | Effect of Polymer Chain Folding on the Transition from Hto J-Aggregate Behavior in P3HT Nanofibers[END_REF][START_REF] Spano | The Spectral Signatures of Frenkel Polarons in H-and J-Aggregates[END_REF]. As for the mesoparticle suspension [START_REF] Virgili | Size-Dependent Photophysical Behavior of Low Bandgap Semiconducting Polymer Particles[END_REF] but also for the mesoparticle film absorption spectrum, the π-π*-transition of the first singlet state S1 is red-shifted with respect to the nanoparticle one while the PL peak is blue-shifted as expected for Jaggregate-like behavior, which is generally associated with stronger intrachain than interchain coupling [START_REF] Portone | Tailoring optical properties and stimulated emission in nanostructured polythiophene[END_REF][START_REF] Spano | The Spectral Signatures of Frenkel Polarons in H-and J-Aggregates[END_REF]. The weak interchain interactions of the mesoparticle films are responsible for the limited charge formation observed here [START_REF] Banerji | Exciton Formation, Relaxation, and Decay in PCDTBT[END_REF].

Stronger interchain interactions evident in the continuous and nanoparticle films explains the pump-probe spectra. The fast red shift in the bleaching, due to excitonic energy transfer [START_REF] Virgili | Size-Dependent Photophysical Behavior of Low Bandgap Semiconducting Polymer Particles[END_REF], reveals the presence of molecular disorder within the continuous polymer and nanoparticle films. This red shift is not present in the mesoparticle films. due to reduced exciton migration and/or greater molecular order in the mesoparticles.

Furthermore, charge formation in the mesoparticle films is suppressed at the longer (610 nm) excitation wavelength. Charge formation is very rapid for both the nanoparticle and continuous films. However, at the higher energy excitation, charge formation is observed in all films, although it is slowest in the mesoparticle film, where formation takes place on the order of 50 ps, compared with 1 ps (continuous film) and 200 fs (nanoparticle film). The reduced efficiency of charge formation in the mesoparticle films indicates of a more ordered system in these films, whereas the rapid charge formation in the other films (nanoparticle and continuous) supports the conclusion of a disordered molecular structure.

Optical gain is absent from the continuous polymer films at both excitation wavelengths, and also from the nanoparticle films at 400 nm. However, long-lived stimulated emission is observed following excitation at 610 nm for the nanoparticle films. This longer wavelength was not tested in the suspension. The mesoparticle films exhibit optical gain in the infrared at both wavelengths.

We cannot exclude optical gain existing in all films, but the nanoparticle and continuous polymer films all exhibit strong photoinduced absorption, which may overlap with the stimulated emission signal. Certainly, it has been possible to detect stimulated emission in the infrared in PCDTBT films spin-cast from chloroform, after excitation at 400 nm, and this signal strongly overlapped with that due to photobleaching [START_REF] Etzold | Ultrafast Exciton Dissociation Followed by Nongeminate Charge Recombination in PCDTBT:PCBM Photovoltaic Blends[END_REF].

It is unlikely that the particle size is the direct reason for the observed behavior, but rather the effect of the particle size on the chain conformation within the particles. This may also be a consequence of the synthetic route used to prepare the films. In the mesoparticle films, the chains are able to order in a way that gives rise to fewer interparticle interactions competing with intraparticle interactions and consequently less charge formation, which allows for an increased optical gain.

The phenomenological explanation of the optical properties of these films also applies to those of the nanoparticle and mesoparticle suspensions, where it was concluded that the mesoparticles have weaker interchain coupling and the optical properties, including stimulated emission, are likely to be due to a more J-aggregate-like behavior that is not observed for the nanoparticle suspension. Although the nanoparticle films revealed some stimulated emission when fewer charges were created after excitation at 610 nm, the longer wavelength excitation was not performed in the experiments on the suspensions, precluding a direct comparison.

Materials and Methods

Two routes to create PCDTBT particles were used. High molar mass nanoparticles were created as described previously [START_REF] Parrenin | Photoactive Donor-Acceptor Composite Nanoparticles Dispersed in Water[END_REF]. PCDTBT (number average molar mass, Mn = 20.2 kDa with dispersity 2.2) was first synthesized by a Suzuki cross-coupling route [START_REF] Wakim | Highly efficient organic solar cells based on a poly(2,7-carbazole) derivative[END_REF] before being dissolved in chloroform. This was added to an aqueous solution of sodium dodecyl sulfate and sonicated to form a mini emulsion. The chloroform evaporated after heating at 70 °C, leaving an aqueous nanoparticle dispersion. The mean particle diameter was determined by transmission electron microscopy to be 50 nm.

PCDTBT mesoparticles were created as previously described [START_REF] Parrenin | Low Bandgap Semiconducting Copolymer Nanoparticles by Suzuki Cross-Coupling Polymerization in Alcoholic Dispersed Media[END_REF]. PCDTBT (Mn = 4.5 kDa with dispersity 2.1) was synthesized by Suzuki cross-coupling polymerization in propanol solution with poly(vinyl pyrrolidone) added as a surfactant. By adjusting the quantity of surfactant, it was possible to adjust the particle size in the range 0.33-1.3 μm, with a virtually uniform size distribution. The films created in this work were made of 0.45 ± 0.05 μm particles.

For the continuous film, chloroform solutions (40 mg/mL) were prepared by dissolving the mesoparticles and then a film of the polymer was prepared by spin coating.

Film thickness was not a determining factor in these experiments but, in order to get a good optical signal during the transient absorption measurements, the mesoparticle films needed to be less than 1 μm thick. Nanoparticle and continuous films were substantially thinner, and they can be assumed to be of the order of ~100 nm that is typical for transient absorption measurements.

Absorption and emission spectra were acquired using a Shimadzu UV-3600 spectrophotometer (Shimadzu, Marne la-Vallée, France) and a Horiba Scientific Fluoromax-4 spectrofluorometer (Horiba, Palaissau, France), respectively. The excitation wavelength for the emission spectra was 390 nm.

Time-resolved measurements were performed using a homebuilt femtosecond pump-probe setup [START_REF] Virgili | Size-Dependent Photophysical Behavior of Low Bandgap Semiconducting Polymer Particles[END_REF][START_REF] Virgili | Direct Evidence of Torsional Motion in an Aggregation-Induced Emissive Chromophore[END_REF]. A Ti:sapphire regenerative amplifier (Libra, Coherent, CA, USA) was used as a laser source, delivering 100 fs pulses at a central wavelength of 800 nm with 4 mJ pulse energy at a repetition rate of 1 kHz. The second harmonic of the fundamental wavelength was used for the excitation at λ = 400 nm and for the excitation pulse at 610 nm, a single stage non-linear optical parametric amplifier, pumped at 400 nm was used. In order to minimize bimolecular effects, the excitation density was kept very low. White light generated with a 2 mm-thick sapphire plate was used as a probe in the visiblenear infrared range from 450 to 780 nm. For spectrally resolved detection of the probe light, spectrographs and CCD arrays were used. The chirp in the white light pulse was taken into account during the analysis and evaluation of the two-dimensional (wavelength and time) ΔT(λ,τ)/T maps before extraction of the spectral and temporal data using homemade software. Overall, a temporal resolution of at least 150 fs was achieved for all excitation wavelengths.

Conclusions

Stimulated emission is a phenomenon which occurs in many conjugated polymers where photoluminescence is observed. However, more often than not, it competes with other processes such as charge creation or photobleaching and its detection and consequent utility is impeded. To this extent, researchers have gone to considerable lengths to create films which have limited interchain optical behavior. In this work, UV-visible absorption and photoluminescence spectroscopy, coupled with time-resolved transient absorption spectroscopy was used to demonstrate unambiguous optical gain in hard latex (450 nm) films of the conjugated polymer PCDTBT. Optical gain was weaker in nanoparticle films. It is concluded that conjugated polymer mini emulsions may be a viable route to creating systems with increased conjugation length. However, a full understanding of why such a preparation route increases the polymer conjugation length, or even whether it does or not, is lacking.

Figure 1 .

 1 Figure 1. Absorption and emission (rightmost peaks) spectra of the continuous polymer (red line), mesoparticle (black line), and nanoparticle (blue line) spin coated films.

Figure 2 .

 2 Figure 2. Left panel: Transient absorption spectra at different probe delays for the continuous polymer film at the different wavelengths of pump excitation: 400 nm (a) and 610 nm (b). The arrows correspond to different wavelengths (black, 550; red, 650; and blue 750 nm) at which the temporal decays are shown in (b) and (d) Right panel: Time decays at the two different pump excitation wavelengths: (b) 400 and (d) 610 nm. The shape of the 610 nm pump beam is also shown in (c).

Figure 3 .

 3 Figure 3. Transient absorption spectra at different probe delays for nanoparticle, (a) and (b), and mesoparticle, (c) and (d), films at two different wavelengths of pump excitation: 400 nm, (a) and (c) and 610 nm, (b) and (d). The ordinates for all axes are ΔT/T and the abscissae are wavelength.

Figure 4 .

 4 Figure 4. Time decays at different wavelengths for nanoparticle, (a) and (c), and for mesoparticle, (b) and (d), spin coated films at two different pump excitation wavelengths of 400 nm, (a) and (b), and 610 nm, (c) and (d), are displayed. The ordinates and abscissae for all axes are normalized ΔT/T and probe delay respectively.

Figure 5 .

 5 Figure 5. Schematic energy diagrams for the nanoparticle (a) and mesoparticle (b) polymer films, where CS indicates the charged state and IC is an internal transition.

Acknowledgments: M.G. is grateful to the Politecnico di Milano and CNR-SCITEC for hosting him in Milan.

Data Availability Statement: Data are contained within the article.

Author Contributions: Conceptualization, M.G. and T.V.; methodology, T.V. and E.C.; formal analysis, T.V.; investigation, T.V., M.M.M., C.B. (Chiara Botta), L.P., E.P., C.B. (Cyril Brochon), and E.C.; writing-original draft preparation, M.G. and T.V.; writing-review and editing, all authors; visualization, M.G. and T.V.; supervision, T.V., M.G., E.C., and E.P.; project administration, M.G.; funding acquisition, T.V. and G.H. All authors have read and agreed to the published version of the manuscript.

Funding: T.V. and M.M.M. acknowledge the project TIMES from Regione Lombardia. In Bordeaux, L.P., C.B., E.C., E.P., and G.H. acknowledge the financial support from the ADEME Project ISOCEL number 1182C0212 and from Region Nouvelle Aquitaine. This work was performed within the framework of the LCPO/Arkema/ANR Industrial Chair HOMERIC ANR-13-CHIN-0002-01 with grant number AC-2013-365.

Conflicts of Interest:

The authors declare no conflict of interest.


Lorène Julia Marchand
email: lorene.marchand@univ-rennes1.fr

Inge Dox

Jožica Gričar

Peter Prislan

Jan Van Den Bulcke

Patrick Fonti

Matteo Campioli

Timing of spring xylogenesis in temperate deciduous tree species relates to tree growth characteristics and previous autumn phenology

Keywords: Fagus sylvatica, Quercus robur, Betula pendula, cambium, seasonal growth, hardwood species

We explored the timing of spring xylogenesis and its potential drivers in homogeneous mature forest stands in a temperate European region. Three species with contrasting leaf development dynamics and wood anatomy were studied: European beech, silver birch and pedunculate oak.

Detailed phenological observations of xylogenesis and leaf phenology were performed from summer 2017 till spring 2018. Cambium reactivation (CR) occurred before the buds of oak and birch were swollen, whereas these two phenological phases were concurrent for beech.

On the other hand, initial earlywood vessels were fully differentiated (FDIEV) after leaf unfolding for all three species. Timing of CR was correlated to average ring-width of the last 10 years , tree diameter, and, partially, with tree age. In addition, the timing of FDIEV was correlated to tree age and previous" year autumn phenology i.e. timing of wood growth cessation and onset of leaf senescence. Multivariate models could explain up to 68% of the variability of CR and 55% of the variability of FDIEV. In addition to the "species" factor, the variability could be explained by ca. 30% by tree characteristics and previous" years autumn phenology for both CR and FDIEV. These findings are important to better identify which factors (other than environment) can be driving the onset of the growing season and highlight the influence of tree growth characteristics and previous" year phenology on spring wood phenology, wood formation and, potentially, forest production.

Introduction

Phenology is the study of event"s occurrence in plant and animal"s life cycles in interaction with the environment [START_REF] Forrest | Toward a synthetic understanding of the role of phenology in ecology and evolution[END_REF]. For trees, phenological observations are typically focussed on primary growth, for example, on the leaf phenophases [START_REF] Koch | Guidelines for plant phenological observations[END_REF]. In temperate deciduous forests, the leaf life cycle starts in spring with bud-burst and ends in fall with leaf senescence. For secondary growth (indicated also as wood formation or xylogenesis), phenology is much less studied because of the complexity of its assessment [START_REF] Zuidema | A Wood Biology Agenda to Support Global Vegetation Modelling[END_REF][START_REF] Friend | On the need to consider wood formation processes in global vegetation models and a suggested approach[END_REF]. However, thanks to its seasonal pattern, secondary growth can also provide key information about tree functioning. In fact, the timing of the increasing number of cambial cells i.e. cambium reactivation (CR from here on), and of the full differentiation of the first row of earlywood vessels i.e. fully differentiated initial earlywood vessels (FDIEV from here on), are key events to characterize the timing of reactivation of the main carbon sink in the tree and the formation of its new hydraulic system [START_REF] Begum | Regulation of cambial activity in relation to environmental conditions: understanding the role of temperature in wood formation of trees[END_REF]. For ring-porous species, CR preeceds bud-burst [START_REF] González-González | Comparative cambial dynamics and phenology of Quercus robur L. and Q. pyrenaica Willd. in an Atlantic forest of the northwestern Iberian Peninsula[END_REF], whereas for diffuse-porous species the two events are concurrent [START_REF] Čufar | Tree-ring variation, wood formation and phenology of beech (Fagus sylvatica) from a representative site in Slovenia, SE Central Europe[END_REF][START_REF] Prislan | Phenological variation in xylem and phloem formation in Fagus sylvatica from two contrasting sites[END_REF]. However, the relationships between spring xylogenesis and the other leaf phenophases (e.g. swollen buds, complete leaf unfolding) are less investigated. Accurate studies on CR and FDIEV (secondary growth phenology) and leaf development (primary growth phenology) are essential to better understand the drivers and mechanisms of the onset of the growing season and important tree (xylem) functional traits [START_REF] Rathgeber | Biological Basis of Tree-Ring Formation: A Crash Course[END_REF].

Secondary growth represents one-third of the annual biomass production in temperate forest ecosystems [START_REF] Campioli | Evaluating the convergence between eddy-covariance and biometric methods for assessing carbon budgets of forests[END_REF]. Nevertheless, models of annual biomass production are largely based on primary growth variables [START_REF] Friend | On the need to consider wood formation processes in global vegetation models and a suggested approach[END_REF]. Therefore, the identification of determinants of secondary growth is important to predict vegetation growth and to increase A c c e p t e d m a n u s c r i p t the accuracy of dynamic global vegetation models [START_REF] Zuidema | A Wood Biology Agenda to Support Global Vegetation Modelling[END_REF][START_REF] Friend | On the need to consider wood formation processes in global vegetation models and a suggested approach[END_REF].

For the main tree species in temperate European forests (e.g. European beech, pedunculate oak, sessile oak and Scots pine) secondary growth is influenced more by the growth duration (controlled by spring and autumn phenology) than the growth rate [START_REF] Michelot | Comparing the intra-annual wood formation of three European species (Fagus sylvatica, Quercus petraea and Pinus sylvestris) as related to leaf phenology and non-structural carbohydrate dynamics[END_REF][START_REF] Pérez-De-Lis | Environmental conditions and vascular cambium regulate carbon allocation to xylem growth in deciduous oaks[END_REF]. Therefore it is important to gain more insights on CR and FDIEV.

For coniferous species, the onset of primary and secondary growth across the northern hemisphere is driven mainly by photoperiod and temperature [START_REF] Piao | Plant phenology and global climate change: Current progresses and challenges[END_REF][START_REF] Huang | Photoperiod and temperature as dominant environmental drivers triggering secondary growth resumption in Northern Hemisphere conifers[END_REF]. Concerning temperature, chilling (amount of cold experienced in late autumn-early winter) and forcing (amount of warmth experienced in late winter-early spring) are both needed to explain the beginning of xylogenesis [START_REF] Delpierre | Chilling and forcing temperatures interact to predict the onset of wood formation in Northern Hemisphere conifers[END_REF]. For angiosperm species, spring xylogenesis is also driven by climate, especially warm temperatures in late winter-early spring [START_REF] Begum | Regulation of cambial activity in relation to environmental conditions: understanding the role of temperature in wood formation of trees[END_REF], 2018, Delpierre, et al. 2016). However, even with this knowledge, inter-individual variability of xylogenesis within a homogenous forest stand cannot be explained [START_REF] Prislan | Phenological variation in xylem and phloem formation in Fagus sylvatica from two contrasting sites[END_REF]. Comprehensive phenological studies of xylogenesis are therefore needed to understand tree functioning at individual scale within a stand.

At individual level, the timing of spring xylogenesis could be related to tree characteristics (e.g. age, size) or previous" year autumn phenology of leaves and wood [START_REF] Rossi | Age-dependent xylogenesis in timberline conifers[END_REF], Pérez-de-Lis et al. 2016). Multiple studies on coniferous species show that age and size affect spring xylogenesis. For example, younger individuals of Larix decidua, Pinus cembra and Picea abies have CR earlier than mature ones [START_REF] Rossi | Age-dependent xylogenesis in timberline conifers[END_REF][START_REF] Rossi | Causes and correlations in cambium phenology: towards an integrated framework of xylogenesis[END_REF]. Furthermore, taller silver-firs have CR earlier the smaller individuals [START_REF] Rathgeber | Cambial activity related to tree size in a mature silver-fir plantation[END_REF]. However, less is known about the drivers of intra-stand variability of spring xylogenesis for angiosperm species as many of the existing studies on spring xylogenesis for angiosperms are mainly focused on hydraulic properties and wood anatomy (e.g. Pérez- de-Lis et al. 2016). Therefore, we hypothesize that, in a homogeneous mature stand, younger trees present earlier CR and A c c e p t e d m a n u s c r i p t FDIEV than older trees. We expect for diameter the same trend as for age because of the consistent relationship between the two variables in homogenous stands. Timing of leaf autumn phenology in the previous year might affect the timing of current wood growth throught the content of non-structural carbohydrates. In fact, later senescence likely assures more carbon uptake and more storage of reserves (Bréda et al. 2006, Estiarte and[START_REF] Estiarte | Alteration of the phenology of leaf senescence and fall in winter deciduous species by climate change: effects on nutrient proficiency[END_REF] in a period (autumn) when growth is modest [START_REF] Dox | Timeline of autumn phenology in temperate deciduous trees[END_REF]. More reserves might trigger fast spring reactivation because of greater fitness and the potential to mitigate risks such as, for example, late cold spells [START_REF] Morin | Variation in cold hardiness and carbohydrate concentration from dormancy induction to bud burst among provenances of three European oak species[END_REF][START_REF] Hartmann | Understanding the roles of nonstructural carbohydrates in forest trees -from what we can measure to what we want to know[END_REF][START_REF] Tixier | Non-structural Carbohydrates in Dormant Woody Perennials; The Tale of Winter Survival and Spring Arrival[END_REF]. In fact, the negative relationship between carbon reserves and spring phenology (i.e. the larger the carbon reserves, the earlier the spring reactivation) has been already suggested for primary growth [START_REF] Delpierre | Tree phenological ranks repeat from year to year and correlate with growth in temperate deciduous forests[END_REF][START_REF] Pérez-De-Lis | Environmental conditions and vascular cambium regulate carbon allocation to xylem growth in deciduous oaks[END_REF]).

Therefore we hypothesize that trees with later onset of leaf senescence also present earlier CR and FDIEV. Our hypothesis is further strengthened by the correlation between carbon storage in the previous year and width-ring variation in the following year, as observed for beech [START_REF] Skomarkova | Inter-annual and seasonal variability of radial growth, wood density and carbon isotope ratios in tree rings of beech (Fagus sylvatica) growing in Germany and Italy[END_REF]. As in temperate angiosperms, autumn wood growth cessation is generally concurrent (or shortly before) onset of leaf senescence [START_REF] Dox | Timeline of autumn phenology in temperate deciduous trees[END_REF], we expect spring xylogenesis to be influenced in the same way by both the onset of leaf senescence and wood growth cessation. To our knowledge, the interaction between previous autumn phenology (leaf senescence and wood growth cessation) and the timing of spring xylogenesis (CR and FDIEV) has never been studied for angiosperm or coniferous tree species.

Here we postulate that under favorable climatic conditions (thus without stressors or extreme events), tree growth characteristics and previous" year autumn phenology can influence the timing of spring wood growth in addition to wood anatomy. The objectives of our study are (i) to determine the timing of CR and FDIEV in relation to spring leaf phenology for common deciduous forest tree species of the European temperate zone, (ii) study the timing of CR and A c c e p t e d m a n u s c r i p t FDIEV in relation to tree growth characteristics and previous autumn phenological events, and (iii) build statistical models to determine the importance of such variables on the timing of CR and FDIEV. The study was performed on pedunculate oak (Quercus robur L.), European beech (Fagus sylvatica L.) and silver birch (Betula pendula Roth.), in two forests nearby Antwerp, Belgium, between summer 2017 and spring 2018. Whereas, (though with fewer details), previous studies have already compared leaf-and wood spring phenology of beech and oak, to the best of our knowledge, this is also the first study presenting a detailed description of leaf-and wood spring phenology for birch.

Material and methods

Species and sites

The species selected for this study (oak, beech and birch) are among the most common temperate forest tree species in Europe. Oak and beech are late-successional species with a determinate growth (i.e. they produce leaves in 1-2 flushes) whereas birch is a pioneer species with an indeterminate growth (i.e. it produces leaves continuously throughout the season).

Regarding the wood type, beech and birch are diffuse-porous species, whereas oak is ringporous. The study took place on two nearby (< 10 km) mature temperate forest sites in northern Belgium, specifically the park of Brasschaat (PB, 51°12" N, 4°26"E and 18m a.s.l) and the "Klein Schietveld" (KS, 51°21"N, 4°37"E and 22m a.s.l). The region is characterized by a maritime climate with regular precipitation (778 mm year -1 ) and a monthly mean temperature ranging from +3°C (January) to 17°C (July) [START_REF] Campioli | Current status and predicted impact of climate change on forest production and biogeochemistry in the temperate oceanic European zone: review and prospects for Belgium as a case study[END_REF]. Both sites are characterized by sandy soil, but with a thicker organic layer (60-90cm) at PB compared to KS.

On each site, we selected one homogeneous stand of each species. These were monospecific A c c e p t e d m a n u s c r i p t for beech PB, oak PB, beech KS and birch KS, but with the presence of other species for birch PB (mainly Quercus rubra) and oak KS (mainly Betula pendula).

Tree selection and their characteristics

Dominant or codominant healthy trees of each species were selected. In total we selected 28 trees (4 oaks, 4 beeches and 4 birches) at each site, except for KS, where we selected 8 beeches trees. For each tree we determined the stem diameter at breast height and a tree density index to consider the effect of competition (Supp. Data Table I). The latter was calculated considering the total basal area of all (co)dominant trees within a circle centered at our target tree. Additionally, to assess the age and the average growth rate of the last decade, selected trees were cored on the north and south part with a Pressler corer. Tree-ring width measurements were performed with a Lintab measuring stage (0.01 precision) and TSAP-Win software (Table Supp. Data. 1).

Wood phenology

Wood xylogenesis was assessed from late August till late November 2017 and from early March till mid-June 2018. Weekly to biweekly stem micro-cores were collected using a Trephor [START_REF] Rossi | Trephor: A New Tool for Sampling Microcores from tree stems[END_REF]). To avoid any impact of previous sampling, sampling followed an upward spiral at a stem height between 1.3 m to 2 m. Sample processing to obtain microsections was performed at the Slovenian Forestry Institute (Ljubljana, Slovenia) according to the protocol described in [START_REF] Prislan | Wood Sample Preparation for Microscopic Analysis[END_REF]. In short, immediately after removal from the trees, the samples were stored in 70 % of ethanol solution and further processed in the laboratory for preparation of permanent cross-sections (5-10 µm in thickness) stained with safranin and Astra blue and embedded in Euparal. Histometric observations and image capturing were performed using the LAS image analysis software (LEICA Application Suite version 4.9.0, copyright [2003][2004][2005][2006] on images of the microsections collected with a LEICA A c c e p t e d m a n u s c r i p t DMC 4500 camera connected to a LEICA DM 4000B light microscope (Leica Microsystems, Wetzlar, Germany) using both transmission and polarised light. The observation included five xylem phenophases [START_REF] Puchałka | Does tree-ring formation follow leaf phenology in Pedunculate oak (Quercus robur L.)?[END_REF], Gričar et al. 2017), specifically: 1: Dormant cambium (DC, cambium is non-productive); 2: Cambium reactivation (CR, increase of the number of cambial cells; production of new xylem and phloem cells); 3: Onset of xylem growth (OG, enlargement of newly formed xylem cells); 4: Appearance of initial earlywood vessels (IEV, enlargement of the first row of vessels of the current year-ring), and 5: Fully differentiated initial earlywood vessels (FDIEV, vessels and neighboring fibers are lignified, sap can be conducted through those vessels). During autumn, two additional phenophases were considered: 6: xylem cells in the wall-thickening phase and 7: last formed cells are mature except for very rare solitary terminal cells. Wood growth cessation dates were determined when latewood differentiation was completed which was visible thanks to the red colour of the cells under Safranin-Astra blue staining [START_REF] Dox | Timeline of autumn phenology in temperate deciduous trees[END_REF].

Leaf phenology

Leaf phenology was assessed in autumn 2017 and spring 2018. Autumn phenology was monitored visually, following the method of [START_REF] Vitasse | Responses of canopy duration to temperature changes in four temperate tree species: relative contributions of spring and autumn leaf phenology[END_REF], from late August till late November, considering first only the tree canopy coloration and, later in the season, a combination of canopy coloration and leaf fall (decrease in canopy greenness). The date of onset of leaf senescence was determined at the appearance of a sharp decrease in canopy greenness [START_REF] Mariën | Detecting the onset of autumn leaf senescence in deciduous forest trees of the temperate zone[END_REF].

Spring phenology was monitored weekly or twice a week from early March to late May 2018 on 10 randomly selected buds, equally distributed from the upper and lower part of the crown, to account for variation within the crown [START_REF] Davi | Diversity of leaf unfolding dynamics among tree species: New insights from a study along an altitudinal gradient[END_REF]. The considered phenophases included bud swelling, bud-burst, leave unfolding and were determined when 50% of the

Statistical analyses

Comparisons of timing between wood and leaf phenophases were performed for each stand and species with a permuted two-sample T-test (999 permutations; n= 4 except for KS-beech where n=8). Relationships of CR and FDIEV along with tree characteristics and previous autumn phenology were quantified using univariate linear models. When the necessary assumptions of normality of the residuals (tested with Shapiro test, [START_REF] Royston | Algorithm AS181: the W test for Normality[END_REF]) and homogeneity of variance (tested with ncvTest, package car, Fox and Weisberg 2019) were not met, even after a box-cox transformation (using package MASS, Venable and Ripley 2002), a generalized linear model with quasi-poisson distribution was performed. The relationships between CR and FDIEV and categorical variables (site and species) were evaluated with oneway ANOVA, type III (package car). To consider nonlinear relationships as well, the Random Forest method (package randomForest, [START_REF] Liaw | Classification and Regression by randomForest[END_REF] was applied to rank the potential explanatory variables. Random Forest is a partitioning method that produces a large ensemble of regression trees considering the full dataset but with random subset of predictor variables. We set 5000 trees to have stable results and 3 (for FDIEV) or 4 variables (for CR) for each split, following the floor function (package randomForest). Finally, we used multivariate linear models (that can be applied with both continuous and categorical variables) to determine the most relevant explanatory variables and to build multivariate models. Because of the relatively small dataset, these models were built starting with the three most important variables identified with the univariate models. Interactions between variables were checked and were reported if significant. Tests for normality of residuals and variance homogeneity were done as reported above. Because these conditions were not met for FDIEV, instead of a linear model, we applied, in this case, a generalized linear model, with quasi-

Results

Wood and leaf phenology

CR (mean date per stand) took place between DOY 77 -87 for birch, DOY 108 -109 for beech and on DOY 87 for oak (Fig. 1). FDIEV took place on DOY 162 for birch, DOY 160 -162 for beech and DOY 134 -143 for oak (Fig. 1). Both CR and FDIEV showed significant inter-species variation (p < 0.01) but not significant inter-site variation (Fig. 1). Leaf phenology (mean date per stand) started with the "swollen bud" stage on DOY 96 -107, then bud-burst occurred on DOY 101 -110 and finally, leaf unfolding took place on DOY 106 -112 (Fig. 2). CR occurred 10 -20 days before buds were swollen for birch and oak (Fig. 2, Table I), but after that (around budburst) for beech (Fig. 2, Table I). FDIEV occurred around a month or two after leaf unfolding (stand average between 22 -55 days; p < 0.05; Table I).

Concerning previous autumn phenology, wood growth cessation varied between species and sites and onset of leaf senescence varied between species only (Fig. 1). Both events were concurrent (p > 0.10) except for KS-beech (p = 0.02) where an early wood growth cessation was observed.

Univariate relationships

Cambium reactivation date was significantly related to species and also to 10 years average ring-width, with a later CR in case of narrower ring-widths (Fig. 3a, Table II). Furthermore, a weak positive relation (p = 0.060) was found between CR and tree diameter (Fig. 3c, Table II). The Random Forest analysis confirmed that species was the most important factor and that tree diameter had also an important influence (Table III). Opposed to the linear models, the

A c c e p t e d m a n u s c r i p t

Random Forest analysis ranked "age" as a key variable and 10 years growth rate among the least important variables (Table III). Age was not significantly linearly related to CR date mainly because of the old trees in the oak-PB stand (Fig. 3e). Indeed, the average age of this stand was 115 years against an average age of 60 years for the rest of the studied trees.

Without these 4 old trees the relationship between CR date and age is significant (T-value = 5.119; p-value < 0.01; Fig. 3e), with later CR for older trees until a certain age. The fact that the 10 years growth rate was ranked low by the Random Forest analysis only means that other variables (e.g. growth cessation and onset of senescence) had stronger non-linear relationships with CR (Table III).

FDIEV was highly related to species and also to age and previous autumn phenology (i.e. onset of leaf senescence and wood growth cessation) (Fig. 3b,d,f, Table II). Age was significantly related to FDIEV when all trees were considered (Fig. 3f, Table II). On the other hand, when trees older than 75 years were removed, the relation became non-significant (Fig. 3f). The relationships between FDIEV and onset of leaf senescence and wood growth cessation, were negative, meaning that a later FDIEV was found for trees with an early leaf onset of senescence or previous" year early cessation of wood growth (Fig. 3b,d Table II).

Moreover, a weak negative linear relation (p = 0.063) was found between FDIEV and tree diameter (Table II). For FDIEV, Random Forest confirmed that species, age, onset of leaf senescence and wood growth cessation were the most important variables affecting FDIEV even when non-linear relationships were also considered (Table III).

Statistical models of wood spring phenology

For both spring wood phenological events, CR and FDIEV, a linear model was developed with the three most important variables obtained with univariate model analysis (Table IV andSI Table III; correlations between variables are shown in SI Table II). For CR, a multivariate model including "species", "diameter" and "age" explained 68% of the variability (Table IV).

A c c e p t e d m a n u s c r i p t

A model without the "species" factor could still explain 28% of the variability (Table IV). For FDIEV, a model including "species", "age" and "onset of leaf senescence" explained 55% of the variability (Table IV). Also in this case, a linear model without the species factor could explain ca. 30% of the variability (Table IV). An alternative model with "wood growth cessation" instead of "onset of leaf senescence", provided similar results (SI Table III).

Discussion

Timeline of wood and leaf spring phenology

For oak, our study aligns with previous observations of xylogenesis and bud-burst for this species and other comparable oak species [START_REF] Sass-Klaassen | Vessel formation in relation to leaf phenology in pedunculate oak and European ash[END_REF][START_REF] González-González | Comparative cambial dynamics and phenology of Quercus robur L. and Q. pyrenaica Willd. in an Atlantic forest of the northwestern Iberian Peninsula[END_REF][START_REF] Pérez-De-Lis | Do changes in spring phenology affect earlywood vessels? Perspective from the xylogenesis monitoring of two sympatric ring-porous oaks[END_REF][START_REF] Pérez-De-Lis | Environmental conditions and vascular cambium regulate carbon allocation to xylem growth in deciduous oaks[END_REF], Lavrič et al. 2017, Gričar et al. 2017). For example, the timing of CR was similar to the one reported for pedunculate oak in The Netherlands [START_REF] Sass-Klaassen | Vessel formation in relation to leaf phenology in pedunculate oak and European ash[END_REF]. On the other hand, we could also show that CR occurs not only before bud-burst but also before the buds are swollen (the first stage of leaf spring phenology). This detail has not been reported before. In beech, CR occurred almost at the same time as bud-burst. Our data are in agreement with other studies on beech [START_REF] Čufar | Tree-ring variation, wood formation and phenology of beech (Fagus sylvatica) from a representative site in Slovenia, SE Central Europe[END_REF][START_REF] Prislan | Phenological variation in xylem and phloem formation in Fagus sylvatica from two contrasting sites[END_REF]. Furthermore, our observation that earlywood vessels were fully formed in beech at the beginning of June is comparable to data from other sites in Germany (in the Bavarian Alps) and Slovenia [START_REF] Čufar | Tree-ring variation, wood formation and phenology of beech (Fagus sylvatica) from a representative site in Slovenia, SE Central Europe[END_REF](Čufar et al. , Kraus et al. 2016)). For birch, although diffuse porous as beech, CR already occured in the second half of March, significantly before bud swelling. We speculate that for this pioneer species the cambial reactivation, uncoupled to leaf phenology, could be due to the temperature conditions in the spring of 2018 with an unusually cooler March that delayed bud-burst (see temperature and phenology comparison with 2019, when bud-burst occurred slightly before CR, SI Fig. 1).

A c c e p t e d m a n u s c r i p t

Only few (and less detailed) studies on spring xylogenesis and wood phenology for birch exist. Therefore, such lag between CR and bud-burst for a pioneer diffuse-porous species is not supported by other research so far. It might indicate that, at least for the birch population investigated here, spring cold spells tend to block the leaf unfolding rather than the cambial activity.

Factors influencing spring wood phenology

In addition to the factor "species", CR was related to the average ring-width of the last 10 years, with an earlier start of the growing season for trees with larger ring widths. A positive relationship between vessel density and ring width was already found for oak and beech (Delpierre et al. 2016[START_REF] Prislan | Precipitation is not limiting for xylem formation dynamics and vessel development in European beech from two temperate forest sites[END_REF], suggesting that larger annual ring width is associated with an earlier start of the wood growth. Cambium reactivation was also related to tree age (at least till the age of 75 years), with younger trees starting CR earlier than older trees. Although the age range was limited (40-75 years), this follows the principle that young trees need to compete for light against the older larger trees as it was also reported for various conifers [START_REF] Rossi | Age-dependent xylogenesis in timberline conifers[END_REF], 2008[START_REF] Li | Age dependence of xylogenesis and its climatic sensitivity in Smith fir on the south-eastern Tibetan Plateau[END_REF]. Therefore, the first hypothesis is confirmed for CR but only for trees up to 75 years old. The fact that beech trees of 45 and 130 years in the same stands were reported to have the same CR [START_REF] Giagli | Monitoring of seasonal dynamics in two age-different european beech stands[END_REF], might confirm that the positive relationship between CR and age is not relevant for older trees. Average ring width is more closely related than age to CR in homogenous stands. The relationship between tree-ring width and CR supports the idea that trees with better fitness and growth are the ones with earlier start of the secondary growth.

FDIEV was not related to age for trees up to 75 years old, but it occurred earlier for the oldest trees (oaks) considered in our study. The first hypothesis is therefore not valid for FDIEV that shows for old trees even the opposite pattern than expected. However, the findings indirectly agrees with Pérez-de-Lis et al. ( 2016), who found a negative relationship between FDIEV and A c c e p t e d m a n u s c r i p t tree diameter for pedunculate oaks. The latter can be an indication that older, larger, trees (thus with more resources and reserves) can have a fast spring reactivation. We also found negative relationships between FDIEV and onset of leaf senescence and wood growth cessation in the previous autumn. To our knowledge, this is the first study to relate earlywood vessel differentiation with previous autumn phenology. The second hypothesis (i.e. that later autumn phenology might trigger earlier spring wood xylogenesis) is therefore confirmed for FDIEV but not for CR. Maybe the surplus of reserves gained by later onset of senescence is more important when the reserve storage is close to its minimum in June (the timing of FDIEV) rather than in March-April (the timing of CR) when the carbon reserves are still rather full [START_REF] Barbaroux | Contrasting distribution and seasonal dynamics of carbohydrate reserves in stem wood of adult ring-porous sessile oak and diffuse-porous beech trees[END_REF].

Statistical models of wood spring phenology

Our models, which are based on simple correlations, captured about 70% and 55% of the variability of CR and FDIEV, respectively. Letting the species factor aside, about 30% of the variability was explained by growth characteristics and previous"s year autumn phenology.

This indicates that the difficult task to fully understand the inter-tree variability in spring xylogenesis still needs to be investigated and more potential drivers (e.g. non-structural carbon of reserves, nutrient content) need to be taken into account. In any case, our study is one of the few trying to elucidate the inter-individual variability in xylogenesis in homogenous stands of angiosperms. The model results, even without being mechanistic, are statistically significant and help us to better understand the tree functioning. For example, we showed that CR and FDIEV are not influenced by the same variables. Cambium reactivation is mainly related to growth characteristics (age, 10 years average ring-width and tree diameter), whereas FDIEV is also affected by previous autumn dynamics (onset of leaf senescence and wood cessation). Thus, FDIEV could be more sensitive to annual variations compared to CR. This is important as FDIEV can be a predictor of latewood formation A c c e p t e d m a n u s c r i p t (Pérez-de-Lis, García-González, et al. 2016). As annual wood growth and forest ecosystem services are depending on the timing of both CR and FDIEV [START_REF] Rathgeber | Biological Basis of Tree-Ring Formation: A Crash Course[END_REF], their relationships with tree growth characteristics and previous autumn phenology should not be overlooked.

Conclusions

Our results demonstrate that CR can occur before any signs of spring bud activity (swollen buds) for birch and oak, while CR and the swollen bud phase occur concurrently for beech. Initial earlywood vessels were fully differentiated always after leaf unfolding. We show that the timing of CR correlates to tree diameter, ring-width and age. On the other hand, the timing of FDIEV was significantly correlated to age and previous" year wood growth cessation and onset of leaf senescence. The hypothesis that younger, smaller, trees present earlier CR and FDIEV has been proved true only for CR of trees up to 75 years old. The hypothesis that trees with later autumn phenology present earlier CR and FDIEV has been proved true only for FDIEV but not for CR. Our study might offer new perspectives in thinning strategies in very homogenous stands. In fact, in absence of major dendrometric differences, trees with earlier leaf senescence might be preferentially removed as they appear to have a later start of the wood growing season. However, our study is mainly fundamental, and the novel insights on inter-individual variability of xylogenesis open the door for further research. The observed relationship between previous year" autumn phenology and spring xylogenesis is particularly interesting as it indicates that stressors affecting autumn phenology during the current growing season (e.g. heat waves, severe drought, etc., all expected to increase in frequencies and intensity in the future) can have an impact on next year wood phenology. A c c e p t e d m a n u s c r i p t

  A c c e p t e d m a n u s c r i p t observed buds or leaves have reached the corresponding phenophase. For detail see Marchand et al. (2020).

  s c r i p t poisson distribution. ANOVA, type III analyses, were done on each model to test categorical variables.

Figure 1 :

 1 Figure 1: Timing of observed phenophases (DOY ± standard error) grouped by species and site. OLS = Onset of leaf senescence in 2017 (grey), CWG = cessation wood growth in 2017 (orange), CR = cambium reactivation in 2018 (yellow), FDIEV = fully differentiated initial earlywood vessels in 2018 (green). P-values are the result of the ANOVA test for "site" and "species" effect.

Figure 2 :

 2 Figure 2: Progress of 2018 spring cambium activity (black dots; stand mean DOY ± standard error) in comparison to the timing of the observed leaf phenophases (coloured bands, mean DOY ± standard error) grouped by species and site. Blue = bud swelling, green = bud-burst, orange = leaf unfolding, DC = dormant cambium, CR = cambium reactivation, OG = onset of xylem growth, IEV = Appearance of initial earlywood vessels, FDIEV = fully differentiated initial earlywood vessels A c c e p t e d m a n u s c r i p t

Figure 3 :

 3 Figure 3: Right panel: Relevant relationships between day of year (DOY) of 2018 cambium reactivation (CR, orange points) with dendrometric parameters or 2017 autumn phenophases: (a) 10 years average ring-widths, (c) tree diameter and (e) tree age. Left panel: relevant relationships between DOY of 2018 fully differentiated initial earlywood vessel (FDIEV, blue points) with dendrometric parameters or 2017 autumn phenophases: (b) wood growth cessation in previous autumn, (d) onset of senescence in previous autumn, and (f) tree age. Each point corresponds to a tree (n = 28). Black lines and grey bands represent the linear models with the 95% confidence interval. Model p-values are reported on each panel. For panel (e) and (f) the red line (and red p-value) represents the linear model only considering trees < 75 years old.

Table I :

 I The p-value of the student test with 999 permutations comparing the 2018 DOY of leaf-and wood spring phenological events for each species within each stand (see Figure1).

	Bold values indicate whether the difference is significant (p-value <0.05)
	site	specie	CR and swollen	CR and bud-burst	FDIEV and leaf
			bud		unfolding
		Beech	0.296	0.808	0.024
	Park of Brasshaat (PB)	Birch	0.04	0.032	0.03
		Oak	0.022	0.032	0.036
	Klein Schietveld (KS)	Beech Birch Oak	0.03 0.03 0.022 m a n u s c r i p t 0.194 0.002 0.03 0.04 0.026 0.032
	A c c e p t e d

Table II :

 II Characteristics of univariate linear models between 2018 cambium reactivation (CR) and fully differentiated initial earlywood vessels (FDIEV) with tree dendrometric parameters or 2017 autumn phenophases for all monitored trees at both sites (see Figure3). Oak (n=8), beech (n=12) and birch (n=8). Significant models (p<0.05) are in bold; "*" indicates the use of a generalized linear model quasi-poisson instead of a linear model. The "¤" indicates the results of a one-way ANOVA.

	Potential predictors Onset of senescence Wood cessation Age 10 years average ring-width Tree diameter Tree density species site	CR T-value or F value (¤) P-value 0.531 0.6 -0.426 0.674 0.358 0.723 -2,198 0.037 1.966 0.060 0.235 0.816 30.296¤ <0.01 2.269¤ 0.144 m a n u s c r i p t FDIEV T-value or F-value (¤) P-value -2.388* 0.024 -2.927* 0.007 -2.81* 0.009 -0.199* 0.844 -1.942* 0.063 -0.984* 0.334 32.744¤ <0.01 0.239*¤ 0.625
	A c c e p t e d

Table III :

 III Importance of predictors (indicated with an increase of mean square error IncMSE) obtained with Random Forest analysis for cambium reactivation (CR) and fully differentiated initial earlywood vessels (FDIEV) in 2018 for all observed trees (n = 28).

	Potential predictors	CR	% IncMSE	FDIEV	% IncMSE
	species		133.61		61.94
	Age Tree diameter Cambium cessation 2017 Onset of senescence 2017 Tree density index 10 years average ring-width Sites		30.44 19.81 8.48 8.43 2.73 1.95 0.75 m a n u s c r i p t 4.22 -0.96 6.26 4.47 2.17 -3.24 -0.41
	A c c e p t e d

Table IV :

 IV Characteristics of multivariate multi-species linear models of cambium reactivation (CR) and multivariate multi-species generalized linear models for fully differentiated initial earlywood vessels (FDIEV); models are presented with or without the "specie" factor (CR n=28; FDIEV n=27). na: not applicable.

		CR model with age		FDIEV model with onset of
						senescence
		F-	P-	F-	P-	LR	P-	T-	P-
		value	value	value	value	Chisq	value	value	value
	species	23.06 <0.01 na	na	14.97	<0.01 na	na
	Tree diameter	3.19	0.09	11.24 0.03	na	na	na	na
	Age	1.32	0.26	4.67	0.04	0.684	0.41 -2.12	0.04
	Tree diameter : Age na	na	6.61	0.017 na	na	na	na
	Onset of senescence na	na	na	na	0.001	0.97 -1.58	0.13
	R² adjusted	68%		28%		55%		29%

Acknowledgments

We acknowledge the funding support from ERC Starting Grant LEAF-FALL (714916). We thank also the Agency for Forest and Nature of the Flemish Government (ANB), the Military Defence of Belgium (Defensie), and the City of Brasschaat, that granted permission to conduct research in the study areas. We are also grateful to the Flemish Institute for Nature and Forest (INBO) and the Integrated Carbon Observation System (ICOS) for providing meteorological data. Special thanks are due to Stijn Willen (UGent) for analysis of the tree cores and to Dirk Leyssens (ANB) for the support.


Jiafa Zeng

Zhujun Yu

Majda Mekic

Jiangping Liu

Sheng Li

Gwendal Loisel

Wei Gao

Adrien Gandolfo

Zhen Zhou

Xinming Wang

Hartmut Herrmann

Sasho Gligorovski
email: gligorovski@gig.ac.cn

Xue Li

intensities of CO 2 , and estimated OH radical concentrations

Cooking emissions represent a major source of air pollution in the indoor environment and exhibit adverse health effects caused by particulate matter together with volatile organic compounds (VOCs). A multitude of unknown compounds are released during cooking, some of which play important roles as precursors of more hazardous secondary organic aerosols in indoor air. Here, we applied secondary electrospray ionization highresolution mass spectrometry for real-time measurements of VOCs and particles from cooking peanut oil in the presence of 300 ppbv nitrogen oxides (NO x ) generated by a gas stove in an indoor environment. More than 600 compounds have been found during and after cooking, including N-heterocyclic compounds, Oheterocyclic compounds, aldehydes, fatty acids, and oxidation products. Approximately 200 compounds appeared after cooking and were hence secondarily formed products. The most abundant compound was 9-oxononanoic acid (C 9 H 16 O 3 ), which is likely the product formed during the heterogeneous hydroxyl (OH) radical oxidation of oleic acid (C 18 H 34 O 2 ) or linoleic acid (C 18 H 32 O 2 ). Real-time detection of an important number of organic compounds in indoor air poses a challenge to indoor air quality and models, which do not account for this extremely large range of compounds.

INTRODUCTION

Combustion processes such as gas stove cooking represent some of the most important temporal sources of air pollution in the indoor environment. A large number of organic particles and volatile organic compounds (VOCs) are released, 1 many of which exhibit adverse health effects. [START_REF] Gligorovski | Indoor (Photo)Chemistry in China and Resulting Health Effects[END_REF][START_REF] Wang | Role of Chinese Cooking Emissions on Ambient Air Quality and Human Health[END_REF][START_REF] Huang | Characteristics and Health Impacts of VOCs and Carbonyls Associated with Residential Cooking Activities in Hong Kong[END_REF][START_REF] Kim | Home Kitchen Ventilation, Cooking Fuels, and Lung Cancer Risk in a Prospective Cohort of Never Smoking Women in Shanghai, China[END_REF] Cooking fumes contain a wide range of compounds, such as alkanes, alkenes, aldehydes, and fatty acids. [START_REF] Schauer | Measurement of Emissions from Air Pollution Sources. 4. C1-C27 Organic Compounds from Cooking with Seed Oils[END_REF][START_REF] Klein | Characterization of Gas-Phase Organics Using Proton Transfer Reaction Time-of-Flight Mass Spectrometry: Cooking Emissions[END_REF][START_REF] Cheng | Characterization of Volatile Organic Compounds from Different Cooking Emissions[END_REF][START_REF] Liu | Emission of Volatile Organic Compounds and Production of Secondary Organic Aerosol from Stir-Frying Spices[END_REF][START_REF] Zhang | Speciation of VOCs in the Cooking Fumes from Five Edible Oils and Their Corresponding Health Risk Assessments[END_REF] Cooking also releases large amounts of nitrogen oxides (NO x ) and nitrous acid (HONO) into the indoor air. [START_REF] Bartolomei | Combustion Processes as a Source of High Levels of Indoor Hydroxyl Radicals through the Photolysis of Nitrous Acid[END_REF][START_REF] Zhou | Time-Resolved Measurements of Nitric Oxide, Nitrogen Dioxide, and Nitrous Acid in an Occupied New York Home[END_REF][START_REF] Collins | Evidence for Gas-Surface Equilibrium Control of Indoor Nitrous Acid[END_REF][START_REF] Liu | Assessing Indoor Gas Phase Oxidation Capacity through Real-Time Measurements of HONO and NO x in Guangzhou, China[END_REF] Recently, it has been demonstrated that HONO can be easily photodissociated in indoor air by the ultraviolet (UV) fraction of sunlight filtered through the windows [START_REF] Gandolfo | Spectrally Resolved Actinic Flux and Photolysis Frequencies of Key Species within an Indoor Environment[END_REF][START_REF] Kowal | Wavelength-Resolved Photon Fluxes of Indoor Light Sources: Implications for HOx Production[END_REF] and generate hydroxyl (OH) radicals at the level of ∼10 6 cm -3 . [START_REF] Liu | Assessing Indoor Gas Phase Oxidation Capacity through Real-Time Measurements of HONO and NO x in Guangzhou, China[END_REF][START_REF] Gomez Alvarez | Unexpectedly High Indoor Hydroxyl Radical Concentrations Associated with Nitrous Acid[END_REF] The compounds in cooking fumes can react with highly reactive OH radicals, forming secondary organic compounds. This is assumed to be similar to the formation of secondary products by reactions of OH radicals with VOCs and organic particles in the outdoor atmosphere, which is known and has been experimentally verified. [START_REF] Emmerson | Urban Atmospheric Chemistry during the PUMA Campaign 2: Radical Budgets for OH, HO2 and RO2[END_REF][START_REF] Ren | Behavior of OH and HO2 in the Winter Atmosphere in New York City[END_REF][START_REF] Elshorbany | Seasonal Dependence of the Oxidation Capacity of the City of Santiago de Chile[END_REF][START_REF] Berndt | Hydroxyl Radical-Induced Formation of Highly Oxidized Organic Compounds[END_REF] Some of these secondary products have lower volatility and are more toxic than the initially emitted VOCs. [START_REF] Gligorovski | The Oxidative Capacity of Indoor Atmospheres[END_REF] Real-time measurements of compounds formed upon reaction of OH radicals with primary emissions from cooking have been attracting more attention in recent years, [START_REF] Farmer | Overview of HOMEChem: House Observations of Microbial and Environmental Chemistry[END_REF][START_REF] Liu | Significant Production of Secondary Organic Aerosol from Emissions of Heated Cooking Oils[END_REF] but only a few studies have been reported to date. [START_REF] Lu | The Emission of PM 2.5 in Respiratory Zone from Chinese Family Cooking and Its Health Effect[END_REF] This can be partially explained by the complexity of these reaction processes and also by the limitations in the analytical tools.

Secondary electrospray ionization high-resolution mass spectrometry (SESI-HRMS) is a novel ambient MS technique and shows great potential in capturing the evolution of compounds in real time. [START_REF] Tejero Rioseras | Comprehensive Real-Time Analysis of the Yeast Volatilome[END_REF] The application of SESI especially favors the detection of reactive compounds, such as fatty acids, aldehydes, ketones, etc.; 27 furthermore, accurate identification of compounds can be realized by HRMS. Recently, via application of SESI-HRMS to real-time breath analysis, around 1000 breath compounds per subject have been measured and unambiguously identified. [START_REF] Gaugg | Expanding Metabolite Coverage of Real-Time Breath Analysis by Coupling a Universal Secondary Electrospray Ionization Source and High Resolution Mass Spectrometry  a Pilot Study on Tobacco Smokers[END_REF] In this study, SESI-HRMS was applied for the first time to measure the primary emissions from cooking fumes and the real-time evolutions of primary and secondary organic compounds in an indoor environment. The literature suggests that high levels of OH radicals can be produced upon Evolution of Indoor Cooking Emissions Captured by Using Secondary Electrospray Ionization High-Resolution Mass Spectrometry photodissociation of HONO, which determine the oxidation capacity of indoor air. [START_REF] Liu | Assessing Indoor Gas Phase Oxidation Capacity through Real-Time Measurements of HONO and NO x in Guangzhou, China[END_REF] Thus, the ultimate goal of this study is to identify the secondarily formed products upon reaction of OH radicals with the primarily emitted organic compounds in the gas phase and particulate phase in the indoor air.

MATERIALS AND METHODS

2.1. Experimental Setup. Measurements were conducted in a 23 m 3 room with a window (0.99 m 2 ) facing southeast at our institute (Figure S1). The sunlight irradiated the room through the window during certain periods of the day. Peanut oil was heated in a pan on a natural gas stove (ZB-16M, Iwatani) for ∼10 min, and then the gas stove was turned off. Two fans were placed at two corners of the room (Figure S1) to ensure the indoor air was homogeneously mixed during the experiments. An air exchange rate of ∼0.6 h -1 was estimated from the evolution of the carbon dioxide (CO 2 ) from combustion (details in the Supporting Information) applying the procedure described elsewhere. [START_REF] Liu | Assessing Indoor Gas Phase Oxidation Capacity through Real-Time Measurements of HONO and NO x in Guangzhou, China[END_REF] The photostationary state (PSS) approach [START_REF] Gomez Alvarez | Unexpectedly High Indoor Hydroxyl Radical Concentrations Associated with Nitrous Acid[END_REF][START_REF] Elshorbany | Oxidation Capacity of the City Air of Santiago, Chile[END_REF] was used to estimate the time profiles of OH radical values. Details about the estimations of OH radical values can refer to our previous work14, including the measurements of spectral irradiance, NO x , HONO, and ozone (O 3 ). The measured O 3 values during the campaign were below 200 parts per trillion by volume (pptv), and therefore, ozonolysis of alkenes as an OH source can be neglected. A control experiment was performed under the same condition but without cooking; each set of measurements was repeated at least twice.

2.2. SESI-HRMS Analysis of Emissions. A homemade SESI source was coupled to a high-resolution quadrupole Orbitrap mass spectrometer (Thermo Scientific). [START_REF] Li | Real-Time Breath Analysis by Using Secondary Nanoelectrospray Ionization Coupled to High Resolution Mass Spectrometry[END_REF] Cooking emissions were delivered through a 6.7 m long Teflon tube (4 mm inside diameter) at a rate of 1.4 L min -1 . Samples without pretreatment were introduced into the SESI-HRMS instrument in real time for direct ionization (Figure S1). [START_REF] Elshorbany | Oxidation Capacity of the City Air of Santiago, Chile[END_REF] Compounds were detected at a time resolution of 1 s in fast polarity-switching mode, and ESI voltages were set to 2.5 and -2.5 kV. The scan range is m/z 50-500, and the mass resolution is 140000 at m/z 255. The limit of detection (LOD) and sensitivity of SESI-HRMS, sampling efficiency, and delay of the 6.7 m sampling line were evaluated by using a set of gas standards (details in Text S1 and Table S1). The LODs range from pptv to parts per brillion by volume (ppbv), which agree with the previous studies. [START_REF] Wolf | Direct Quantification of Chemical Warfare Agents and Related Compounds at Low Ppt Levels: Comparing Active Capillary Dielectric Barrier Discharge Plasma Ionization and Secondary Electrospray Ionization Mass Spectrometry[END_REF] The sampling efficiencies are 4-28% and 14-49% for 2 and 20 ppbv standard compounds, respectively. The sampling efficiencies vary greatly from compound to compound, and thus we compared the signal intensities only within individual compounds instead of between compounds. Because <200 pptv O 3 was detected in this indoor environment, O 3 -induced heterogeneous chemistry in the Teflon tubing can be neglected. Prior to measurements, the mass spectrometer was calibrated with the commercial standard calibration solutions (Thermo Scientific; details in Text S1).

2.3. Data Analysis. The raw data from SESI-HRMS were pretreated with a Python program named BreathFinder to obtain a data matrix, including m/z values, signal intensities, and time points of all of the ions detected (details in the Supporting Information). Then, the data matrix was further investigated by using cluster analysis in Matlab R2018b and Kendrick Mass Defect (KMD) analysis. [START_REF] Hughey | Kendrick Mass Defect Spectrum: A Compact Visual Analysis for Ultrahigh-Resolution Broadband Mass Spectra[END_REF] KMD analysis is especially suitable for the data acquired by HRMS to show visualized trends in the elemental compositions of complex organic compounds. [START_REF] Kourtchev | Effects of Anthropogenic Emissions on the Molecular Composition of Urban Organic Aerosols: An Ultrahigh Resolution Mass Spectrometry Study[END_REF] The Kendrick mass (KM) and KMD are calculated as follows:

KM IUPAC mass (14/14.01565) = ×

(1)

KMD nominal KM exact KM = - (2)
where 14.01565 is the exact mass of the CH 2 unit and 14 is the nominal mass. Thus, compounds that have the same general molecular formula (e.g., C n H 2n O 2 ) but differ in the number of the base units (such as CH 2 groups) have identical KMD values and thus lie on the same horizontal line in the plot of KM versus KMD. The elemental composition of each compoupnd used in KMD analysis was obtained using Xcalibur (version 4.2.28.14, Thermo).

The enhancement ratio (ER) was used to evaluate the enhancement of the signal intensity of a compound during cooking; ER is calculated on the basis of the maximum signal intensity of a compound and its signal intensity at the beginning of the cooking experiment.

RESULTS AND DISCUSSION

3.1. Evolution of Cooking Emissions. By using SESI-HRMS, an important number of reported compounds have been detected, such as aldehydes, fatty acids, N-heterocyclic compounds, and O-heterocyclic compounds (Figure S2). Obvious differences were found between mass spectral fingerprints before, during, and after cooking (Figure S3). To determine how many groups of compounds were emitted from cooking, 864 ions with intensities of ≥1 × 10 5 au captured in real time were further investigated using cluster analysis. Six groups (657 ions in total) characterized with different chronological variations are presented in Figure 1, which have not been observed in the control experiment (Figure S5). In addition, to explore the evolution of primary emissions and the formation of secondarily formed products, molecular formulas were assigned to the ions detected (Tables S2-S7), and ER values were provided for some of these compounds (Table S8). For example, compounds in group 1 (93 ions) are characterized by reaching the peak level during cooking. Most of the N-heterocyclic compounds, e.g., trimethylpyrazine (ER = 3.0) and 3-methoxy-4-pyridine (ER = 3.2), were identified in this group, which are consistent with the flavor components with low boiling points in peanut oil during cooking. [START_REF] Liu | Changes in Volatile Compounds of Peanut Oil during the Roasting Process for Production of Aromatic Roasted Peanut Oil[END_REF][START_REF] Hu | Characterization of Volatile Components in Four Vegetable Oils by Headspace Two-Dimensional Comprehensive Chromatography Time-of-Fl Ight Mass Spectrometry[END_REF] Compounds in groups 2 (97 ions) and 3 (208 ions) achieved the maximum signal intensity around the time when the gas stove was turned off. These two groups mainly contain unsaturated aldehydes and O-heterocyclic compounds that are formed by chemical processes involved in the oil heating, inducing the breaking bonds of the fatty acids. [START_REF] Huang | Characteristics and Health Impacts of VOCs and Carbonyls Associated with Residential Cooking Activities in Hong Kong[END_REF] Polyunsaturated fatty acids such as linoleic acid (LA) can readily undergo autoxidation [START_REF] Gardner | Oxygen Radical Chemistry of Polyunsaturated Fatty Acids[END_REF] and generate the unsaturated aldehyde decadienal in peanut oil. [START_REF] Choe | Comprehensive Reviews in Food Science and Food Safety Mechanisms and Factors for Edible Oil Oxidation[END_REF] 2,4-Decadienal (ER = 155.1), which exhibits the second highest signal intensity, is genotoxic to human lung cells and increases the carcinogenic risk of inhabitants. [START_REF] Wang | Pulmonary Changes Induced by Trans, Trans -2,4-Decadienal, a Component of Cooking Oil Fumes[END_REF][START_REF] Young | DNA Damage Induced by Trans, Trans-2,4-Decadienal (Tt-DDE), a Component of Cooking Oil Fume, in Human Bronchial Epithelial Cells[END_REF][START_REF] Mu | Indoor Air Pollution and Risk of Lung Cancer among Chinese Female Non-Smokers[END_REF] O-Heterocyclic compounds like 2-acetyl-5methylfuran (ER = 4.9) and 5-propyldihydro-2(3H)-furanone (ER = 5.4) are flavor compounds found during thermal processing of foods [START_REF] Vranova | Furan in Fooda Review[END_REF] and thus can be used as tracers for cooking emission. [START_REF] Rogge | Sources of Fine Organic Aerosol. 1. Charbroilers and Meat Cooking Operations[END_REF] Most of the unsaturated and saturated fatty acids were classified into groups 4 (101 ions) and 5 (82 ions). Oleic acid (OA, C 18 H 34 O 2 ) and palmitic acid (C 16 H 32 O 2 ; ER = 106.8) are the most abundant unsaturated fatty acid and saturated fatty acid, respectively. It is noteworthy that OA was detected in the particle phase and its peak level is ∼60% of the palmitic acid signal intensity. Furthermore, the signal intensities of unsaturated fatty acids decreased faster than those of saturated fatty acids after the gas stove had been turned off (Figure S4). This is reasonable because compared to saturated fatty acids, unsaturated fatty acids are more easily oxidized by OH radicals or O 3 due to the reactions with the CC bond. The addition of OH to the CC bond is found to be the dominant reaction pathway for unsaturated compounds. [START_REF] Atkinson | Atmospheric Degradation of Volatile Organic Compounds[END_REF] The uptake coefficients for the reactions of OH radicals with unsaturated compounds are also much larger than those of the OH reactions with saturated fatty acids. [START_REF] Atkinson | Atmospheric Degradation of Volatile Organic Compounds[END_REF] Moreover, both unsaturated and saturated fatty acids contribute to the formation of SOA in the indoor environment as they are established precursors of SOA in the atmosphere. [START_REF] Liu | Emission of Volatile Organic Compounds and Production of Secondary Organic Aerosol from Stir-Frying Spices[END_REF][START_REF] Ho | Dicarboxylic Acids, Ketocarboxylic Acids, α -Dicarbonyls, Fatty Acids, and Benzoic Acid in Urban Aerosols Collected during the 2006 Campaign of Air Quality Research in Beijing[END_REF] The compounds in group 6 (76 ions) represent a mixture of oxidation products and "background" compounds emerging from typical indoor materials (furniture, walls, etc.) that were present in the air before the experiments started. For example, the peak at m/z 60.0816 is assigned to trimethylamine (TMA, C 3 H 9 N; ER = 54.7) and might be the product of choline in the peanut oil through thermal decomposition. [START_REF] Ahn | Odor and VOC Emissions from Pan Frying of Mackerel at Three Stages: Raw, Well-Done, and Charred[END_REF] TMA can potentially form SOA during gas-phase reactions with OH radicals, O 3 , and NO 3 . [START_REF] Murphy | Secondary Aerosol Formation from Atmospheric Reactions of Aliphatic Amines[END_REF][START_REF] Silva | Trimethylamine as Precursor to Secondary Organic Aerosol Formation via Nitrate Radical Reaction in the Atmosphere[END_REF][START_REF] Tang | NO3 Radical, OH Radical and O3-Initiated Secondary Aerosol Formation from Aliphatic Amines[END_REF] The peak at m/z 133.1335 (C 6 H 16 ON 2 ; ER = 95.3) shows a trend very similar to that of the TMA profile (Figure S6), which might be a product from TMA. The NO x mixing ratio reached ∼300 ppbv during cooking, [START_REF] Liu | Assessing Indoor Gas Phase Oxidation Capacity through Real-Time Measurements of HONO and NO x in Guangzhou, China[END_REF] which is also comparable with the previous reports. [START_REF] Bartolomei | Combustion Processes as a Source of High Levels of Indoor Hydroxyl Radicals through the Photolysis of Nitrous Acid[END_REF][START_REF] Zhou | Time-Resolved Measurements of Nitric Oxide, Nitrogen Dioxide, and Nitrous Acid in an Occupied New York Home[END_REF] For the background compounds, the intensities decreased during cooking and then returned to the initial level after cooking.

3.2. Formation of Secondary Compounds. A total of 526 ions were detected as new products resulting from cooking, among which 229 ions were detected after cooking. By using KMD analysis, a series of homologous series were distinguished (Figure 2 and Figures S6-S8). Taking the ions with odd m/z values for example, we identified 14 homologous series of compounds containing a carbon atom (C), a hydrogen atom (H), and an oxygen atom (O) (Figure 2). Most of these compounds were found as fatty acids and their oxidation products. For example, compounds on line 1 are saturated fatty acids (C n H 2n O 2 ), and those on line 2 are monounsaturated fatty acids (C n H 2n-2 O 2 ).

Because the mechanism of heterogeneous OH oxidation of OA has been well studied, [START_REF] Atkinson | Atmospheric Degradation of Volatile Organic Compounds[END_REF][START_REF] Zhang | Time Resolved Study of Hydroxyl Radical Oxidation of Oleic Acid at the Air-Water Interface[END_REF][START_REF] Nah | OH-Initiated Oxidation of Sub-Micron Unsaturated Fatty Acid Particles[END_REF] 2 are likely produced by addition of OH to the CC bond of polyunsaturated fatty acids like LA. [START_REF] Nah | OH-Initiated Oxidation of Sub-Micron Unsaturated Fatty Acid Particles[END_REF] Temporal evolution profiles of these compounds have also been captured in real time (Figure 3); among all of the secondary products, 9-ON shows the highest signal intensity. The temporal profiles of OA, LA, and OH radicals are similar to each other, and meanwhile, the signal intensities of 9-ON and AA increase continuously after cooking, implying that 9-ON and AA are formed via the heterogeneous reaction of OH with OA or LA. However, the possibility that AA and 9-ON can also be formed by the reaction of OA or LA with O 3 in indoor air that exhibits higher ozone concentrations (in this study, the O 3 mixing ratio was <200 pptv) cannot be excluded. [START_REF] He | Kinetics Study of Heterogeneous Reactions of Ozone with Unsaturated Fatty Acid Single Droplets Using Micro-FTIR Spectroscopy[END_REF][START_REF] Hoffmann | On-Line Characterization of Gaseous and Particulate Organic Analytes Using Atmospheric Pressure Chemical Ionization Mass[END_REF] To illustrate the importance of OH radical reactions, we performed a simple calculation based on the reactivity of O 3 and OH radical on the particle surface [START_REF] Nah | OH-Initiated Oxidation of Sub-Micron Unsaturated Fatty Acid Particles[END_REF][START_REF] Zahardis | The Oleic Acid-Ozone Heterogeneous Reaction System: Products, Kinetics, Secondary Chemistry, and Atmospheric Implications of a Model System -A Review[END_REF] to compare the lifetimes of OA by reactions with OH radical and O 3 and removal of OA by air exchange with outdoor air. The calculations show that the most important removal process for OA is the reaction with the OH radical, because the effective uptake coefficient was reported to be 1.72 ± 0.08, indicating particle-phase secondary chain chemistry. [START_REF] Nah | OH-Initiated Oxidation of Sub-Micron Unsaturated Fatty Acid Particles[END_REF] Removal by reactions with O 3 would take ∼13 min (k first = 1.31 × 10 -3 s -1 ), [START_REF] He | Kinetics Study of Heterogeneous Reactions of Ozone with Unsaturated Fatty Acid Single Droplets Using Micro-FTIR Spectroscopy[END_REF] and removal by air exchange would take 1.6 h (taking the average exchange rate calculated in this work of 1.7 × 10 -4 s -1 ). The uptake coefficient of OH on LA is even higher than that on OA (3.75 ± 0.18), and the reaction of O 3 with LA (k first = 4.4 × 10 -4 s -1 ) is slower than that with OA, indicating that the removal of OA and LA is driven by OH radical chemistry. 9-ON and AA are most likely formed through the OH oxidation of OA, LA, or both simultaneously.

In conclusion, by using a novel real-time MS technique SESI-HRMS, a number of primary cooking emissions and secondary products formed upon reactions with OH radicals in indoor air have been detected and identified. High-temporal resolution data make it possible to monitor primary compounds and their intermediate and final products during and after cooking. In terms of OA, for example, 9-ON and AA, two well-known SOA materials, [START_REF] Ho | Dicarboxylic Acids, Ketocarboxylic Acids, α -Dicarbonyls, Fatty Acids, and Benzoic Acid in Urban Aerosols Collected during the 2006 Campaign of Air Quality Research in Beijing[END_REF][START_REF] Kawamura | Water Soluble Dicarboxylic Acids and Related Compounds in Antarctic Aerosols[END_REF][START_REF] Kawamura | Source and Reaction Pathways of Dicarboxylic Acids, Ketoacids and Dicarbonyls in Arctic Aerosols: One Year of Observations[END_REF] were identified as secondary products of OA due to OH reaction. We have demonstrated that OH radical reactions occur indoors as predicted in previous studies. [START_REF] Gomez Alvarez | Unexpectedly High Indoor Hydroxyl Radical Concentrations Associated with Nitrous Acid[END_REF][START_REF] Gligorovski | The Oxidative Capacity of Indoor Atmospheres[END_REF] We strongly recommend a future model study to better predict the formation of secondary products through reactions of O 3 and OH with primarily emitted organic compounds in indoor air; many secondary products formed through OH radical reactions can be even more toxic than the primarily emitted organic compounds. [START_REF] Gligorovski | Indoor (Photo)Chemistry in China and Resulting Health Effects[END_REF] This is especially important considering that a recent model study [START_REF] Won | Understanding the Spatial Heterogeneity of Indoor OH and HO2 Due to Photolysis of HONO Using Computational Fluid Dynamics (CFD) Simulation[END_REF] predicted that although OH radicals produced by HONO photolysis are located only in the sunlit regions of the indoor environment, the oxidation products are relatively well distributed in the entire room.

■ ASSOCIATED CONTENT

* sı Supporting Information

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.estlett.0c00044.

Detailed information about the cooking experiments and data analysis (Text S1); detection limit and sensitivity of SESI-HRMS, and sampling efficiency of the 6.7-m sampling tube (Table S1); detailed information about ions classified in six different groups (Tables S2-S7); detailed information about the enhancement ratio (ER) of identified compounds (Table S8); graphical illustration of the experimental room (Figure S1); temporal profiles of unsaturated aldehydes, unsaturated fatty acids, saturated fatty acids, and N-heterocyclic compounds from cooking (Figure S2); mass spectral fingerprints of the peanut oil heating experiment (Figure S3); timesignal intensity profiles of saturated and unsaturated fatty acids (Figure S4); cluster plot of ions detected during the control experiment (Figure S5); temporal profiles of ions detected at m/z 60.0816 and 133.1335 (Figure S6); and KMD plots of even and odd m/z values of ions detected in negative and positive ion detection mode (Figures S7-S9) (PDF)

Figure 1 .

 1 Figure 1. Cluster analysis of 864 ions with signal intensities of ≥1 × 10 5 au during and after cooking. Six groups distinguished by cluster analysis are numbered according to the order at which they reached their peak levels and represented by different colors: group 1, green; group 2, yellow; group 3, wine red; group 4, purple; group 5, blue; group 6, orange; background, black. The normalized signal intensity is presented by a color-coded scale; i.e., the signal intensity increases from dark blue (normalized value of -3) to wine red (normalized value of 3).

  all of the known intermediate and final products were carefully examined. On the basis of the accurate mass and isotopic pattern, LA, OA, and their intermediate and final products have been distinguished unambiguously, including C 18 H 34 O 2 [m/z 281.2483, OA, line 2 (Figure 2)], C 18 H 32 O 2 [m/z 279.2326, LA, line 4 (Figure 2)], C 18 H 35 O 5 (m/z 330.2408), C 18 H 35 O 4 (m/z 314.2459), C 18 H 34 O 3 [m/z 297.2429, line 5 (Figure 2)], C 9 H 16 O 4 [m/z 187.0967, azelaic acid, AA, line 8; ER = 45.5 (Figure 2)], and C 9 H 16 O 3 [m/z 171.1017, 9-oxononanoic acid, 9-ON, line 5; ER = 207.7 (Figure 2)]. Some compounds from line 10 (e.g., C 18 H 34 O 5 , m/z 329.2329) to line 14 (e.g., C 18 H 32 O 6 , m/z 343.2120) in Figure

Figure 2 .

 2 Figure 2. CH 2 KMD vs nominal KM of odd m/z values of ions detected in the negative ion detection mode during cooking (orange dots) and after cooking (blue dots). Ions found in the background (gray dots) represent compounds observed in the control experiments. The compounds that have the same general molecular formula but different numbers of CH 2 group appear on the same horizontal line: line 1, C n H 2n O 2 ; line 2, C n H 2n-2 O 2 ; line 3, C n H 2n O 3 ; line 4, C n H 2n-4 O 2 ; line 5, C n H 2n-2 O 3 ; line 6, C n H 2n O 4 ; line 7, C n H 2n-8 O 2 ; line 8, C n H 2n-2 O 4 ; line 9, C n H 2n-4 O 4 ; line 10, C n H 2n-2 O 5 ; line 11, C n H 2n-4 O 5 ; line 12, C n H 2n-2 O 6 ; line 13, C n H 2n-6 O 5 ; line 14, C n H 2n-4 O 6 .

Author Contributions

J.Z., Z.Y., and M.M. contributed equally to this work.

Notes

The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS