domain: phys chim info

text:
Florent Calvo
email: florent.calvo@univ-grenoble-alpes.fr

Aude Simon

Pascal Parneix

Cyril Falvo

Clément Dubosq

Infrared Spectroscopy of Chemically Diverse Carbon Clusters: A Data-Driven Approach

Keywords: Infrared spectrum, carbon clusters, statistical analysis, machine learning, density-functional tight-binding

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

The interest in the formation of fullerenes in the gas phase dates back from their experimental discovery. [START_REF] Kroto | C 60 : Buckminsterfullerene[END_REF][START_REF] Curl | On the formation of the fullerenes[END_REF] Recently, the astronomical detection of neutral [START_REF] Cami | Detection of C 60 and C 70 in na Young Planetary Nebula[END_REF] and cationic [START_REF] Berné | Interstellar C + 60[END_REF] fullerenes in the interstellar medium has revived this topic and notably led to several proposed formation mechanisms involving either smaller or larger building blocks. [START_REF] Dunk | Metallofullerene and fullerene formation from condensing carbon gas under conditions of stellar outflows and implication to stardust[END_REF][START_REF] Chuvilin | Direct transformation of graphene to fullerene[END_REF] Such astrophysical observations could be achieved by spectroscopy, for which the convenient signature of the highly symmetric buckminsterfullerene was the key to its identification. Spectroscopy is undoubtedly the method of choice on which astrochemists rely to better understand the fullerene road under such harsh astrophysical conditions. Until now, the traditional strategy for interpreting astrophysical spectra of complex molecules, in both the infrared (IR) or optical ranges, has mostly rested on the comparison with spectra usually obtained from electronic structure calculations offering a sufficient trade off between accuracy and efficiency, density-functional theory (DFT) and its time-dependent version representing the most popular methods of choice. [START_REF] Bauschlicher | Infrared spectra of pilycyclic aromatic hydrocarbons (PAHs)[END_REF][START_REF] Malloci | Electronic absorption spectra of PAHs up to vacuum UV[END_REF] In some cases such as polycyclic aromatic hydrocarbons (PAHs), further refinements have been attempted in order to account for anharmonicity or finite temperature effects in vibrational spectra, [START_REF] Calvo | Temperature effects on the rovibrational spectra of pyrene-based PAHs[END_REF][START_REF] Simon | Molecular dynamics simulations on [FePAH] + π-complexes of astrophysical interest: anharmonic infrared spectroscopy[END_REF][START_REF] Mackie | The anharmonic quartic force field infrared spectra of three polycyclic aromatic hydrocarbons: Naphthalene, anthracene, and tetracene[END_REF] or to model the IR emission cascade spectrum originating from a UV excitation. [START_REF] Basire | Modeling the anharmonicinfrared emission spectra of pahs: Application to the pyrene cation; in PAHs and the Universe[END_REF] Unraveling the molecular compounds responsible for the spectroscopic features observed in astrophysical objects like planetary nebulae can also be tackled from a more statistical perspective, without assuming specific molecular structures but focusing instead on broad relevant samples generated by intuition or, in a much less biased way, by means of high-throughput atomistic simulations. [START_REF] Bonnin | Simulating the structural diversity of carbon clusters across the planar-to-fullerene transition[END_REF] In particular, upon categorizing conformers into different structural families, their generic spectral features can be inferred from the individual contributions of each member. Such an approach was found fruitful in assisting the interpretation of astronomical spectra of planetary and protoplanetary nebulae, [START_REF] Dubosq | Mapping the structural diversity of C 60 carbon clusters and their infrared spectra[END_REF] highlighting in particular the potential role of compounds with a significant aromatic content such as fullerenes. While such a statistical averaging washes away the detailed spectroscopic information pertaining to individual conformers, it is consistent with the relatively poor knowledge of such astrophysical regions and the likely presence of many chemical species. However, extending such a bottom-up strategy to broader classes of compounds having a larger size or containing heteroatoms is not straightforward, because of the very fast increase in the number of stable conformers and the difficulty of sampling their potential energy landscape as exhaustively as possible. Moreover, for each conformer a suitable method is needed to determine their spectrum, which for chemical accuracy requires either advanced (polarizable) force fields or schemes offering an explicit description of electronic structure, even simplified.

While the two tasks of sampling the energy landscape and determining the individual spectra can be addressed successively along the lines of a multiscale description, the computational effort associated with the spectral determination can be particularly heavy for large samples, even with efficient methods such as density-functional based tight-binding (DFTB). However, extracting the spectral features from a statistical sample seems also naturally suited to be tackled by machine learning (ML) techniques, a broad range of computational approaches that have become spectacularly popular in chemical physics and physical chemistry in the recent years. [START_REF] Ceriotti | Machine learning meets chemical physics[END_REF][START_REF] Prezhdo | Advancing Physical Chemistry with Machine Learning[END_REF] Within the general context of relating properties to structure, several groups have recently shown the benefits of employing ML for vibrational spectroscopy [START_REF] Gastegger | Machine learning molecular dynamics for the simulation of infrared spectra[END_REF][START_REF] Hu | Machine learning protocol for surface-enhanced Raman spectroscopy[END_REF][START_REF] Härkönen | Bayesian quantification for coherent Anti-Stokes Raman scattering spectroscopy[END_REF] through a variety of approaches aiming to represent potential energy and electric dipole moment surfaces within perturbative frameworks, [START_REF] Lam | Combining Quantum mechanics and machinelearning calculations for anharmonic corrections to vibrational frequencies[END_REF] to condense molecular information into topological descriptors, [START_REF] Kovács | Machine-learning prediction of infrared spectra of Interstellar Polycyclic Aromatic Hydrocarbons[END_REF] or to numerically solve the quantum nuclear dynamics problem by better partitioning the various degrees of freedom. [START_REF] Gandolfi | Machine learning for vibrational spectroscopy via divide-and-conquer semiclassical intial value representation molecular dynamics with application to N -methylacetamide[END_REF] In the present contribution, we explore several ML ideas to reconstruct the infrared spectrum of carbon clusters in a statistical sense, from a limited sample and using interpolation techniques in a multidimensional feature space, supervision being introduced through metric learning. We notably show how such an interpolation scheme can predict the spectral trends for clusters that are not members of the initially chosen sample.

We also use a clustering algorithm to further address the classification of conformers into structural families.

In the next section, we briefly describe how the samples of conformers were obtained and the infrared spectra reconstructed from these samples. In Sec. 3 we present several statistical analyses of the spectra in terms of their convergence within the sample, their reconstruction by interpolation assisted by metric learning, and we discuss the typical error of the method on individual and collective spectra. Sec. 5 discusses further the classification of the structures into families based on iterative label spreading [START_REF] Parker | Selecting appropriate clustering methods for materials sciences applications of machine learning[END_REF] rather than ad hoc order parameters alone. Finally, some concluding remarks are given in Sec. 6.

Sample generation

Our systems of interest are carbon clusters C n of selected sizes n, which exhibit a great diversity of structures that have been the specific subject of earlier contributions. [START_REF] Bonnin | Simulating the structural diversity of carbon clusters across the planar-to-fullerene transition[END_REF][START_REF] Kim | Melting the fullerenes: A molecular dynamics study[END_REF] Here we summarize the main steps of the sample generation, and how the infrared spectra were obtained from it.

Force field exploration

The samples of isomers for C 24 , C 42 , and C 60 were originally produced by replica-exchange molecular dynamics simulations based on the REBO potential, [START_REF] Brenner | A second-generation reactive empirical bond order (REBO) potential energy expression for hydrocarbons[END_REF] with details given in Ref. 13. For these simulations, the clusters were enclosed in spherical containers to prevent dissociation at high temperature, and various simulations using different radii were performed to enhance sampling, fixing the densities at ρ = 0.025, 0.15, 0.4, and 1.7 g.cm -3 for the three cluster sizes. These densities were chosen on a trial-and-error basis in order to generate clusters with different structural trends. In particular, lower densities are needed to favor the very extended branched conformers, while higher densities tend to produce cages more efficiently. New simulations were also performed for C 33 and C 52 to assess the performance of the interpolation scheme introduced below, away from the testing set that only includes data for C 24 , C 42 and C 60 . For these additional clusters, we used a more efficient criterion based on bond connectivity in a parallel tempering Monte Carlo framework: a configuration was rejected if two subclusters are separated from each other by more than 3 Å. Except for this difference, the parameters were otherwise similar, with 28 temperatures distributed in a geometrical fashion in the range 500-5500 K, and 4 additional temperatures around the melting point.

