Direct Observation of Hierarchic Molecular Interactions Critical to Biogenic Aerosol Formation

Small clusters consisting of sulfuric acid/bisulfate and oxidized organics have been identified in both aerosol field measurements and laboratory experiments, and their formation is suggested to be the rate-limiting step in the formation of new particles. However, the underlying mechanism for cluster formation is still largely unclear. Here we show, through an integrated negative ion photoelectron spectroscopy and quantum chemical study on a series of (HSO4−)(organic molecule) surrogate binary clusters, that the functional groups are more important in determining the extent of the enhanced role of the organics in aerosol formation process than the average carbon oxidation states or O/C ratios. This extent is quantified explicitly for specific functional groups, revealing highly hierarchic intermolecular interactions critical to aerosol formation. Born–Oppenheimer molecular dynamics simulations are employed to probe the water-binding abilities of these clusters under ambient conditions, and their statistical hydrogen-bonding networks. Studying the early stages of aerosol formation is a challenge in physical and environmental chemistry. Here photoelectron spectroscopy, quantum chemical calculations, and molecular dynamics simulations quantify how specific functional group interactions stabilize clusters of bisulfate anions and organic molecules.

A tmospheric aerosol particles can influence the cloud properties and Earth's radiative balance both directly by interacting with the solar radiation and indirectly by acting as cloud condensation nuclei (CCN) 1 , thus impacting the climate change, air quality, and human health [2][3][4] . New particle formation (NPF) has been recognized as an important pathway contributing to atmospheric aerosols 5 . It is generally accepted that NPF proceeds via two distinct stages-nucleation from gaseous precursors to form critical nuclei in the sub-2-nm size, and subsequent growth to larger size (>2-3 nm) clusters and particles 5,6 . Recent field measurements and laboratory experiments have demonstrated that both N-containing compounds like NH 3 and amines [6][7][8] , and the oxidation products of volatile organic compounds (VOCs) [9][10][11] can greatly enhance the nucleation and initial growth of atmospheric aerosol particles, and have further identified abundant small clusters of H 2 SO 4 /HSO 4 − and NH 3 / amines or oxidized organic molecules existed in the~2 nm size aerosol particles. Most of the molecular clusters involved in the NPF process in the atmosphere are electrically neutral 5 , although the importance of ion-induced nucleation (IIN) [12][13][14][15] has been progressively recognized. We are, in particular, interested in elucidating the roles of oxidized organics in the initial aerosol nucleation, as recent kinetic measurements have indicated that small clusters consisting of one sulfuric acid molecule and one or two organic molecules from oxidation products of biogenic emissions are involved in the rate-limiting steps [9][10][11]16 , strongly suggesting potentially important roles of such binary clusters in the initial NPF processes.
To date, no systematic direct experimental characterizations on the structures, charge locations, interaction strengths, thermodynamics, and kinetic properties of these small clusters beyond the chemical compositions have been reported. This lack of information precludes a fundamental molecular-level understanding of the underlying mechanisms for initial nucleation in NPF processes, and hampers the efforts in developing reliable and predictive nucleation and climate models.
Here we report a combined spectroscopic and theoretical investigation on a series of (HSO 4 − )(organic molecule) binary clusters that can be regarded as surrogate clusters involved in the formation of critical nuclei. Negatively charged clusters instead of the corresponding (H 2 SO 4 )(VOC) neutral clusters are chosen, enabling precisely selecting and controlling the sizes and compositions of the clusters of interest. The fundamental molecular interactions derived from these charged clusters should be helpful in elucidating the critical forces in the neutral clusters as well, as recent studies suggested similar trends in structural evolution and interaction for clusters involving VOCs and sulfuric acids in both neutral and negatively charged states 11,17 . Because α-pinene is one of the most abundant emissions from natural plants and accounts for nearly 50% of the global monoterpene budget 18 , we are focusing on oxidized organics produced in different generations of oxidation processes from α-pinene, including pinanediol, cis-pinonic acid (cPA), pinic acid (PA), and 3-methyl-1,2,3butanetricarboxylic acid (MBTCA). Extensive research has identified the critical roles of these oxidation products in contributing to the NPF events over forests, urban, and rural areas [18][19][20][21][22][23] . We also select a variety of oxidation products from anthropogenic emissions, such as aromatic acids, monocarboxylic acids (MCAs), and dicarboxylic acids (DCAs), for comparative studies. These compounds derived from both biogenic and anthropogenic emissions all together span a wide range of carbon numbers from 1 to 10, with their average carbon oxidation states OSc À Á 24 from −2 to +3 and O/C ratios from 0 to 2, and contain multifunctional groups, such as hydroxyl, carbonyl, and carboxylic groups. Therefore, we anticipate the organic compounds studied here to be important surrogate molecules representing a significant portion of organics with distinctly different physiochemical properties in the atmosphere.