The total numbers of distinct conformers obtained with the REBO potential are listed in Table 1 for the five clusters.

Reference infrared spectra

Because it lacks any information about the dipole moment surface, the REBO potential energy surface used to generate the structures is not suited to infrared spectroscopy.

Instead we resorted to the self-consistent-charge density-functional-based tight-binding method [START_REF] Elstner | Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties[END_REF] for this purpose, uusing the dispersion correction parameters of Ref. 14. The individual absorption spectrum associated with each structure was determined in the harmonic approximation after local reoptimization. Throughout this article (and in the supplementary material), this method will be denoted simply as DFTB. Structures saved periodically and locally reoptimized using the REBO potential were further refined using the DFTB method to calculate their IR absorption spectrum. The numbers of distinct structures obtained after this reoptimization step, listed in Table 1, are always lower than the initial sample size, because occasionally several starting points endded up into the same minimum (a reduction in information and pool size would also naturally occur by conducting the optimizations in the opposite way). Here we should also note that the loss of a significant fraction in the initial sample size could be partly alleviated using ML techniques that screen the structure before deciding whether they are worth considering for further optimization. [START_REF] Jesus | Modeling microsolvation clusters with electronic-structure calculations guided by analytical potentials and predictive machine learning techniques[END_REF] A standard double harmonic approximation was employed to calculate the IR absorption spectrum with the DFTB method, without any scaling factor for the frequencies. This choice is motivated by the methodological nature of the present study, rather than on any intention of a direct comparison with existing experiments or observations. It should be noted that, for the present systems, different scaling factors should also be applied in different spectral ranges for such a comparison, [START_REF] Dubosq | Mapping the structural diversity of C 60 carbon clusters and their infrared spectra[END_REF] which also leads to narrow spectral ranges being uncovered and appearing as spurious holes.

The spectra of entire structural families were determined as in Ref. 14 by simple addition over all contributions from individual members and no particular weighting. Individual spectra were made smoother by Gaussian broadening with a 5 cm -1 width over the 0-3000 cm -1 range with a 3 cm -1 bin size. Such parameters were necessary in order to satisfactorily describe both the broad plateaus in the 500-1500 cm -1 range, as well as the few narrow peaks occuring at various frequencies depending on the structural family.

Finally, and following our earlier work, [START_REF] Dubosq | Mapping the structural diversity of C 60 carbon clusters and their infrared spectra[END_REF] the high-frequency part of all IR absorption spectra was further attenuated using the blackbody radiation law at 300 K in order to make comparison with IR emission astrononomical spectra more realistic. All DFTB calculations were performed with the deMonNano software package. [START_REF] Heine | deMonNano[END_REF]

Structural classification

At the sizes n covered here of n = 24, 42, and 60, carbon clusters exhibit a broad variety of structures that can be categorized using order parameters. Following our earlier analysis [START_REF] Bonnin | Simulating the structural diversity of carbon clusters across the planar-to-fullerene transition[END_REF][START_REF] Dubosq | Mapping the structural diversity of C 60 carbon clusters and their infrared spectra[END_REF] we use the fraction of sp 2 atoms and a dimensionless asphericity shape parameter β to sort these structures into four main families: (i) the cages family, which includes the archetypal fullerenes, consists of a high sp 2 content and a low asphericity; (ii) the flakes family, notably including graphene cuts, also has a high fraction of sp 2 atoms but a significant asphericity; (iii) pretzel-like structures, [START_REF] Kim | Melting the fullerenes: A molecular dynamics study[END_REF] containing several chains of sp 1 atoms but not necessarily associated with marked deformations from the sphere; (iv) branched structures also dominated by long chains of sp 1 atoms but significantly aspherical.

The shape parameter β is one of the three quantities that can be obtained by assigning the three principal momenta of inertia I k of the conformer following the Hill-Wheeler representation as

I k = 2 3 r 2 c [ 1 + β sin ( γ + (4k -3)π 6 )] , k = 1, 2, 3, (1)
where r c is the gyration radius and the angle γ measures the cluster triaxiality. With such a definition, β lies in the range 0-1, small values indicating nearly spherical shapes while high values are typical of prolate ellipsoids. The fraction of sp 2 atoms is straightforwardly obtained in the DFTB scheme from the orbital populations. [START_REF] Dubosq | Mapping the structural diversity of C 60 carbon clusters and their infrared spectra[END_REF] In the REBO model, it is simply estimated from the relative number of carbon atoms with exactly three nearest neighbors.

With such dimensionless order parameters at hand, all conformers were categorized into either of the four families, with sizes that are given in Table 2 for each of the five clusters. This partitioning of isomers between structural classes was initially introduced without reference to their spectroscopic response, and in the following we challenge this clas-sification using additional tools.

3 Statistical analyses

Convergence of global spectra

For each cluster size, the IR spectra associated with each of the four structural families result from the combination of all individual spectra from each member of the family, and we first discuss how the resulting spectra are sensitive to the diversity among the corresponding sample. Individual IR spectra are generally highly resolved, because the clusters are not very large and display at most 172 fundamental modes, of which only a subset are IR active. However, the combination of many such spectra into an average spectrum for the corresponding family will likely yield much smoother features all the more than the family is large and diverse.

Here we introduce an error measure to quantify the discrepancy between a reference spectrum I ref (ω) and an estimated spectrum I(ω), from the average sum of the integrated difference between them on the relevant interval of interest:

E = 1 ω max -ω min ∫ ωmax ω min |G ref (ω) -G(ω)| dω, ( 2
)
G ref (ω) = ∫ ω ω min I ref (x)dx G(ω) = ∫ ω ω min I(x)dx
with ω min = 0, ω max = 2500 cm -1 . This error measure is not as sensitive to spectral shifts as the more standard least square error in which the first integral operates directly on |I ref -I| or its square, and this was notably recognized recently by Kovacs and coworkers in the ML approach to the IR spectroscopy of polyaromatic hydrocarbons. [START_REF] Kovács | Machine-learning prediction of infrared spectra of Interstellar Polycyclic Aromatic Hydrocarbons[END_REF] Here the reference spectrum is simply defined for each family as the spectrum obtained from the full corresponding sample.

To analyse the sensitivity of the spectra towards the structural diversity among the set, Monte Carlo simulations have been performed in which the spectra were reconstructed from a random limited subset. For each structural family, 1000 random subsets of conformers were chosen and the spectra averaged from their specific contributions.

The average errors obtained for the four families of C 60 isomers are shown in Fig. 1 as a function of sample size. The converged spectra obtained from the full samples are themselves also shown as insets for subsequent reference, together with representative conformers of each family for visual identification.

Clearly the IR spectra obtained for the cages and flakes do not depend as much on the underlying sample as the more disordered pretzels and branched conformers, for which the average error is about one order of magnitude larger. The spectral features of cages and flakes are reasonably converged already with a few thousands of the corresponding structures, while this number is closer to tens of thousands for pretzels, and exceeds it for branched structures.

Similar results are obtained for C 42 , but for C 24 the numbers of cages and flakes are much smaller and the resulting spectra show a much greater dependence on the subset selection. The corresponding spectra and their statistical convergence are given as Supplementary information.

An interpolation scheme for IR spectra

The previous statistical analysis revealed that the spectral features associated with structural families are differently robust depending on the diversity within each family. We now attempt to use this information to reduce the amount of conformers needed to reproduce the global spectra in a more systematic and automated fashion that could also make it transferable to other systems.

More precisely, we use interpolation to relate the spectrum I R (ω) of an arbitrary conformer R to the known spectra I k (ω) from reference conformers R k of a fixed sample, much smaller than the complete set of conformers generated by simulation. In a first natural approach, linear interpolation is employed as

I R (ω) = ∑ k g R,k I k (ω). (3)
where the weight g R,k measures the degree of similarity between structures R and R k , and is a simple but nonlinear function of a distance d R,k . Here we considered two kernels, namely Gaussian regression (GR) g(d) = g 0 exp(-γd 2 ) with g 0 a normalization factor, as well as inverse distance weighting (IDW)

g(d) = 1/d 2 if all d > 0, or, if d = 0
for some member of the set, then g(d) = 1 for this member and g(d) = 0 for all other members.

The distance d R,k is defined by a metric based on N descriptors that we collectively denote as q(R) = {q j , j = 1, . . . , N }. These descriptors were chosen to cover various molecular features of relevance in vibrational spectroscopy, and they are detailed below.

Next we select a reduced sample of conformers on which to interpolate the spectra for arbitrary structures. In a first approach, we mesh the order parameter space in β and sp members only, or 13314 distinct structures.

Once descriptors are known, a metric can be chosen for d and the interpolation method can be applied straightforwardly to predict the IR spectra from the reference spectra on the reduced sample. However, the quality of the prediction can be improved using ideas from machine learning methods, and here metric learning will be employed.

Descriptors encoding molecular structures

Twenty features or descriptors were chosen to characterize the various isomers generated by molecular simulation. In the context of vibrational spectroscopy we have primarily used descriptors involved in the local bond properties involving between 2 and 4 con-nected atoms, as well as more global properties related to the overall shape of the entire configuration or its bond topology seen as a graph. The descriptors are detailed and numbered as follows. As they should, they are all invariant over global translations, rotations, and permutations among atoms.

For a N -atom cluster, and from the gyration tensor S = {S αβ },

S αβ = 1 N ∑ i r α i r β i ,
where r α i is the component α = x, y, z of atom i, we can extract its trace, giving the square gyration radius R 2 g , while the higher moments give insight into the so-called asphericity and prolateness parameters A and P . [START_REF] Wei | On shape asymmetry of Gaussian molecules[END_REF] Introducing the traceless tensor D we thus define the three first descriptors as

R 2 g = Tr S (D1) A = 3 Tr (D 2 ) 2R 4
g (D2)

P = 3 Tr (D 3 ) 2R 6 g (D3) D = S - R 2 g 3 I,
where I is the 3×3 identity matrix. Note that only R 2 g is size sensitive, A and P being both dimensionless quantities. Also note that the asphericity parameter A is not strictly identical to the parameter β used to classify the structures into the four families, for which Eq. ( 1) was used.

From the structure we then define bond connectivities using a cut-off radius of 1.85 Å, the average bond length and the fluctuations around this average providing our next descriptors

r = ⟨r ij ⟩ nn , (D4) δr = [ ⟨r 2 ij ⟩ nn -⟨r ij ⟩ 2 nn ] 1/2 , ( D5
)
where the subscript nn indicates that average is taken over nearest neighbors. Likewise, for any triplet i, j, k of connected atoms we define the average angle θ and the fluctuation around this average as

θ = ⟨θ ijk ⟩ nn (D6) δθ = [ ⟨θ 2 ijk ⟩ nn -⟨θ ijk ⟩ 2 nn ] 1/2 , ( D7
)
and similarly torsion angles ϕ ijkℓ are identified and included in the set for all quadruplets i, j, k, ℓ of connected atoms:

φ = ⟨ϕ ijkℓ ⟩ nn (D8) δϕ = [ ⟨ϕ 2 ijkℓ ⟩ nn -⟨ϕ ijkℓ ⟩ 2 nn ] 1/2 . ( D9
)
From the global connectivity we also calculate the adjacency matrix A ij , which we diagonalize into the set of eigenvalues {α k } from which the next descriptors are extracted as

α min = min k α k (D10) α max = max k α k (D11) ᾱ+ = ⟨α + k ⟩ (D12)
In the above equations, ⟨α + k ⟩ denotes the average over all strictly positive eigenvalues (the trace and average over all eigenvalues being strictly zero for any adjacency matrix).

The four next descriptors are also related to bond connectivity but in more direct connection with the chemical nature of the carbon bonds. More specifically, the numbers of atoms that are singly (N 1 ), doubly (N 2 ), triply (N 3 ), and quadruply (N 4 ) connected to other carbon atoms are identified and their fractions define the descriptors as

f 1 = N 1 /N (D13) f 2 = N 2 /N (D14) f 3 = N 3 /N (D15) f 4 = N 4 /N (D16)
In a first approximation, these fractions measure the relative amounts of carbon atoms that terminate a chain or are hybridized as sp 1 , sp 2 , and sp 3 , respectively. In particular, f 3 is used together with β to assign the REBO structures into the four families.

The next descriptor is also topological and defined as the normalized meshedness χ of the graph made by carbon atoms. [START_REF] Buhl | Efficiency and robustness in ant networks of galleries[END_REF] It is defined from the numbers of edges N and bonds

N b as χ = N b -N + 1 2N -5 (D17)
Finally, in relation with the aromatic nature of many carbon nanostructures, we have considered the numbers N k of cycles of length k = 5-7 and normalized them to yield three more descriptors

N5 = N 5 /N (D18) N6 = N 6 /N (D19) N7 = N 7 /N. (D20)

Principal components analysis

Each structure in the samples generated for C 24 , C 42 and C 60 was assigned a point in 20dimensional space using the above listed descriptors. Because the various dimensions possibly cover quite different numerical ranges, we further standardize the dataset by shifting each value with respect to the average, and dividing it by the mean square fluctuation. To assess the quality of the dataset and remove possible redundancies among descriptors, we next analyze the data in terms of their principal components (PCs).

After diagonalization of the PC matrix, the eigenvalues are ordered by decreasing value, and the eigenvectors corresponding to the highest eigenvalues are employed to determine the number of dimensions needed to capture the greatest amount of the data and their dispersion. Such an analysis is performed independently for the three datasets corresponding to C 24 , C 42 and C 60 samples. The percentage of explained variance obtained from the eigenvectors of the PC matrix is given as an inset in Fig. 3 for the lowest 10 dimensions that correspond to the highest eigenvalues. Most of the variance can thus be explained with only 5 principal components, especially for the samples corresponding to the two larger clusters, the data being somewhat more widespread for C 24 . We interpret this difference as reflecting the poorer ability of this small cluster to form regular cages, which contribute significantly to the other samples but are only very few (eleven) in this case.

The quality of the representation of the data by the 5 principal components associ-ated with the highest eigenvalues of the PC matrix was determined from the sum of the corresponding square eigenvector coefficients. It is shown in the main part of Fig. 3 for the three samples against descriptor number, and reveals the different contributions of these descriptors to the set of 5 principal components. While the results for C 24 are again slightly different from the two other systems, we find some common patterns that were not anticipated when the descriptors were originally designed: the angular terms for triplets of connected atoms essentially do not contribute, as is the case for the normalized numbers of 6-and 7-atom rings. In contrast, we find the shape parameters of asphericity and prolateness, and the fractions of atoms with 2 and 3 neighbors to be quite well represented in those principal components. These results are consistent with our initial but empirical choice of order parameters to sort the structures into the four families of cages, flakes, pretzels, and branched isomers. [START_REF] Bonnin | Simulating the structural diversity of carbon clusters across the planar-to-fullerene transition[END_REF]

Metric learning

The previous PC analysis indicates that 10 descriptors are sufficient to capture the structural diversity amont the available set, and we thus discard those contributing the least to the 5 highest components, namely the average bond length and its fluctuation, the average bond angle and its fluctuation, the fluctuation in torsion angles, the maximum and average of positive eigenvalues of the adjacency matrix, and the fractions of 5-, 6-, and 7rings.

Distances are thus evaluated in the remaining 10-dimensional set, and for each kernel g(d) we introduce a specific metric d(R) as a set of strictly positive numbers a j such that

d 2 R,k = N ∑ j=1 a j [q j -q (k) j ] 2 ,
where q j and q (k) j denote the corresponding descriptors for conformers R and R k , respectively.

To improve the quality of the interpolation scheme, we then introduce some degree of supervision to the algorithm by training the metric so the spectra predicted for a specific set of structures (a testing set) mimics as best as possible their true IR absorption spectra determined by the DFTB method. The testing set was chosen by randomly selecting 300 conformers among the totally available sample, 100 for each of the C 24 , C 42 and C 60 systems, only ensuring none of them also belongs to the interpolation sample as this would make error evaluation non differentiable with the IDW kernel. The testing sample is also depicted on the density plots in Fig. 2.

For a given kernel, the metric is then optimized by minimizing a global error E over the whole testing set using again a Monte Carlo procedure, following here a zero temperature Metropolis acceptance rule in which random moves that do not decrease the error are rejected, alternating with steps of gradient-based local minimizations.