Results
The VOCs studied in this work. A total of twelve different organic compounds are studied. Their distributions in OSc or O/C ratios vs. the number of carbon atoms (n c ) spaces are presented in Fig. 1 25,26 . The oxidized α-pinene product with the highest oxidation state, MBTCA, falls into the extremely lowvolatility organic compound (ELVOC), and thus can also be considered as a ELVOC 27 . Other ELVOCs containing several peroxy acid groups have not been included in the current study due to the unavailability of these chemicals 11 . Figure 1 reveals that the distributions of these compounds in either the OSc-n c or O/Cn c spaces are overall similar, i.e., the OSc and O/C of these compounds in general decrease with increasing n c . However, the volatility of the compounds varies appreciably for each given n c = This can also be seen quantitatively from the plots of logC* in OSc-n c or O/C-n c spaces (Supplementary Figure 1). To facilitate the evaluation of the experimental data, we grouped these VOCs into five categories based on their functional groups-alcohol (pinanediol; IVOC), MCAs (p-toluic acid, benzoic acid, formic acid, acetic acid, and pentanoic acid; IVOC), ketone-carboxylic acid (cis-pinonic acid; SVOC), DCAs (pinic acid, oxalic acid, and succinic acid; SVOC), and tricarboxylic acid (MBTCA; LVOC).
NIPE spectra of (HSO 4 − )(VOC) binary clusters. The gaseous 1:1 (HSO 4 − )(VOC) binary clusters were generated via electrospray ionization and characterized by cryogenic negative ion photoelectron spectroscopy (NIPES) 28 . The experimental threshold and vertical detachment energies, i.e., TDEs and VDEs, are determined from the onset and peak maximum of each spectrum, respectively, and are summarized in Supplementary Table 2. Supplementary Figure 2 compares five NIPE spectra of (HSO 4 − )(MCA) (pentanoic, p-toluic, benzoic, acetic, and formic acids) binary clusters, as well as three spectra of (HSO 4 − )(DCA) (pinic, succinic, and oxalic acids). It can be seen that all (HSO 4 − ) (MCA) clusters have a similar electron binding energy increase (ΔEBE) of 1.05 ± 0.05 eV relative to that of HSO 4 − , while those of (HSO 4 − )(DCA) clusters have a much higher ΔEBE of 1.5 ± 0.1 eV. The (HSO 4 − )(succinic acid) cluster has an extra~0.2 eV ΔEBE increase compared to the other two dicarboxylic acids, which may be due to the formation of a strong intramolecular hydrogen bond in the succinic acid part of the cluster 17 . Here it should be mentioned that the ΔEBE increase corresponds to the energy difference between HSO 4 − and HSO 4 • clustering with the organic compound, well reflecting the lower limits of the intermolecular binding or interaction strength in the anionic clusters (Supplementary Figure 3) 17,29 . Thus, the larger the ΔEBE is, the stronger intermolecular binding of the organics with bisulfate in the corresponding anionic clusters will be. The experimental data shown in Supplementary Figure 2 suggests that when the compounds have the same functional groups, they have similar binding/interaction strengths with the bisulfate ion in forming clusters even though they have very different O/C ratios and oxidation states (O/C ratio spans from 0.2 to 2.0, 0.4 to 2.0; and OSc from −1.2 to +2.0 and −0.67 to +3.0 for the five MCAs, and three DCAs, respectively. Supplementary Table 1 . This implies that the functional group is much more important than the O/C ratio or OSc when considering the binding ability of a specific organic molecule interacting with bisulfate. Therefore, we selected ptoluic acid and pinic acid as representative molecules for MCAs and DCAs, and compared their spectra with other biogenic organic compounds that possess different functional groups from MCA and DCA, but with similar carbon numbers. Figure 2 presents the 20 K NIPE spectra measured for (HSO 4 − ) (VOC) binary clusters with the VOCs having different functional groups but similar n c . It can be seen that the ΔEBEs generally increase as O/C ratio or OSc becomes large; specifically they are 0.85, 1.25, 1.50, 1.10, and 1.85 eV for VOC = pinanediol, cispinonic acid, pinic acid, p-toluic acid, and MBTCA, respectively. (See Supplementary Table 2 for detailed experimental values and  Supplementary Figures 8-9 for the ΔEBE plots with O/C ratio or OSc). This is actually consistent with the change of the VOC functional groups along this series of compounds, as more extensive functional groups are expected to be present in the molecules with increasing O/C ratio or OSc when they have similar carbon numbers. One exception is that the ΔEBE of (HSO 4 − )(p-toluic acid) lies between those of (HSO 4 − )(pinanediol) and (HSO 4 − )(cis-pinonic acid), although p-toluic acid has a much higher OSc than cis-pinonic acid (Supplementary Figure 9). However, the ΔEBE trend in Fig. 2 follows the O/C ratio as also shown in Supplementary Figure 8. This may be due to the fact that the extra negative charge is more delocalized over the (HSO 4 − )(p-toluic acid) cluster than in (HSO 4 − )(cis-pinonic acid), leading to this electron binding energy variation from the OSc perspective on one hand, and on the other hand suggests that cautions must be practiced when considering subtle changes in correlating the enhanced roles of organics with O/C ratio or OSc. The data shown in Fig. 2 and Supplementary Figure 2 suggests that the specific types of functional groups of organics play overwhelming roles in determining intermolecular interaction strength for formation of these small clusters, e.g., the carboxylic group interacts much stronger with bisulfate than the hydroxyl group (for example, ΔEBE = 1.85 eV for MBTCA with three COOH groups, and 1.50 eV for pinic acid with two COOH groups over 0.85 eV for pinanediol with two OH groups). This enormous difference of intermolecular interactions allows us to track the specific molecular properties that are critical during the NPF processes by ranking highly hierarchic intermolecular interactions according to the different functional groups. This ranking is quite useful to predict whether the interaction between bisulfate (sulfuric acid) and one organic molecule is thermodynamically more favorable than the other one, which, in addition, will be complementary to the classical nucleation theory (CNT) that relies on the macroscopic properties of precursor molecules to study the nucleation mechanisms. Next, we will turn to analyses of the optimized cluster structures and energetics to capture the molecular characteristic signatures of the VOCs with different functional groups upon interacting with the bisulfate ion.
Structures of (HSO 4 − )(VOC) binary clusters. The optimized low-lying structures of (HSO 4 − )(VOC) binary clusters with different functional groups are presented in Fig. 3 Figure 10), involving both proton-transferred (PT) and non-protontransferred (nPT) binding schemes with nearly degenerate stabilities (<1.6 kcal mol −1 ). For example, cis-pinonic acid forms two low-lying isomers when interacting with bisulfate-a nPT structure (pinonic acid)(HSO 4 − ), and a PT one with the proton transferred from the organic acid to bisulfate (pinonic acid-H + ) − (H 2 SO 4 ) that is in fact more stable by 1.60 kcal mol −1 . The calculated vertical detachment energy (VDE) of 6.30 eV of the nPT structure agrees well with the maximum of the first spectral band, while the calculated VDE of the PT isomer (6.63 eV) aligns with the second spectral feature. Therefore, proton transfer from monocarboxylic acids to bisulfate in forming (carboxylate)(sulfuric acid) clusters should be quite general. Considering sulfuric acid as being a much stronger acid compared to monocarboxylic acids, the formation of PT isomers is counter intuitive to gasphase acidity predictions. The underlying mechanism has been recently studied in detail in the (formic acid)(bisulfate) case 30 .
Dicarboxylic acids, however, are found much more favorable in transferring one proton to bisulfate, leading to (sulfuric acid) (carboxylate) double hydrogen bonds, along with formation of an extra intramolecular hydrogen bond in the dicarboxylic acid moiety (Supplementary Figure 11). This is because the formation of intramolecular hydrogen bond slightly increases the acidity of the corresponding organics, thus facilitating the proton transfer from organics to bisulfate ion 17 . Further analysis reveals that the strength of the intramolecular hydrogen bond in these dicarboxylic acids increases in the order of pinic acid < oxalic acid < succinic acid (Supplementary Figure 12), which is consistent with the correspondingly observed and calculated ΔEBE increases (Supplementary Table 2) and QTAIM analysis (Supplementary Table 3, Supplementary Note 1). A nPT isomer of (HSO 4 − )(pinic acid) with triple hydrogen bonds has been found to be slightly more stable by~0.6 kcal mol −1 . This structure has relatively lower calculated VDE compared to the experimental value, and may contribute to the rising edge in the threshold region of the experimental spectrum. The PT structure for the pinic acid case should be the dominant one observed in the experiments and contributes to the main spectral features. For MBTCA, the PT and nPT structures are almost degenerate in energy, only differing by 0.24 kcal mol −1 . In summary, both types of structures, i.e., the nPT−(HSO 4 − )(organic acid), and PT −(H 2 SO 4 )(carboxylate) clusters, coexist in the experiments and contribute to the experimental spectra. However, only the nPT type makes significant contribution to the onset band in each spectrum, as the VDE of this structural motif is always lower than the corresponding PT one (Fig. 3 and Supplementary Figures 10-13), with one exception, that is MBTCA, for which the calculated VDE of the nPT form is slightly higher than the PT one, rendering both structures contributing to the threshold rising edge in the spectrum (Fig. 2). Further experimental evidence and insights on the structures and proton transfer mechanism can be complementally obtained by measuring their infrared photodissociation (IRPD) spectra [31][32][33][34] , and we hope that the publication of the current study will stimulate such efforts.
Thermochemical properties of cluster formation. The thermochemical parameters, e.g., binding energy (BE), enthalpy, entropy, and Gibbs free energy (ΔG), associated with the (HSO 4 − )(VOC) cluster formation under ambient conditions (1 atm and 298.15 K) for the characteristic structures shown in Fig. 3   same functional groups have similar BE and ΔG when clustering with HSO 4 − -the higher the BE is, the more favorable the ΔG will be. The calculated BEs and ΔGs for biogenic VOCs interacting with bisulfate increase along with their oxidation pathways, from α-pinene → pinanediol → cis-pinonic and pinic acids, and → MBTCA. The calculated values of p-toluic acid are in between those of pinanediol and cis-pinonic acid, exactly following the ΔEBE trend derived from the experimental spectra (Supplementary Figures 8, 9).
In the atmosphere, the formation of stable clusters is governed by the competition between intermolecular stabilization, i.e., the thermodynamic stability of the clusters with respect to their noninteracting constitutes, and decay of the clusters due to evaporation or coagulation by pre-existing condensation sink particles 16,35 . Small clusters must be stable enough, and must remain suspended in the atmosphere for a sufficient time to interact with other atmospheric molecules in order to grow into large particles. Recent studies have found that the growth rate observed for aerosol particles smaller than 1.2-nm mobility diameter is of the order 0.2 nm h −1 ref. 6 , which requires that the evaporation rate of such clusters must be lower than 5.56 * 10 −5 s −1 . To note, Kulmala and co-workers have reported a median value of condensation sink during NPF events to be 1.3 * 10 −3 s −1 from a 11-year field measurement 36 . In order to have a direct comparison, we converted the (HSO 4 − )(VOC) cluster mass to mobility-equivalent diameters 11,37 and found that they all fall into the size range of 1.0-1.2 nm (Supplementary Table 1 Table 5), rendering them high possibilities to serve as embryo clusters to grow into large particles from a thermodynamic point of view. It should be mentioned that the kinetic factors due to curvature (Kelvin) effect and concentrations of precursor ions and molecules, which are also important in understanding the nucleation and growth of particles, are out of scope in the present work.
As reported previously 5,9 , critical nucleus of atmospheric particles is generally of a size of a few nanometers, often, including water molecules. In addition, temperature and pressure are important factors in determining the stability of these nuclei. Therefore, we carried out Born-Oppenheimer molecular dynamics (BOMD) simulations on (HSO 4 − )(VOC)(H 2 O) 10 ternary systems under ambient temperature and pressure that can better mimic the atmospheric conditions to further investigate the stability of these (HSO 4 − )(VOC) binary clusters within water clusters, to probe their dynamic behaviors and water-binding abilities, and to obtain statistic average of hydrogen-bonding network (HBN) information. Figure 4a-e shows the snapshots of the equilibrium structures of (HSO 4 − ) (VOC) binary clusters each interacting with 10 water molecules and Fig. 4f-j presents the total hydrogen bonds of clusters as well as the hydrogen bonds formed by VOCs. The VOCs with more extensive functional groups (pinic acid and MBTCA) have a higher possibility to form hydrogen bonds with HSO 4 − or water. In particular, MBCTA can form 2-3 times more hydrogen bonds with bisulfate and water than other VOCs being capable of (see the average of VOC H-bonds in Fig. 4f-j), and there are 2-3 more hydrogen bonds for the total number of hydrogen bonds formed in the (HSO 4 − )(MBTCA)(H 2 O) 10 cluster than other ternary clusters (Supplementary Table 6 and Supplementary Figure 14). The simulations further indicate that these (HSO 4 − )(VOC) (water) clusters, albeit with fluidic HBN and total hydrogen bond numbers fluctuated, do possess dynamic stability with average hydrogen bond numbers remaining constant during the entire simulation timeframe (30 ps). Interestingly, we found that the water molecules mainly accumulate on the hydrophilic carboxylic/bisulfate part, which is consistent with the finding reported by Zhang and co-workers in their aerosol chamber experiments and theoretical study on neutral cis-pinonic acid/ sulfuric acid/water systems 9 .