From a practical perspective, all statistical computations were performed with homemade codes.

4 Performance of the interpolation scheme

Individual IR spectra

The general behavior of the interpolation scheme can be first discussed on the example of individual structures. Figure 4 shows the reference and interpolated IR spectra obtained for members of the testing set that contribute the most and the least to the final error obtained after metric learning, using the GR and IDW kernels. The corresponding conformers are also shown on this figure. As already mentioned, the IR spectra of specific conformers are highly resolved for such relatively small clusters. By combining the contribution of various individual spectra, the interpolated spectrum is expected a priori to be much smoother and convey the statistically dominant features. For the two kernels, the best prediction is obtained for cage isomers of C 60 , whose spectra exhibit rather large absorption bands near 700 cm -1 and a main peak, as well as in the 1000-1700 cm -1 range with a broader bump. Here the interpolated spectrum manages to reproduce these features and especially the low intensity parts.

In contrast, the most dissimilar spectra are obtained for much more disordered, branchedtype structures exhibitting fewer but highly intense IR active modes. The same conformer is responsible for this highest contribution to the error for both kernels. Interestingly, the Gaussian weight manages to preserve this highly resolved character, but the most intense peak near 1150 cm -1 is not reproduced. This suggests that none of the reference spectra in the interpolation sample matches this specific spectrum well enough to be able to produce anything similar in output. In turn, it suggests excessive structural dissimilarities between this specific conformer and all members from the sample. Inverse distance weighting produces an excessively smooth spectrum lacking intense features, as the result of multiple conformers from the sample producing similar weights.

These results indicate that the interpolation scheme can produce very different types of spectra, typically smoother but not necessarily always smooth either since the weighting scheme is highly nonlinear. For comparison, we have also shown in Fig. 4 its prediction without carrying the metric learning stage, i.e. including all 20 descriptors in the metric and not optimizing its coefficients a j . The most similar spectra are weakly altered, except at high frequency where spurious peaks are found near 2000 cm -1 . With inverse distance weighting, the interpolated spectrum is also even smoother and the error to the reference spectrum is increased by about 25%. Concerning the most dissimilar spectra, the absence of supervision produces an even worst prediction for Gaussian regression, with spurious peaks near 550 cm -1 and 2000 cm -1 , without improving on the most intense peak at 1150 cm -1 . However, the much broader averaging achieved by inverse distance weighting is preserved without metric learning, the predictions being very similar.

Global IR spectra

The global spectra predicted for the four structural families of C 60 after averaging on all their members are shown in Fig. 5, in comparison with the reference DFTB spectra and for the two interpolation schemes employing the GR or IDW kernels. The corresponding spectra predicted for the families of C 24 and C 42 are given as supplementary information.

Overall, the spectral features are fairly well reproduced by both interpolation schemes, in the entire relevant spectral range. In agreement with the purely statistical analysis of Fig. 1, the greater discrepancies are found for the pretzels and branched families, for which convergence of the global spectrum is the slowest with set size. The most significant deviation occurs for the flakes, in the 1000-1200 cm -1 range, where interpolation with inverse distance weighting notably underestimates the IR intensity by a few percents, while the narrow peak at 2150 cm -1 is overestimated.

Without metric learning, the spectrum obtained for cages is barely affected, but more significant differences are found for the other families, especially with inverse distance weighting where the broad region 1200-2200 cm -1 is notably underestimated. It is also quite striking that the spectral features of the pretzels and branched families are very similar, except for minor intensities differences in the broad band at 600 cm -1 .

This spectroscopic similarity obviously results from structural similarities that can be ascribed to the predominance of sp 1 carbon chains in both cases, as discussed further below.

The performance of the interpolation scheme was also tested on a more challenging case, namely a carbon cluster different from any member of the sample and testing sets. The energy landscape of C 52 explored using a combined REBO-DFTB methodology provided broad samples of cages, flakes, and branched structures, although pretzels are fewer, their numbers being given in Table 2. Their IR spectra were determined for each conformer individually using DFTB, and accumulated to yield a global spectrum for each of the four families, serving again as reference.

The interpolation schemes were also applied to predict the IR spectra from the sole information of the conformers provided by their descriptors. These spectra are depicted in Fig. 6. For this cluster the main features of the IR spectra are reproduced very satisfactorily for the four families, the main deviations being found again in the 1000-1200 cm -1 Figure 6: IR spectra predicted for structural families of C 52 by sample interpolation using the GR and IDW kernels (red and blue lines, respectively), compared to the DFTB reference spectra (thick brown lines).

range for the flake conformers, Gaussian regression performing somewhat better than inverse distance weighting. The spectrum for the pretzels is more resolved owing to the much smaller size of this family, but the trends are again correctly described by both interpolation schemes. For the cages, the reference spectrum is mostly in error near 1000 cm -1 , while the spectrum for the branched family is remarkably described.

The interpolation schemes were also applied to the structures directly produced at the initial REBO level of modeling, without refining them using DFTB. The results, given as supplementary information, remain satisfactory for the more ordered cages and flakes families but for the disordered structures several spurious peaks are produced upon interpolation, especially when employing Gaussian regression, indicating that some conformers are not correctly recognized from their descriptors. This issue could probably be addressed by training the metric on a set made from REBO structures, rather than the DFTB conformers presently used. For the smaller cluster C 33 , the number of cages and pretzels is rather reduced (see Table 2) and the resulting spectra are not particularly

smooth. Yet the trends predicted by the interpolation schemes are again rather good, the main discrepancy appearing once again for the main band of the flakes family spectrum near 1000-1200 cm -1 , which further indicates that not enough representative members of these specific conformers are included in the sample. As was the case for C 52 , the interpolation scheme for this smaller cluster performs not as well when the descriptors are those of the REBO structures (see supplementary information), although it is worth noting that the smoothest spectrum, obtained for the flakes family, is not affected to such an extent.

Clustering analysis by iterative label spreading

The above analysis has shown that interpolation on a much reduced sample of representative data can reproduce the structural trends of collective but also individual spectra, especially with metric learning. It was also found that the IR spectra obtained for two of the four families are very similar, namely the pretzels and the branched structures.

They were arbitrarily distinguished from one another based on the asphericity parameter, pretzel conformers being more spherical, hence having a lesser tendency for exhibitting terminating carbon chains. Here were return to this issue of structural partitioning but ignoring now such intuitive criteria and turning instead to the purely statistical analysis of iterative label spreading (ILS). [START_REF] Parker | Selecting appropriate clustering methods for materials sciences applications of machine learning[END_REF] ILS is a semisupervised clustering method chosen for its ability to handle dense sets of points. Briefly, ILS proceeds by reordering all members of the set, labelling them each after the other, picking the next labelled member as the closest from the already labelled set. Once all members are labelled, groups of points within the set forming families are expected to be separated by peaks in the minimum distance. ILS requires a metric for distance evaluation between two arbitrary members of the set, and we naturally use the metrics optimized for the interpolation scheme for either Gaussian regression or inverse distance weighting. As in Ref. 23 we chose the starting points as the set barycenters. Fig. 7 shows the resulting partitioning for the structural database of C 60 , as emerging from ILS with the metric optimized for Gaussian regression. Similar results are obtained with the other metric optimized for inverse distance weighting. To relate these partitionings with our initial sorting into the four families of cages, flakes, pretzels and branched structures, we further tag each labelled member with the value of its sp 2 fraction, which is the most discriminating of the two order parameters with the asphericity β. The corresponding graph tagged with the value of β is given as supplementary information. While there is no obvious jump in the minimum distance, the structures are clearly ordered within the four families identified empirically. The first ones lying closer to the set barycenter are the cages (highest sp 2 , lowest β content) and are followed by the flakes (medium sp 2 , higher β). The more disordered structures follow, first with the branched structures (lowest sp 2 values, highest β) and finally the pretzels (low to medium sp 2 , low β).

These results, and especially the proximity between pretzels and branched structures, are fully consistent with their spectroscopic similarities already noted in Figs. 5 and6 which themselves originate from the common chemical motifs of long sp 1 carbon chains.

However, the ILS analysis also suggests the possibility of subgroups inside some of the four main families, particularly for the more disordered pretzels and especially branched cases. Further partitioning among these families would require additional order parameters for a more complete characterization.

Unravelling the relation between structure and vibrational or electronic spectra is one of the important issues in modern chemical physics in the gas phase. The interpolation scheme explored in this work attempts to address this issue from an automated statistical perspective. In this respect it is part of the broad current effort to use machine learning techniques to assist spectral determination. The interpolation scheme developed here is relevant for extracting spectroscopic trends over a statistical set of structures, and cannot claim chemical accuracy for individual conformers. Its main appeal resides in its systematic nature and in not being limited to chemically similar structures such as pure aromatics. While it can be used as a nonsupervised method, optimization of the metric was found to improve the performance especially in the cases where the structural diversity was the highest among the set. As with conventional machine learning approaches, the performance of the present method is naturally limited by the existence of sufficiently similar members in the sample on which interpolation is performed.

With respect to neural networks, the present approach involves far fewer parameters but requires a conformational sample. Its main computational interest was motivated in the need to bypass the heaviest numerical effort associated with the determination of numerous individual IR spectra, structural optimization itself being mostly achieved at a lower level of theory (here, the REBO force field).

Classification of the structures into the four families of cages, flakes, pretzels and branched conformers was also confirmed by performing independently a clustering analysis based on iterative label spreading. Here it could be interesting to compare the predictions of this method with deep learning approaches that are also becoming more widespread in physical chemistry as a classification tool.

In the future the interpolation scheme could be extended and improved along several directions. For pure computational efficiency, the k-nearest neighbor algorithm could be used to limit the number of members in the sample to be included for each newly pre-dicted spectrum, even though the number of neighbors and the metric should be both optimized themselves self-consistently as hyperparameters. Another natural improvement could be in the selection of the interpolating sample, which was here taken from a systematic mesh of conformers based on dedicated order parameters. A Monte Carlo procedure could be introduced to further optimize the members of this sample in their capability to describe the spectra of the testing set. Alternatively, Bayesian inference could be envisaged as well for even smaller samples.

From the physical chemistry perspective, the interpolation method could also be applied to systems exhibiting even greater chemical diversity, such as hydrogen-, oxygenor nitrogen-containing compounds, for which REBO or other reactive potentials [START_REF] Ni | A reactive empirical bond order (REBO) potential for hydrocarbon-oxygen interactions[END_REF][START_REF] Senftle | The ReaxFF reactive forcefield: development, applications, and future directions[END_REF] could be employed. The inverse problem of finding pools of structures that match observational data could also be tackled stochastically, under given constraints of available samples in terms of sizes, compositions, and individual conformations. Here also, ML techniques appear as naturally suited to such a purpose.

Figure 1 :

 1 Figure 1: Average error in describing the IR spectra of the four structural families of C 60 clusters, by selecting only a subset of the entire corresponding sample, with error bars that indicate the standard fluctuations around the average. The insets show the reference IR spectra obtained for the full samples for each family. Next to each panel corresponding to a specific family, a typical structure is also depicted.

Figure 2 :

 2 Figure 2: Distributions of locally optimized structures of carbon clusters C n projected onto the plane of sp 2 fraction and asphericity shape parameter β, in logarithmic scale. The full blue and empty green circles mark the location of the interpolation sample and the testing sets members, respectively.

Figure 3 :

 3 Figure 3: Quality of the description provided by using only the principal components corresponding to the 5 highest eigenvalues of the PC matrix, for the three samples of C 24 , C 42 , and C 60 . The inset shows the percentage of variance explained by the successive principal components ordered by decreasing eigenvalue, up to dimension 10.

Figure 4 :

 4 Figure 4: Most similar ('best') and dissimilar ('worst') individual IR absorption spectra, as defined from their lowest and highest contribution to the global error among the testing set, obtained for the GR (left panels) and IDW (right panels) kernels after metric learning (blue lines), or without (red lines). The DFTB reference spectra are shown as thick brown lines, and the corresponding conformers are depicted next to the spectra.

Figure 5 :

 5 Figure5: IR spectra of the four structural families of C 60 clusters, obtained from DFTB electronic structure calculations (thick brown lines) and from sample interpolation employing the GR and IDW kernels (red and blue lines, respectively). The dashed lines show the corresponding predictions of the interpolation method without metric learning.

Figure 7 :

 7 Figure7: Clustering analysis for the isomers of C 60 by iterative label spreading, from the metric optimized for the GR kernel. The various minima are numbered (abscissa) by increasing minimum distance (ordinate) to the set of already labelled members. Each line is colored according to its fraction of sp 2 atoms.

33

 33

Table 1 :

 1 Numbers of distinct isomers in the samples produced by the simulations with the REBO potential, and after local reoptimization at the DFTB level.

	Cluster REBO	DFTB
	size	sample sample
	24	51 901 44 341
	33	62 974 52 682
	42	240 305 196 519
	52	260 746 150 475
	60	656 438 309 167

Table 2 :

 2 Sizes of the configurational samples for all carbon clusters considered in this work, as described by the DFTB method.

	System cages flakes pretzels branched
	C 24	11	740	6312	37282
	C 42	27350 40740	36013	92416
	C 60	82111 78182	38141	110733
	C 33	13	4764	202	47703
	C 52	20331 56067	741	73336

  [START_REF] Curl | On the formation of the fullerenes[END_REF] fractions and only include those conformers that fall near the regular positions on this 2D grid, within the grid resolution for each system, namely 10×10 for C 24 and C 42 , and 8×8 for C 60 .Structures selected for C 24 , C 42 and C 60 are all depicted onto the corresponding maps in Fig.2. This sample of the totally available set of conformers amounts to about 2.4%

		C 24		C 42	C 60
		1		1	1
		0.6 0.8	set testing sample	0.6 0.8	0.6 0.8
	β	0.4		0.4	0.4
		0.2		0.2	0.2
		0		0	0
		0 0.2 0.4 0.6 0.8 1 %sp 2	0 0.2 0.4 0.6 0.8 1 %sp 2	0 0.2 0.4 0.6 0.8 1 %sp 2

Acknowledgments

The authors gratefully acknowledge financial support by the Agence Nationale de la Recherche (ANR) Grant No. ANR-16-CE29-0025, and the GDR EMIE 3533.

Supporting Information Available

Statistical convergence of the global IR spectra for C 24 and C 42 families. Interpolated spectra for the families of C 33 with isomers described at the DFTB level. Interpolated spectra for the families of C 33 and C 52 with isomers described at the force field level. Iterative label spreading classification for C 60 isomers, tagged by their asphericity β. Optimized metrics for the GR and IDW kernels.


Adrià Descals
email: a.descals@creaf.uab.cat

Aleixandre Verger
email: verger@creaf.uab.cat

Iolanda Filella

Dennis Baldocchi
email: baldocchi@berkeley.edu

Ivan A Janssens
email: ivan.janssens@uantwerpen.be

Yongshuo H Fu
email: yongshuo.fu@uantwerpen.be

Shilong Piao
email: slpiao@pku.edu.cn

Marc Peaucelle

Philippe Ciais
email: philippe.ciais@lsce.ipsl.fr

Josep Peñuelas
email: josep.penuelas@uab.cat

Soil thawing regulates the spring growth onset in tundra and alpine biomes

Keywords: Boreal forests, deciduous forests, plant phenology, spring onset, thawing, tundra. 1. Introduction

Soil temperature remains isothermal at 0 °C and water shifts to a liquid phase during soil thawing. Vegetation may receive this process as a signal and a key to restore physiological activity. We aimed to show the relationship between the timing of soil thawing and the spring growth onset. We estimated the delay between the soil thawing and the spring growth onset in 78 sites of the FLUXNET network. We built a soil thawing map derived from modelling for the northern hemisphere and related it to the greenness onset estimated with satellite imagery. Spring onset estimated with GPP time series occurred shortly after soil surface thawing in tundra (1.1 ± 3.5 days) and alpine grasslands (16.6 ± 5.8 days). The association was weaker for deciduous forests (40.3 ± 4.2 days), especially where soils freeze infrequently. Needleaved forests tended to start the growing season before the end of thawing (-17.4 ± 3.6 days), although observations from remote sensing (MODIS Land Cover Dynamics) indicated that the onset of greenness started after the thawing period (26.8 ± 3.2 days). This study highlights the role of soil temperature at the spring growth onset at high latitudes. Soil thawing becomes less relevant in temperate forests, where soil is occasionally frozen and other climate factors become more important.