Discussion
Ion-induced nucleation (IIN) has been recognized as an important alternative nucleation pathway besides the neutral form in the overall aerosol formation processes 6,11-15 . Schobesberger et al. have shown that both the compositions of negatively charged clusters and their NPF formation mechanisms are similar to those of their neutral counterparts, particularly for particles containing oxidized organics and sulfuric acids 11 . Our previous study on (bisulfate/sulfuric acid) (succinic acid) clusters 17 also supported such findings, suggesting similar structural evolution and interaction trends between the neutral and negatively charged states. Several recent theoretical studies have explored the structures of neutral (H 2 SO 4 )(VOC) clusters and their formation thermochemistry [38][39][40][41][42][43][44] . Hence, it would be instructive to compare our studies on the negatively charged clusters with those on the neutral ones. In Supplementary Table 7, we summarized the Gibbs free energy changes (ΔG) for the formation of (H 2 SO 4 ) (VOC) versus (HSO 4 − )(VOC) clusters. It can be seen that the trend across different VOCs is overall in accord with each other, although the absolute thermodynamic quantities are much smaller in the neutral clusters than the negatively charged ones as we studied. In addition, it is also noted that the calculated absolute quantities of ΔG for the neutral clusters vary a lot depending on different theoretical methods employed, highlighting the importance of reliable experimental benchmarks to determine the (relative) energetics. Therefore, the detailed information on the structures, stabilities, water-binding abilities, hierarchic functional group-determined molecular interactions, and the dynamic properties of a wide range of (HSO 4 − )(VOC) binary clusters, investigated in the present work, provides underlying molecular insights to help better understanding the early stages of the aerosol nucleation, both in the neutral and negatively charged states.
In conclusion, we have studied a series of 1:1 clusters with mobility diameters of~0.9-1.2 nm formed by the bisulfate ion and the oxidation products from both biogenic and anthropogenic emissions employing NIPES, ab initio quantum chemical calculations, and molecular dynamics simulations. Both experiment and theory suggest that the exact nature of the organics (the carbon backbone) is less important than the specific functional groups when they form clusters (via hydrogen bonds) with bisulfate ion or sulfuric acid. The stabilization effect of clustering processes has been experimentally quantified for each specific functional group commonly existed in oxidized organics, and it is in good agreement with theoretical results. Our study reveals that the functional group is more important than the OSc or O/C ratio, and needs to be taken into account when considering the binding ability of a specific organic molecule interacting with bisulfate, especially when organics have different carbon numbers. The organics with multi-functional groups like MBTCA have high possibilities to interact with bisulfate ion/sulfuric acid, and the resulting clusters are extremely stable, which will further serve as a core when interacting with other species (like water) to grow into larger particles. We expect the detailed molecular-level information, including structures, energetics, thermochemistry, and dynamic properties obtained in this work, to be important to benchmark theoretical models and to help better understand the underlying NPF nucleation mechanisms involved in both neutral and ion-induced pathways.