Phenology is defined as the timing of periodic life-cycle events [START_REF] Schwartz | Phenology: An integrative environmental science[END_REF]. Phenology has recently gained attention as an indicator of global warming and more broadly for its link with global environmental change [START_REF] Badeck | Responses of spring phenology to climate change[END_REF][START_REF] Cleland | Shifting plant phenology in response to global change[END_REF][START_REF] Peñuelas | Complex spatiotemporal phenological shifts as a response to rainfall changes[END_REF][START_REF] Peñuelas | Responses to a warming world[END_REF]. Plant phenology is primarily modulated by the seasonal variation in climate, which involves a stage of dormancy when the climate is adverse and a period of growth when conditions are favorable for vegetation activity. Temperature, water availability, and day length have been proposed as the main environmental drivers that constrain vegetation activity and regulate plant phenology, acting eventually as confounded restraints for vegetation growth [START_REF] Chuine | Process-Based Models of Phenology for Plants and Animals[END_REF][START_REF] Jolly | A generalized, bioclimatic index to predict foliar phenology in response to climate[END_REF][START_REF] Kramer | The Annual Cycle of Development of Trees and Process-Based Modelling of Growth to Scale Up From the Tree To the Stand[END_REF]. Temperature is the main climatic factor regulating plant phenology for the onset of vegetation in the northern high latitudes [START_REF] Schwartz | Phenology: An integrative environmental science[END_REF]. Evidence for the interaction of temperature with other limiting factors, e.g. photoperiod and water availability, though, is under discussion [START_REF] Fu | Daylength helps temperate deciduous trees to leaf-out at the optimal time[END_REF][START_REF] Körner | Phenology Under Global Warming[END_REF][START_REF] Peñuelas | Complex spatiotemporal phenological shifts as a response to rainfall changes[END_REF] and, the choice of the climatic variables for modeling phenology is an open science question.

Air temperature has been widely used to simulate the transition from dormancy to growth in temperature-limited biomes, and jointly with day length, an indicator of photoperiod, for modeling plant phenology in the Northern Hemisphere [START_REF] Badeck | Responses of spring phenology to climate change[END_REF][START_REF] Peñuelas | Responses to a warming world[END_REF]. Some biome-specific studies, however, have suggested that soil temperature more accurately indicates the start of the growing season [START_REF] Baldocchi | Predicting the onset of net carbon uptake by deciduous forests with soil temperature and climate data: A synthesis of FLUXNET data[END_REF][START_REF] Jiang | The impacts of soil freeze/thaw dynamics on soil water transfer and spring phenology in the Tibetan Plateau[END_REF][START_REF] Lieth | Phenology and seasonality modeling[END_REF][START_REF] Semenchuk | High Arctic plant phenology is determined by snowmelt patterns but duration of phenological periods is fixed: An example of periodicity[END_REF][START_REF] Starr | The Photosynthetic Response of Alaskan Tundra Plants to Increased Season Length and Soil Warming[END_REF], although a study showed a stronger connection with air temperature than soil temperature for boreal forests [START_REF] Tanja | Air temperature triggers the recovery of evergreen boreal forest photosynthesis in spring[END_REF]. Some studies even suggested that the lengthening of the growing season in tundra during recent decades was primarily caused by a change in the seasonal thaw cycle [START_REF] Barichivich | Large-scale variations in the vegetation growing season and annual cycle of atmospheric J o u r n a l P r e -p r o o f CO2 at high northern latitudes from 1950 to 2011[END_REF] J o u r n a l P r e -p r o o f [START_REF] Kimball | Spring Thaw and Its Effect on Terrestrial Vegetation Productivity in the Western Arctic Observed from Satellite Microwave and Optical Remote Sensing[END_REF]. A soil warming experiment in subarctic grasslands, which reduced the occurrence of soil frost, revealed that the growing season lengthened across the entire range of IPCC warming projections for 2100 and reported no indications of photoperiod constraints within that warming range [START_REF] Leblans | Phenological responses of Icelandic subarctic grasslands to short-term and long-term natural soil warming[END_REF].

A frozen soil strictly limits vegetation growth and overrules any impact of eventual favorable air temperatures during spring. Thawing occurs when soil remains isothermal at 0 °C, triggering some processes that may be used as signals for the vegetation to restore physiological activity.

Water only returns to its liquid form and becomes available to roots [START_REF] Noormets | Phenology of ecosystem processes: Applications in global change research[END_REF] during and after thawing. In addition to a direct signal of soil thawing, soil temperature has been found to exert a stronger influence on plants than air temperature. For example, soil temperature has also been associated with physiological processes such as root elongation (Wang et al., 2018) and root respiration [START_REF] Lloyd | On the Temperature Dependence of Soil Respiration[END_REF]. Soil-warming experiments have found that an overall increase in soil temperature affected specific traits; leaf area, leaf expansion rates, and diameter growth in deciduous trees [START_REF] Farnsworth | Phenology and Growth of Three Temperate Forest Life Forms in Response to Artificial Soil Warming[END_REF][START_REF] Wheeler | Red hot maples: Acer rubrum first-year phenology and growth responses to soil warming[END_REF], needle starch content and photochemical efficiency in boreal forests [START_REF] Repo | The effect of soil temperature on the bud phenology, chlorophyll fluorescence, carbohydrate content and cold hardiness of Norway spruce seedlings[END_REF], and flowering frequency in tundra [START_REF] Khorsand | Plant phenological responses to a long-term experimental extension of growing season and soil warming in the tussock tundra of Alaska[END_REF].

These results suggest that soil temperature may play a more deterministic limiting role than air temperature at the onset of spring at high latitudes. We tested the hypothesis that soil thawing regulates spring growth onset at high latitudes but becomes less relevant where soil is only sporadically frozen.

Materials and Methods

J o u r n a l P r e -p r o o f

Phenology derived from in situ and satellite data

We used the same methodology as [START_REF] Zhang | Monitoring vegetation phenology using MODIS[END_REF] to estimate the start of season at the FLUXNET sites. We took the daily mean gross primary productivity (GPP) measurements (daytime partitioning method) and adjusted a logistic function for each year for the period January 1st to July 1st. The start of season at the FLUXNET sites (GPP SoS) was the day when the second derivative of the logistic function was maximal.

We used the MCD12Q2 V5 Land Cover Dynamics product (MODIS SoS) for the period 2001-2014. MCD12Q2 phenological metrics are estimated yearly, with a logistic function fitted over the MODIS EVI time series [START_REF] Zhang | Monitoring vegetation phenology using MODIS[END_REF]. We used the layer 'Onset Greenness Increase', which represents the day of the year on which the positive curvature of the fitted logistic function is maximal.

Estimation of the start and end of thawing

We defined the end of thawing (EoT) as the last day before solstice with daily mean soil temperature <0.5 °C and defined the start of thawing (SoT) as the closest day before EoT with a temperature below -0.5 °C. We estimated the EoT and SoT using the in situ records of soil temperature at the FLUXNET sites in the Profile 1 (TS_1_1_1), which corresponds to soil temperature measured at the shallowest layer, with a different depth depending on the FLUXNET site but commonly 2-5 centimeters. The EoT was also estimated using the layer of soil temperature at 10-centimeter depth available in the Global Land Data Assimilation System (GLDAS) dataset (See section 2.3.). This resulted in global maps of EoT at 0.25° resolution. We used a threshold of 0.5 °C instead of 0 °C to avoid false-positive estimates of EoT due to noise and slight variations in the soil temperature time series during thawing.

J o u r n a l P r e -p r o o f

Global climatic data set

We used the GLDAS version 2.1 as a global climatic data set for the daily representation of soil temperature at a global scale. GLDAS is a land-surface modeling system developed by the Goddard Space Flight Center of the National Aeronautics and Space Administration and by the National Centers for Environmental Prediction of the National Oceanic and Atmospheric Administration. The version 2.1 provides data for the period 2000 up until the present time.

Ground measurements (meteorological data), such as soil temperature, and satellite observations are ingested in order to constrain the land-surface model and restrain unrealistic model states.

GLDAS offers 3-hourly estimates of earth-surface variables such as soil temperature at different depths. We used the soil temperature product at 10-centimeter depth with a spatial resolution of 0.25°.

Resizing and filtering MODIS SoS

The spatial resolutions of GLDAS EoT and MODIS SoS differed, so we resized the resolution of MODIS SoS (500 m) to the resolution of GLDAS EoT (0.25°) to make the comparison between the two datasets consistent. We first masked the yearly MODIS SoS images with the MODIS IGBP Land Cover (MCD12Q1) product. MODIS SoS was masked separately with the IGBP classes 1) deciduous broadleaved forests (deciduous forests), 2) evergreen needleleaved forests (needleleaved forests), 3) shrubland, 4) open shrubland, and 5) grassland. For the tundra biome, we only considered the pixels classified as shrubland, open shrubland, and grassland at >50° latitude. Finally, the masked MODIS SoS images were resized to 0.25° using the mean as the resampling operation. This resulted in six sets of MODIS SoS for each land cover for the period 2001-2014.

J o u r n a l P r e -p r o o f

Experimental setup and statistical analyses

We analyzed 78 sites of the FLUXNET network and compared the timing of spring onset of ecosystem GPP with thawing date. We used all sites north of 30°N latitude that are located in tundra (13 sites) and alpine (4 sites) climates, and the sites categorized in the IGBP classification as deciduous broadleaved forests (DBF) (18 sites in temperate and 1 site boreal climate zones) or evergreen needleaved forests (ENF) (25 sites in temperate and 17 sites boreal climate zones).

Table S1 shows the coordinates and biome of the selected sites. We excluded the sites without in situ soil temperature records.

We investigated at which soil temperature the GPP was significantly positive in each biome. We grouped the half-hourly GPP observations in intervals of 0.5°C of soil temperature. We applied a two-tailed Student's t-test at each interval of soil temperature with a significance of 95%. We considered that GPP was positive at a certain temperature when the lower GPP confidence interval was greater than 0.5 mol CO 2 m -2 s -1 .

The SoT and EoT estimated with the FLUXNET records were compared with two independent observations of the start of the growing season (SoS): 1) SoS estimated with daily in situ observations of the FLUXNET GPP time series (GPP SoS) and 2) SoS extracted from remotesensing sources (MODIS SoS). We applied the two-tailed Student's t-test with a significance of 95% to find whether the day of the EoT and the spring onset (GPP SoS and MODIS SoS) were significantly different. We reported the mean difference, confidence interval, and coefficient of determination (R 2 ) between the end of thawing and the spring onset.

We also compared the spring onset observed with satellite derived data (MODIS SoS) with the thawing period estimated with the GLDAS. The MODIS SoS was filtered and resized to the J o u r n a l P r e -p r o o f resolution of the GLDAS to make the datasets comparable. Similarly to the comparison at the FLUXNET sites, we reported the R 2 and the mean difference with a confidence interval of 95% between the EoT and the MODIS SoS. The analysis covers a time period from 2000 to 2015 for the comparison between FLUXNET and MODIS, and from 2001 to 2017 for the comparison between the GLDAS and MODIS.

Results

The growing seasons were primarily driven by temperature in the 78 FLUXNET sites; GPP and air temperature had similar intra-seasonal patterns. Gross primary productivity at these sites was highest in summer, when the air and soil temperatures were also highest. The Italian site IT-Ro1 was the only exception, where GPP decreased during summer due to water stress. A representation of air and soil temperature, GPP and spring onset is shown in Fig. S1 for five FLUXNET sites. GPP was insignificant with frozen soil in all biomes (Fig. 1). The lowest soil temperature at which GPP was significantly higher than 0.5 mol CO 2 m -2 s -1 was 1°C in tundra, -0.5°C in alpine grasslands, -1 °C in needleaved forests, and 4.5 °C in deciduous forests.

Spring onset occurred during thawing (41% of observations for GPP SoS and 6% of the observations for MODIS SoS) or after EoT (51% of observations for GPP SoS and 59% of the observations for MODIS SoS) (Fig. 2). Only 8% of the observations occurred when soil was frozen in GPP SoS, primarily in needleaved forests. Spring onset for tundra and alpine grasslands tended to start shortly after the EoT (1.1 ± 3.5 d for GPP SoS and 18.5 ± 3.0 d for MODIS SoS for tundra, 16.6 ± 5.8 d for GPP SoS and 32.9 ± 8.3 d for MODIS SoS for alpine grasslands).

Spring onset tended to start with a longer delay after EoT in deciduous forests with seasonally frozen soil (40. The association between the GPP SoS and the EoT SoS is different in deciduous than in needleleaved forests (Fig. 3a). Spring onset occurs with positive temperatures in deciduous forests, while needleleaved forest starts their growing season before the soil has completely thawed. Tundra and alpine grasslands show the highest agreement of spring onset with the EoT (RMSE = 9.03 in GPP SoS and RMSE = 20.45 in MODIS SoS). FLUXNET sites in the boreal biome also showed a better agreement than the sites in the temperate biome in both GPP SoS and MODIS SoS (Fig. S2). In addition, soil did not freeze in some sites of the temperate biome, which explains the EoT values equal to 0 in Fig. 3 while deciduous forests showed the highest delay (34.24 ± 18.59 d). The high delay in deciduous forests can be observed in the map showing the difference mean MODIS SoS -mean EoT (Fig. 5) at the most southerly latitudes of the deciduous forest distribution (eastern North America, central Europe, and northern China). We observed a high interannual variability in the MODIS SoS and EoT also at low latitudes (Fig. S3). Contrarily to the results in the FLUXNET sites, the MODIS SoS in needleleaved forests show a similar pattern as deciduous forests and do not present a negative delay compared to the EoT (20.03 ± 14.15 d).

Discussion

J o u r n a l P r e -p r o o f

Spring onset of vegetation activity and soil thawing were closely associated in tundra and alpine grasslands, and weakly associated where soil froze intermittently. Water remains inaccessible to plants in frozen soils, inhibiting vegetation activity. This process is particularly relevant in regions with deep frost, where the vegetation has maximal root depth <0.5 m in tundra or <2 m in boreal forests [START_REF] Canadell | Maximum rooting depth of vegetation types at the global scale[END_REF]. Likewise, most root biomass in boreal and temperate forests is in shallow soil layers, where most mineralization processes release the greatest amount of nutrients. Frozen soil, even if only the first few centimeters, may therefore imply inhibition of vegetation activity, as suggested by the results from the FLUXNET towers, where the GPP is insignificant when soil temperatures are negative (Fig. 1). In addition, most of the observations of spring onset, both MODIS SoS and GPP SoS, coincided with positive or near freezing soil temperatures (Fig. 2).

Water gradually shifts from the solid to the liquid form during thawing, which plants may use as a signal to start their activity. Our results, however, indicated that the time between EoT and spring onset depended on the biome. These differences can be attributed to different forcing requirements for which vegetation has different adaptive histories. Deciduous forests seem to require greater heat accumulation than tundra and needleleaved forests to initiate spring growth.

Many studies have suggested that the high thermal forcing required by deciduous trees, which delays the vegetation onset date after the EoT, is an adaptation to prevent frost damage in late spring [START_REF] Fu | Understory CO2, sensible heat, and latent heat fluxes in a black spruce forest in interior Alaska[END_REF][START_REF] J O U R N A L P R E -P R O O F Laube | Chilling outweighs photoperiod in preventing precocious spring development[END_REF]. In addition, deciduous forests have adapted their high thermal forcing depending on the background climate [START_REF] Peaucelle | Spatial variance of spring phenology in temperate deciduous forests is constrained by background climatic conditions[END_REF], with lower heat requirements in colder climates. This would explain the high delay between the spring onset and EoT in deciduous forest in the lowest latitudes of the study, while the only FLUXNET site with deciduous forests located in the boreal zone shows a close association between the spring onset J o u r n a l P r e -p r o o f and the thawing period. Similarly to the site with deciduous forest in the boreal biome, tundra and alpine grasslands have a low heat requirements, very likely because their SoS also occurs very late in the season, when light availability is already high, and these plants may therefore have reduced their heat requirement to maximize profit from the high light conditions in late spring.

Our results show that needleaved forests can release from dormancy and start their gross primary production during the thawing period with near freezing soil temperatures, which corroborates the photosynthetic acclimation to cold in boreal needleaved forests [START_REF] Oquist | Photosynthesis of overwintering evergreen plants[END_REF].

However, the results obtained from remote sensing show that the onset in greenness in needleaved forests occurs much closer to the end of the thawing period. This suggests that evergreen needleaved trees start the growing season when temperatures become favorable using the needles preserved during the dormant period, but needle emergence occurs much later than the start of the photosynthetic activity. The delay between GPP SoS and onset of greenness was also observed with remotely-sensed data in previous studies [START_REF] Melaas | Using FLUXNET data to improve models of springtime vegetation activity onset in forest ecosystems[END_REF][START_REF] Walther | Satellite chlorophyll fluorescence measurements reveal large-scale decoupling of photosynthesis and greenness dynamics in boreal evergreen forests[END_REF]). An alternative explanation for the early GPP SoS in boreal forests is that understory might have started the photosynthetic activity before the trees in the overstory and, thus, vegetation in the understory layer accounted for the carbon fluxes before the thawing period. This would be supported by the findings of Ikawa et al., 2015, which showed the contributions of the understory in boreal ecosystems to the carbon and energy balances.