Methods
Experimental methods. The NIPES experiments were performed using the PNNL low-temperature, magnetic-bottle time-of-flight (TOF) photoelectron spectrometer, coupled with an electrospray ionization (ESI) source and a temperaturecontrollable cryogenic ion-trap 28 . The (HSO4 -)(VOC) binary clusters were produced via electrospraying into the gas phase a 0.1 mM solution of sulfuric acid dissolved in water/acetonitrile, mixed with three times of 0.1 mM solution of VOC dissolved also in water/acetonitrile. The produced ions were accumulated and collisionally cooled in the cryogenic 3D ion trap set at 20 K for 20-100 ms to minimize the populations of high energy isomers and eliminate the extra features due to vibrational hot bands in the NIPE spectra. The cold ions were then transferred into the extraction zone of a TOF mass spectrometer for mass and charge analyses.
During each NIPES experiment, the desired (HSO4 -)(VOC) binary clusters were mass-selected and decelerated before being photodetached by a F 2 laser beam (157 nm, 7.867 eV). The laser was operated at a 20 Hz repetition rate with the ion beam off at alternating laser shots to enable shot-by-shot background subtraction. Photoelectrons were collected at nearly 100% efficiency by the magnetic-bottle and analyzed in a 5.2 m long electron flight tube. The TOF NIPE spectra were converted into electron kinetic energy spectra by calibration with the known NIPE spectra of I − and Cu(CN) 2 − . The electron binding energies (EBEs) were obtained by subtracting the electron kinetic energies from the detachment photon energy. The electron energy resolution was about 2%, i.e., 20 meV for 1 eV kinetic energy electrons.
Theoretical methods. Quantum chemical calculations were performed using the NWChem program suite 45 . We have shown previously that the combination of structure optimization using the B3LYP functional followed by single-point energy calculations using the M06-2× functional can give good agreement between calculations and experiments for the atmospheric pre-nucleation clusters 25,29 . Thus we employed the same theoretical approach in this work. The 6-31++ G(d,p) and maug-cc-pVTZ basis sets were used for the structural optimization and singlepoint energy calculation, respectively. Various structures for each cluster have been considered and were optimized without any symmetry constraints. All optimized minima were verified using harmonic vibrational frequency calculations and the calculated frequencies were used to estimate the cluster thermochemistry. Theoretical vertical detachment energies (VDEs) were calculated as the total energy differences between the neutrals and anions both at the optimized anion geometries. All the energies except VDEs reported in this work have been corrected by zero-point vibrational energies (ZPEs) obtained at the B3LYP/6-31++ G(d,p) level of theory. It should be pointed out that our focus in this work is to obtain general views of the structures/interactions formed in between bisulfate and various VOCs, not for obtaining the global minimum for each cluster. Therefore we locate the structures for each cluster by maximizing the hydrogen bonds. Cartesian coordinates are presented in Supplementary Table 8.
Born-Oppenheimer molecular dynamics (BOMD) simulations were performed in a constant volume and temperature (NVT) ensemble using the Amber16 software package 46 . The dispersion and hydrogen bond correction version of semiempirical Hamiltonian PM6, PM6-DH+ 47 , was used to describe the molecular interactions in the BOMD simulations. The (HSO4 -)(VOC) binary clusters with ten water molecules were first optimized at PM6-DH+ level of theory. The systems were then heated up to 298.15 K, followed by an equilibrated run of 20 ps. The produce run was 30 ps in length with a 0.5 fs time step.
Conversion of mass to mobility-equivalent diameter. The mobility-equivalent diameters d p of the sulfuric acid/bisulfate-organic clusters are converted from their masses m p according to the Ehn et al 11,37 .
where d p is the mobility-equivalent diameter, m p is the mass of the cluster ion, ρ is the density of the cluster ion, d g is the diameter of the carrier gas molecule (0.3 nm) 48 , m g is the mass of the carrier gas molecule (28.8 Da).
In the calculation, densities of 1400 kg m −3 were assumed for sulfuric acid/ bisulfate-organic clusters, and 1840 kg m −3 were used for pure sulfuric acid/ bisulfate clusters. The real densities of the organic compounds studied in this work can be found in Supplementary Table 1. It should be pointed out that previous studies found that the sensitivity of d p to the density ρ is relatively small, so the assumption of constant density for these clusters will not induce too much error. The d p of these cluster ions are summarized in Supplementary Table 1.
Details of the evaporation rates calculations. We calculated the evaporation rates of these complexes by employing the method proposed by Ortega et al 49 . The collision rates of ions with neutral molecules were calculated according to the parameterizations from trajectory simulations of collisions between a point charge and a rigidly rotating molecule by Su and Chesnavich 50,51 . In this way, the collision rates β i,j and further the evaporation rates γ i,j are calculated according to the following equations: Among which, β L i;j ¼ q i m À1=2 red ðπα j =ε 0 Þ 1=2 , χ ¼ μ j =ð8πε 0 α j k B TÞ 1=2 and IÃ ¼ μ j I=ðα j qm red Þ; I is the moment of inertia of the neutral molecule. At low values of I* with I Ã <ð0:7 þ χ 2 Þ=ð2 þ 0:6χÞ, the collision rate was noted to be independent of I*. All ion-neutral collisions occurring in this study fall into this low-I* region. q i is the charge of the ion, α j and μ j are the polarizability and dipole moment of the neutral molecule, respectively, ε 0 is the vacuum permittivity. m red is the reduced mass of the collision partners, k B is the Boltzmann constant, and T is the temperature.
ΔG is the Gibbs free energy of formation of the evaporating cluster and the products at temperature T and pressure P ref .
Data availability. All the data generated or analysed during this study are included in this published article (and its electronic supplementary material files). The output files of the quantum chemical calculations associated with this study are available from the authors upon request.