The hypothesis that tundra has adapted to the long period with frozen soils by reducing its heat requirement implies that spring onset is highly responsive to soil thawing. This hypothesis is supported by studies with long time series of snow melt and by soil-warming experiments reporting that spring onset is advanced by an increasingly earlier soil thawing (Barichivich et al., J o u r n a l P r e -p r o o f 2013; [START_REF] Kimball | Spring Thaw and Its Effect on Terrestrial Vegetation Productivity in the Western Arctic Observed from Satellite Microwave and Optical Remote Sensing[END_REF][START_REF] Leblans | Phenological responses of Icelandic subarctic grasslands to short-term and long-term natural soil warming[END_REF]. Spring onset in arctic and subarctic vegetation may thus continue to advance as soil thawing occurs earlier with climate warming. The possible restriction in vegetation growth because of insufficient levels of incoming radiation or hours of sunshine is, however, still undetermined.

Our results suggest that the role of soil temperature is secondary at southerly latitudes for deciduous and needleaved forests. Sites with favorable soil and air temperatures throughout the year may be constrained by photoperiod [START_REF] Zohner | Day length unlikely to constrain climate-driven shifts in leaf-out times of northern woody plants[END_REF], by chilling requirements [START_REF] Chuine | Warming, Photoperiods, and Tree Phenology[END_REF], or by both [START_REF] Fu | Daylength helps temperate deciduous trees to leaf-out at the optimal time[END_REF]. These regions without soil temperature regulation had the largest differences between MODIS SoS and EoT, as observed in the difference map in Figure 5.

The differences observed between the EoT and the GPP SoS at the FLUXNET towers with needleleaved forests cannot be depicted with the MODIS SoS and, thus, are not reflected in the map in Figure 5.

In summary, we demonstrate that spring onset rarely occurs with negative soil temperatures, indicating that frozen soil acted as a major constraint on vegetation activity. We also demonstrate that vegetation activity began shortly after thawing in tundra, after soil temperatures at 2centimeter depth surpassed the 0 °C threshold. This correlation became weaker for lower latitudes. Needleleaved forests can start the gross primary production during the thawing period, while deciduous forests delayed the spring onset especially at more southerly latitudes. These results thus provide insight in the onset of spring that can be used for improving phenological models, especially in tundra and boreal forests where soil thawing plays a major role in controlling the start of vegetation activity. Further studies may take advantage of recent advances in mapping soil thawing using microwave remote sensing to estimate the start of the growing season at northern latitudes.

J o u r n a l P r e -p r o o f

Figure 1 | Gross primary productivity (GPP) in relation to soil temperature in the FLUXNET sites. The x-axis shows soil temperature measured at the shallowest layer. The distribution of GPP is presented with the 50 th percentile (continuous line), 75 th percentile (dashed line), and 90 th percentile (dotted line). The red line shows the lowest temperature in which GPP is significantly higher than 0°C. The temporal resolution of the GPP and temperature data is half-hourly.

J o u r n a l P r e -p r o o f

  3 ± 4.2 d for GPP SoS and 46.4 ± 3.8 d for MODIS SoS), while needleleaved J o u r n a l P r e -p r o o f Journal Pre-proof forests showed discrepant results between the two sources of data (-17.3 ± 3.6 d for GPP SoS and 26.8 ± 3.2 d for MODIS SoS).

  and Fig. S2. The comparison of MODIS SoS with the global estimates of EoT (Fig. 4) led to similar results as the observations at the FLUXNET sites. The onset of greenness estimated with MODIS occurred near EoT in the land cover types that represent tundra (8.63 ± 7.96 d in shrubland IGBP class, 16.97 ± 7.72 d in open shrubland IGBP class, and 18.80 ± 14.55 d in grassland IGBP class),

Figure 2 |Figure 3 |

 23 Figure 2 | Boxplots of the normalized delay ((a) and (c)) and absolute delay ((b) and (d)) of spring onset after the end of thawing (EoT) for tundra, alpine grasslands, needleaved forests, and deciduous forests. The absolute delay is the number of days between spring onset and EoT, and the normalized delay represents spring onset relative to thawing, with start of thawing (SoT) = 0 and EoT = 1. Spring onset is presented as FLUXNET gross primary productivity (GPP) start of season (SoS) in (a) and (b) and as MODIS SoS ('MidGreenup1' band in MCD12Q2v6 product) in (c) and (d). The central mark of the boxplot indicates the median, the bottom and top edges indicate the 25 th and 75 th percentiles. The whiskers cover until the most extreme data points not considered outliers, and outliers are plotted with the '+' symbol. J o u r n a l P r e -p r o o f

Figure 4 |

 4 Figure 4 | Relationships between MODIS start of season (SoS) and end of thawing (EoT) estimated using the GLDAS dataset. Each panel represents MODIS SoS ('MidGreenup1' band in MCD12Q2v6 product) masked with the MODIS IGBP Land Cover product for cover types (a) grassland, (b) open shrubland, (c) shrubland, (b) evergreen needleleaved forests, and (e) deciduous broadleaved forests. The land cover types grassland, open shrubland, and shrubland for latitudes higher than 50° represent predominantly the tundra biome. The maximum density (md) varies depending on the land cover.

Figure 5 |

 5 Figure 5 | Maps of (a) mean end of thawing (EoT), (b) mean MODIS start of season (SoS), and (c) mean difference (MODIS SoS -EoT) for the period 2001-2017. EoT was estimated using the GLDAS data set as the last day of the dormant period with soil temperatures <0.5 °C. MODIS SoS is the 'MidGreenup1' band in MCD12Q2v6 product. We resized MODIS SoS to the spatial resolution of the GLDAS dataset. The map projection is the Lambert Azimuthal Equal-Area.

Figure 2 Figure 3

 23 Figure 2

Figure 4

 4 Figure4

Figure 5

 5 Figure5

Acknowledgements

The authors would like to acknowledge the financial support from the European Research Council Synergy grant ERC-SyG-2013-610028 IMBALANCE-P, the Spanish Government grant PID2019-110521GB-I00 and the Catalan Government grant SGR 2017-1005.

Data availability statement

The data that support the findings of the study are openly available from FLUXNET at [https://fluxnet.fluxdata.org/], GLDAS version 2.1 at [https://ldas.gsfc.nasa.gov/gldas/], MODIS Land Cover Type MMCD12Q1v6 at [https://doi.org/10.5067/MODIS/MCD12Q1.006], and MODIS Land Cover Dynamics MCD12Q2v6 at [https://doi.org/10.5067/MODIS/MCD12Q2.006].

Code availability

The code is available from the authors upon request.

J o u r n a l P r e -p r o o f

Journal Pre-proof Author contributions: AD and JP conceived the research idea. AD, AV, and JP designed the study. AD performed the analyses and wrote the first version of the manuscript. AD, AV, IF, DB, IJ, YF, SP, MP, PC, and JP contributed to the interpretation of the results and to revisions of the manuscript.

J o u r n a l P r e -p r o o f

Declaration of interests

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

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:  We estimated the delay between the soil thawing and the spring growth onset in 78 sites of the FLUXNET network and with remote sensing data.

 Soil thawing plays a role at the onset of spring growth in alpine biomes and at high latitudes, and the further north the vegetation grows the stronger is the role.

J o u r n a l P r e -p r o o f