Highly oxygenated organic molecule cluster decomposition in atmospheric pressure interface time-of-flight mass spectrometers

Identification of atmospheric molecular clusters and measurement of their concentrations by atmospheric pressure interface time-of-flight (APi-TOF) mass spectrometers may be affected by systematic error due to possible decomposition of clusters inside the instrument. Here, we perform numerical simulations of decomposition in an APi-TOF mass spectrometers and formation in the atmosphere of a set of clusters which involve a representative kind of highly oxygenated organic molecule (HOM), with the molecular formula C₁₀H₁₆O₈. This elemental composition corresponds to one of the most common mass peaks observed in experiments on ozone-initiated autoxidation of α-pinene. Our results show that decomposition is highly unlikely for the considered clusters, provided their bonding energy is large enough to allow formation in the atmosphere in the first place.

1 Introduction

Recent developments in mass spectrometry have brought huge advancements to the field of atmospheric science. For example, mass spectrometers are now able to detect parts-per-quadrillion-level concentrations of both clusters and precursor vapors in atmospheric samples (Junninen et al., 2010; Jokinen et al., 2012), as well as directly explore the chemistry of new-particle formation (NPF) in the atmosphere (Kulmala et al., 2014; Almeida et al., 2013; Bianchi et al., 2016; Ehn et al., 2014; Thomas et al., 2016; Hogan and de la Mora, 2010).

One of the most common mass spectrometers used to measure online cluster composition and concentration in the atmosphere is the atmospheric pressure interface time-of-flight mass spectrometer (APi-TOF MS). It can detect naturally charged molecules and clusters, or it can be used in combination with a chemical ionization chamber (CI-APi-TOF) to detect neutral molecules and clusters. Unfortunately, the detection process in the APi-TOF MS involves energetic interactions between the carrier gas and the clusters, possibly leading to their decomposition and thus altering the measurement results. Specifically, the ions are guided and focused by an electric field inside the atmospheric pressure interface (APi) through a series of three vacuum chambers before arriving to the time-of-flight mass spectrometer itself. The pressure decreases between successive chambers, until the final value of 10⁻⁶ mbar is reached in the time of flight (TOF). Inside the APi, the clusters can collide with carrier gas molecules and can possibly be decomposed in the process.

Cluster decomposition in the APi-TOF MS is one of the main sources of uncertainty in the measurements. Lack of accuracy for a quantitative estimate of decomposition for example makes it difficult to draw definitive conclusions on the presence (or absence) of certain molecular clusters in the atmosphere. Often the absence of observations of specific clusters by the APi-TOF MS has led to speculation about decomposition inside the mass spectrometer (Olenius et al., 2013). Recently, a numerical model to study decomposition in the APi-TOF MS has been developed by Zapadinsky et al. (2018). The decomposition model has been tested on simple clusters involving one bisulfate anion and two sulfuric acid molecules, giving very good agreement with experimental results (Passananti et al., 2019).

The uncertainties on cluster concentration measurements by the APi-TOF MS and the lack of a comprehensive understanding of decomposition inside the instrument are the motivations of the present work. Here we use a theoretical model to study in detail the decomposition of clusters involving so-called highly oxygenated organic molecules (HOMs), which have recently been identified as a key contributor to NPF (Bianchi et al., 2019). HOMs are molecules formed in the atmosphere from volatile organic compounds (VOCs). Some VOCs with suitable functional groups can undergo an autoxidation process involving peroxy radicals, generating polyfunctional low-volatility vapors (i.e., HOMs) that subsequently condense onto pre-existing particles. HOMs thus contribute to secondary organic aerosol (SOA), which constitutes a significant fraction of the submicron organic aerosol mass and is known to affect the Earth's radiation balance (Jimenez et al., 2009; Donahue et al., 2009; Hallquist et al., 2009). Recent chamber experiments have shown that NPF induced by multicomponent systems, such as sulfuric acid, ammonia and HOMs, could correctly reproduce the NPF events observed in boreal forests (Lehtipalo et al., 2018).

Our study involves a specific kind of representative HOM (C₁₀H₁₆O₈) in the APi. This elemental composition corresponds to one of the most common mass peaks observed in experiments on ozone-initiated autoxidation of α-pinene, which also fulfills the “HOM” definition of Bianchi et al. (2019). The precise molecular structure was adopted from Kurtén et al. (2016) and corresponds to the lowest-volatility structural isomer of the three C₁₀H₁₆O₈ compounds investigated in that study. The structure of the molecule is shown in Fig. 1.

Figure 1 Molecular structure of the model C₁₀H₁₆O₈ HOM used in this study.

The main scope of this work is to determine to what extent we are able to perform measurements of atmospheric cluster concentrations using APi-TOF mass spectrometers. More specifically, we want to determine whether decomposition can possibly be responsible for the lack of observations of some HOM-containing clusters in an APi-TOF. Basically, the formation of clusters in the atmosphere (often referred to as “nucleation”) is driven by the bigger stability of the cluster with respect to its separated molecular (and sometimes ionic) components (fragments). The degree of stability is given by the energy difference between the fragments and the cluster in their (electronic and rotational–vibrational) ground states, which is called either reaction, binding or decomposition energy, and it is equal to the amount of energy necessary to decompose the cluster (or, vice versa, the amount of energy released by the clustering). Higher decomposition energy implies lower decomposition and evaporation rates, and thus higher net formation rates in the atmosphere. On the contrary, decomposition in the APi-TOF MS is enhanced when the decomposition energy decreases, since less energy is required to break the cluster. It is clear at this point that decomposition and new-particle formation (“nucleation”) have opposite dependences on the cluster decomposition energy. For a given temperature and vapor concentration, two situations are possible (Fig. 2):

1. There is a limited range of values for the decomposition energy which allows both decomposition in the APi-TOF and particle formation in the atmosphere.

2. The smallest decomposition energy that allows particle formation in the atmosphere is still too large to permit decomposition in the APi-TOF.

Figure 2 (a) The overlap between the decomposition and particle formation regions defines the range of decomposition energy at which both the events (decomposition in the APi-TOF MS and particle formation in the atmosphere) are allowed. (b) The two events are incompatible.

Here, we predict both an upper bound for decomposition energy necessary for decomposition in the APi-TOF, and a lower bound for new-particle formation in the atmosphere given realistic vapor concentrations. The former can be evaluated using the numerical model developed by Zapadinsky et al. (2018), while the latter is computed using the Atmospheric Cluster Dynamics Code (ACDC) (McGrath et al., 2012). In the present work the simulations on decomposition have been performed by a new C++ version of the decomposition code, which keeps the same basic algorithm but performs faster.

Another purpose of this work is to analyze the dynamics of complex clusters in the APi, showing where and why the decomposition events take place.

In Sect. 2 we present the decomposition model. In Sect. 3 we show the main results on HOM cluster decomposition and we compare them with ACDC simulations. Finally, in Sect. 4, we present the conclusions.

2 Method

The model computes the probability of decomposition of a single negatively charged cluster from a large number (≈10³) of independent stochastic realizations of its dynamics in the APi. The algorithm of the code is presented in the simplified flowchart shown in Fig. 3. We give only a brief description of the algorithm here; for full details see Zapadinsky et al. (2018). The dynamics start with the ionized cluster accelerating in the mass spectrometer under the effect of an electric field generated by the electrodes inside the APi. While moving, a random time interval to the next collision is computed from the cumulative distribution function 𝒫_coll(t), which expresses the probability of encountering a collision after a time t:

$$\begin{array}{}\text{(1)}& {\mathcal{P}}_{\mathrm{coll}}\left(t\right)=\mathrm{1}-e^{-{\int }_{\mathrm{0}}^t\mathrm{{\rm Y}}\left[v\left(t{}^{\prime }\right)\right]\mathrm{d}t{}^{\prime }}.\end{array}$$

Here, Υ is the collision frequency which depends on the cluster velocity v, while t=0 is the moment when the previous collision occurred. Simultaneously, another random time interval is computed for the decomposition event from a Poisson distribution, which corresponds to the time-dependent survival probability 𝒫_surv after collision:

$$\begin{array}{}\text{(2)}& {\mathcal{P}}_{\mathrm{surv}}\left(t\right)=\mathrm{1}-{\mathcal{P}}_{\mathrm{dec}}=e^{-k\left(\mathrm{\Delta }E\right)t},\end{array}$$

where 𝒫_dec is the decomposition probability, and k, the decomposition rate constant, is the inverse of the statistical average of decomposition time, which depends on the cluster excess energy ΔE beyond its decomposition energy threshold. The decomposition time is interpreted as the time the cluster spends intact before decomposing. The decomposition rate constant is derived from phase space theory (PST) of chemical reactions. Basically, it is determined by the ratio of the densities of states of products to the densities of states of reactant, exploiting the detailed balance condition (Zapadinsky et al., 2018).

Figure 3 Simplified flowchart of the decomposition code. The process is repeated up to the required number of realizations.

At this point, if the time required by the cluster to escape from the simulated region of the mass spectrometer is less than both the collision and decomposition times, the single realization is completed and the count of intact clusters is increased by one. A second possibility is that the decomposition time is smaller than collision and escape times, in which case the count of decomposed clusters is increased by 1. Finally, as a third possibility, the collision time may be the smallest, in which case a collision between the cluster and a carrier gas molecule takes place.

The dynamics of the collision are described by a stochastic process, where random velocities for the carrier gas molecules are computed from the Maxwell–Boltzmann distribution. In this case, the motion of the gas molecules is considered as solely translational. During the collision, the kinetic energy of the two colliding objects is partially transferred to the internal degrees of freedom of the cluster (vibrational and rotational modes), according to the microcanonical ensemble approach: any configuration of the system with the same energy is equally probable. The fundamental relation governing these statistics is given by the proportionality between the probability density function of energies for two different interacting degrees of freedom and their densities of states:

$$\begin{array}{}\text{(3)}& f_{a,b}\left({\mathit{ϵ}}_a,{\mathit{ϵ}}_b\right)\propto {\mathit{\rho }}_a\left({\mathit{ϵ}}_a\right){\mathit{\rho }}_b\left({\mathit{ϵ}}_b\right)\mathit{\delta }\left(E-{\mathit{ϵ}}_a-{\mathit{ϵ}}_b\right),\end{array}$$

where ϵ_i is the energy of mode i, ρ_i is the density of states of mode i and the Dirac delta function δ ensures conservation of total energy E. After integrating Eq. (3) over all the possible values that ϵ_b can assume, we get the probability density function for ϵ_a:

$$\begin{array}{}\text{(4)}& f_a\left({\mathit{ϵ}}_a\right)={\displaystyle \frac{{\mathit{\rho }}_a\left({\mathit{ϵ}}_a\right){\mathit{\rho }}_b\left(E-{\mathit{ϵ}}_a\right)}{{\int }_{\mathrm{0}}^E\mathrm{d}{\mathit{ϵ}}_a{\mathit{\rho }}_a\left({\mathit{ϵ}}_a\right){\mathit{\rho }}_b\left(E-{\mathit{ϵ}}_a\right)}}.\end{array}$$

Notice that this formula holds at equilibrium. This assumption is justified from the fact that the energy transfer takes place at timescales much shorter than the timescale between two consecutive collisions. Specifically, vibrational–vibrational energy exchange takes place in about 10⁻¹³–10⁻¹² s, rotational–vibrational exchange in 10⁻¹¹ s and collisions every 10⁻⁹–10⁻⁵ s.

After collision, the energy is then redistributed between rotational and vibrational modes following the same principle, and the dynamics continue with the acceleration of the cluster in the electric field.

In region IV, where the quadrupole is located, the cluster is subjected to additional acceleration in the transverse directions, because of the alternating electric field generated by the quadrupole. The dynamics is determined by the renowned Mathieu equation (Miller and Denton, 1986):

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{d^{\mathrm{2}}y}{d{\mathit{\xi }}^{\mathrm{2}}}}+\left(a-\mathrm{2}q\cos \mathrm{2}\mathit{\xi }\right)y=\mathrm{0},\end{array}$$

where

$$\begin{array}{}\text{(6)}& a={\displaystyle \frac{\mathrm{4}eU}{m{r}_{\mathrm{0}}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}}},q={\displaystyle \frac{\mathrm{2}eV}{m{r}_{\mathrm{0}}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{2}}}}\text{and}\mathit{\xi }={\displaystyle \frac{\mathit{\omega }t}{\mathrm{2}}},\end{array}$$

with y one of the two transverse displacements from the quadrupole axis, t the time, e the elementary charge, U the DC component, V the AC amplitude of the electric potential, m the cluster mass, r₀ the half-distance between quadrupole rods and ω the AC angular frequency.

The computation of the density of vibrational and rotational states requires the knowledge of the corresponding energy levels, which are computed by the quantum chemistry program Gaussian, using the PM7 semi-empirical method, which is the newest and generally best semi-empirical method available in Gaussian (Frisch et al., 2016). We note that while semi-empirical methods are unable to accurately model the energetics of molecular clustering, this is not a problem in the present study, as the binding/decomposition energy is treated here as a freely variable parameter (i.e., the unreliable PM7 binding energy is not actually used). The purpose of the PM7 optimizations and frequency calculations is simply to provide a qualitatively correct distribution of rotational and vibrational energy levels. As vibrational modes with wavenumbers above 2500 cm⁻¹ are strongly underestimated by the PM7 method, we have further rescaled the frequencies following the relation suggested in Rozanska et al. (2014):

$$\begin{array}{}\text{(7)}& \begin{array}{ll}{\mathit{\nu }}_{\mathrm{PM}\mathrm{7}}^{\mathrm{scaled}}& ={\mathit{\nu }}_{\mathrm{PM}\mathrm{7}}^{\mathrm{unscaled}}\\ & \text{for}\mathit{\nu }<\mathrm{125}{\text{ cm}}^{-\mathrm{1}},\\ & \\ {\mathit{\nu }}_{\mathrm{PM}\mathrm{7}}^{\mathrm{scaled}}& =\mathrm{1.006}{\mathit{\nu }}_{\mathrm{PM}\mathrm{7}}^{\mathrm{unscaled}}-\mathrm{27.7}\\ & \text{for}\mathrm{125}{\mathrm{cm}}^{-\mathrm{1}}<\mathit{\nu }<\mathrm{2500}{\mathrm{cm}}^{-\mathrm{1}},\\ & \\ {\mathit{\nu }}_{\mathrm{PM}\mathrm{7}}^{\mathrm{scaled}}& =\mathrm{3.99}{\mathit{\nu }}_{\mathrm{PM}\mathrm{7}}^{\mathrm{unscaled}}-\mathrm{7849.2}\\ & \text{for}\mathit{\nu }>\mathrm{2500}{\mathrm{cm}}^{-\mathrm{1}}.\end{array}\end{array}$$

The vibrational density of states is then computed using the harmonic approximation, which is most likely the biggest source of error in the evaluation of the cluster survival probability. (See Sect. 3 for a sensitivity analysis on the effect of varying the vibrational frequencies.)

3 Results

3.1 Decomposition in the APi

The simulation of decomposition involves only a portion of the total length of the APi-TOF mass spectrometer. Specifically, the simulations take place in the most critical region, ranging from the end of the first to the end of the second chamber of the APi, where the pressure values make decomposition possible: at the beginning of the first chamber the collisions are not energetic enough, while in the third chamber the carrier gas is so sparse that no collision happens (Zapadinsky et al., 2018).

In Fig. 4 we have sketched the simulated section of the instrument. We subdivide the simulated section into five regions with different longitudinal electric field and pressure values. Region I is at the end of the first chamber; region II defines the interface between the first and the second chambers; and regions III, IV and V are located in the second chamber of the APi. The pressure in the first region is P₁, while in regions III, IV and V it is P₂. In region II, where the skimmer is located, the pressure changes continuously from P₁ to P₂. The electric field takes different values in regions I, III, IV and V, while in II it is set to zero. The voltage and pressure configuration inside the APi chambers, used in the CLOUD10 experiments (Lehtipalo et al., 2018), are shown in Table 1, where V_DC and V_AC are the direct and alternating components of the quadrupole electric potential in region IV, while Ω is its radio frequency. We will adopt the same parameter configuration in our simulations. It is important to notice here that the electric field takes very different values in different regions, varying over 3 orders of magnitude, which will greatly diversify the probability of decomposition at different locations.

Figure 4 Simulated region of APi-TOF mass spectrometer. The dashed lines identify the regions with different voltages. The thin solid line illustrates a single realization of a cluster trajectory. The regions' lengths are not to scale.

Table 1 Parameters configuration of APi used in the CLOUD10 experiments and adopted in the simulations.

In this study the cluster decomposition energy E_f is treated as a free parameter. This allows us to explore the behavior of cluster survival probability as a function of E_f over a large energy range. Moreover, keeping E_f as a variable is also useful since its computation by quantum chemistry calculations (especially low-level methods such as PM7 used here) is affected by large errors.

The clusters studied here are formed by one bisulfate anion, one to two sulfuric acid molecules and one to three HOM molecules. The clusters are assumed to decompose by losing one HOM, as follows:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& \left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right){\left({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}\right)}_n{\left({\mathrm{HOM}}_{\mathrm{10}}\right)}_m\\ & \to \left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right){\left({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}\right)}_n{\left({\mathrm{HOM}}_{\mathrm{10}}\right)}_{m-\mathrm{1}}+\left({\mathrm{HOM}}_{\mathrm{10}}\right),\end{array}\end{array}$$

with $n=\mathrm{1},\mathrm{2}$, $m=\mathrm{1},\mathrm{2},\mathrm{3}$ and HOM₁₀ is the structure shown in Fig. 1. Our specific choice for HOM is not unique: this is one among many potential HOM produced in the α-pinene +O₃ reaction. The compound chosen is broadly representative of autoxidation products as it contains both hydroperoxide, ketone and carboxylic acid groups. The clusters have been constructed by first maximizing H-bonds between the ${\mathrm{HSO}}_{\mathrm{4}}^{-}$ core ion and other molecules, and then maximizing other H-bonds without creating too much strain. We note that the cluster conformers obtained in this fashion are unlikely to correspond precisely to the global energy (or free energy) minima. However, this mainly affects the computed binding energy – which (as pointed out earlier) is not actually used in this study. For the purposes of generating a representative ensemble of vibrational and rotational energy levels, the conformer generation approach used here is adequate.

For each kind of cluster we performed several simulations at different E_f. The simulated carrier gas is air, modeled to consist of 80 % nitrogen and 20 % oxygen.

The final survival probabilities P_surv are shown in Fig. 5. As expected, P_surv is monotonic with respect to E_f. From these results we can identify a boundary at E_f≈25 kcal mol⁻¹ beyond which decomposition is highly unlikely. The weakest bonds between non-hydrogen atoms in HOM molecules are likely to be the O−O bonds of peroxide or hydroperoxide groups. These typically have bond dissociation energies around 35–50 kcal mol⁻¹ (Bach et al., 1996; Schweitzer-Chaput et al., 2015, 2017). Vinyl peroxide systems have much lower dissociation energies, but (with the exception of short-lived vinyl hydroperoxides generated in ozonolysis) these are unlikely to form in gas-phase oxidation processes. The threshold cluster decomposition energies reported here are thus much lower than the energies required to dissociate even the weakest covalent bonds in the HOM molecules.

Figure 5 Survival probabilities P_surv for different HOM clusters as a function of their decomposition energy E_f. The error bars are relative only to statistical error.

Figure 6 (a) Distribution of vibrational energy at equilibrium at temperature T=300 K. (b) Decomposition rate constants at E_f=20 kcal mol⁻¹. Collision frequencies in the first and second chambers are shown as reference.

Larger clusters are found to decompose more easily. This can be understood looking at the probability density function (PDF) of vibrational energy at equilibrium at temperature T=300 K (Fig. 6a). The distributions have been computed multiplying the density of states of vibrational modes by the Boltzmann factor, and then renormalized:

$$\begin{array}{}\text{(9)}& \text{PDF}\left(E\right)={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{v}}\left(E\right)e^{-E/k_{\mathrm{B}}T}}{{\int }_{\mathrm{0}}^{+\mathrm{\infty }}{\mathit{\rho }}_{\mathrm{v}}\left(E\right)e^{-E/k_{\mathrm{B}}T}}}.\end{array}$$

Larger clusters contain a higher number of bonds, which, for fixed decomposition energy, will each increase the internal energy of the cluster, because of the equipartition theorem. This leads to a higher probability of exceeding the decomposition energy, even at equilibrium conditions, i.e., in the absence of an electric field. Suppose, for example, E_f = 20 kcal mol⁻¹. Figure 6a shows that the smallest cluster has roughly 50 % probability of having an internal energy higher than E_f, while the biggest cluster is certainly at higher energy than E_f, with an average near 50 kcal mol⁻¹. This effect is in part counterbalanced by the behavior of the decomposition rate constant as a function of cluster size. As we can see in Fig. 6b, the decomposition rate at fixed internal energy tends to be smaller for larger clusters. Indeed, a larger number of bonds reduces the probability of finding a large quantity of energy in a single bond.

The model requires as input parameters the vibrational and rotational frequencies of both the reactant and products from an external quantum chemistry program (in our case Gaussian 16). These input data are affected by errors that depend on the level of accuracy of the computational method. It is therefore interesting to study how errors in the inputs affect the survival probability. For this reason, we generated new frequencies, adding random noise to the original ones. The new frequencies are generated from normal distributions centered at the original frequencies, with a standard deviation equal to 20 % of their values. Figure 7 shows the deviations of the survival probability for 10 different sets of random input frequencies. Here we can see that, even with large errors in the input parameters (20 %), the final results change by no more than 10 %. These results demonstrate the low sensitivity of the model on deviations in the input frequencies, thus validating our use of rather crude approaches for both the quantum chemical calculations, and the conformer generation.

Figure 7 Deviations (thin lines) from original survival probability P_surv (thick line) for 10 different sets of random input frequencies (noisy frequencies), obtained adding random noise to the original frequencies. The new sets of frequencies have been drawn from normal distributions having the original frequency as mean and with a standard deviation equal to 20 % of their mean. The results are relative to the $\left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right)\left({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}\right)\left({\mathrm{HOM}}_{\mathrm{10}}\right)$ cluster.

Let us now analyze the locations where collisions take place in the APi. In particular, we are interested in the highest-energy collisions, which lead to decomposition. We dub them as fatal collisions. For this study we used the voltage configuration and pressures used in CLOUD10 and we collected data from 5000 realizations of $\left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right)\left({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}\right)\left({\mathrm{HOM}}_{\mathrm{10}}\right)$ cluster dynamics at E_f=17.88 kcal mol⁻¹, corresponding to a survival probability P_surv≈50 % (see Fig. 5).

Figure 8 (a) Average number of total and fatal collisions per realization of cluster evolution in the different APi regions. (b) Ratio between total and fatal collisions. The results have been obtained from 5000 realizations of $\left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right)\left({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}\right)\left({\mathrm{HOM}}_{\mathrm{10}}\right)$ cluster dynamics. Notice that regions' lengths of the APi are different; in particular region IV, where the quadrupole is located, is much longer than the other regions, which leads to a higher number of collisions.

Figure 8a shows the average number of total and fatal collisions per single cluster evolution in different regions of the APi. The average number of fatal collisions is computed by dividing the total number of fatal collisions by the number of realizations (dark bars in Fig. 8a), and it can be interpreted as the probability of the cluster to decompose in the correspondent region in a complete realization. Summing up these values of all five regions gives the total probability of decomposition. Figure 8b instead shows the ratio between fatal and total collisions, which can be interpreted as the probability of the cluster to decompose after experiencing a collision in the correspondent region. Therefore, Fig. 8b identifies the regions where the collisions are most energetic (regions III and V).

The cluster dynamics starts in the first region (end of first APi chamber), where the high pressure (201 Pa) generates a large number of collisions (≈200 collisions per realization). Here the ratio between fatal collisions and total collisions is about 0.07 % (Fig. 8b), a sign of low-energy collisions despite the very strong electric field (7189 V m⁻¹). When the cluster moves to region II, in the skimmer, we see a drop in collisions, reflecting the negative pressure gradient. Moreover, the absence of electric field in this region makes fatal collisions very unlikely to happen (≈0.017 % probability). In region III, at the beginning of the second APi chamber, the pressure is low (2.96 Pa), and the cluster experiences few collisions. However, the strong electric field (804.9 V m⁻¹) leads to cluster acceleration and an increase in fatal collisions, about 0.4 % of total collisions. The passage to region IV is characterized by a high speed of the cluster acquired in the previous region. This leads to a large number of decompositions, as shown in Fig. 8a. Subsequently, the cluster decelerates because of the very weak electric field in this region (18.04 V m⁻¹). Since the cluster moves for a long distance at low velocities in region IV, the average probability of decomposition per collision drops to ≈0.05 %. Finally, in region V, the cluster accelerates again due to a strong electric field (1104 V m⁻¹), which is visible from the large increase in the probability of fatal collisions (≈2 %). The cumulative survival probability of clusters as they travel through regions I to V is illustrated by the red dots in Fig. 9. The decomposition probabilities in each respective region are indicated by the histogram bars.

Figure 9 Cumulative survival probability after going through each region (red dots) and decomposition probability (bars) for a single cluster realization. The decomposition probability in each region is given by the ratio of clusters decomposing inside the region, over clusters entering the region. The results are relative to the $\left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right)\left({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}\right)\left({\mathrm{HOM}}_{\mathrm{10}}\right)$ cluster.

Although the decomposition is mainly induced by the voltages applied between regions II and III as shown in Passananti et al. (2019), these results demonstrate that the decomposition takes place mainly in region IV (Fig. 8a). This is due to the large acceleration acquired by the cluster in region III, with subsequent possible collision at high energy in the beginning of the next region. These results are complementary to previous data reported in the literature (Passananti et al., 2019; Lopez-Hilfiker et al., 2016) and help to define the terms “declustering voltages” as the voltages that induce decomposition (between regions II and III) and the “declustering region” as the region of the APi where the decomposition takes place (region IV).

Figure 10 (a) Comparison between survival probabilities in presence (V_AC=200 V) and absence (V_AC=0 V) of the quadrupole electric field. (b) Probability of decomposition in the quadrupole region, with and without quadrupole electric field. The results are relative to the $\left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right)\left({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}\right)\left({\mathrm{HOM}}_{\mathrm{10}}\right)$ cluster.

In region IV a quadrupole is used as ion guide in order to focus the cluster trajectories, through an alternating transverse electric field (Miller and Denton, 1986). For this reason, the cluster is also subjected to a transverse acceleration in addition to the longitudinal one, increasing the velocity and, consequently, the probability of decomposition. In order to separate the contributions to decomposition from the longitudinal and the transverse electric fields, we performed additional simulations in which the quadrupole electric potential was set to V_AC=0. In this situation the motion in the transverse directions is negligible (as in the other regions of the APi). The increase in the final survival probabilities can be seen in Fig. 10a, which shows a minor deviation from the case with V_AC=200 V. More specifically, Fig. 10b shows the slight reduction in decomposition probability within the quadrupole region when V_AC=0. Hence, with the present voltage configuration, the quadrupole voltage does not affect the probability of decomposition significantly.

3.2 Formation in the atmosphere

In parallel to the simulations on decomposition in the APi, we used the ACDC code (McGrath et al., 2012) to compute the concentrations of HOM₁₀ needed to provide a significant enhancement of new-particle formation rate in the atmosphere (details on ACDC simulations provided in Appendix A). Specifically, at given H₂SO₄ monomer concentration and HOM cluster decomposition energy E_f, we computed the concentration of HOM₁₀ at which the new-particle formation rate J is increased by 1 cm⁻³ s⁻¹ with respect to the pure sulfuric acid system (i.e., $J_{\mathrm{sa}+\mathrm{hom}}=J_{\mathrm{sa}}+\mathrm{\Delta }J$, with ΔJ=1 cm⁻³ s⁻¹).

Figure 11a is useful to understand the framework of the simulations: the scheme shows the clusters involved in the ACDC simulations: along the x axis we have clusters with an increasing number of sulfuric acid molecules, while along the y axis we have clusters with an increasing number of HOM₁₀. Clusters grow through collisions with bisulfate anions, sulfuric acid and HOM₁₀ molecules, until they reach a size bigger than three HOM₁₀ or three sulfuric acid (SA) molecules. At this point the clusters exit the simulation region and they are accounted as outgrowing (stable) clusters. ACDC code computes the rate of formation and evaporation of all the clusters in the simulation region. We are interested in the outgoing fluxes (blue and red arrows in Fig. 11a), which are the ones that contribute to new-particle formation. In particular, we are interested in the enhancement of the total outgoing flux (formation rate) when HOM₁₀ molecules are added to a system of sulfuric acid particles. The increment of the formation rate (ΔJ) is represented by the red arrows, and we want to find the concentration of HOM₁₀ that corresponds to ΔJ=1 cm⁻³ s⁻¹. This value is indicative, and it serves as a reference for a reasonable NPF process in the atmosphere (Lehtipalo et al., 2018). The concentration of ${\mathrm{HSO}}_{\mathrm{4}}^{-}$ is kept constant at 700 cm⁻³, which corresponds to the steady-state concentration at the ion formation rate J_ions=4 cm⁻³ s⁻¹, reproducing the galactic cosmic ray (GCR) conditions (Kirkby et al., 2016).

Figure 11 (a) Clusters involved in the ACDC simulations. Out of the box, the clusters are considered to be stable and contributing to new-particle formation. The blue arrow indicates the flux of the outgrowing cluster that contributes to the formation rate of only sulfuric acid particles (J_sa). After adding HOM₁₀ monomers to the system, new outgrowing HOM clusters (red arrows) are participating in new-particle formation, resulting in a total formation rate J_sa+hom. (b) Concentrations of HOM₁₀ needed to increment the new-particle formation rate by 1 cm⁻³ s⁻¹, for different decomposition energies.

Different simulations have been performed for different sulfuric acid monomer concentrations and different decomposition energies E_f. Figure 11b shows the dependence of HOM₁₀ concentration on the SA monomer concentration, in order to maintain ΔJ=1 cm⁻³ s⁻¹, for different decomposition energies.

As we can see, the H₂SO₄ and HOM₁₀ concentrations decrease when the decomposition energy increases, reaching experimental conditions (e.g., non-nitrate HOM₁₀ and H₂SO₄ concentrations measured in the CLOUD experiment; Lehtipalo et al., 2018) at E_f>30 kcal mol⁻¹, where HOM₁₀ and H₂SO₄ concentrations do not exceed 10⁸ cm⁻³. At this decomposition energy, decomposition in APi is negligible (P_surv≈1 in Fig. 5). Thus we can conclude that for the case of SA–HOM clusters, rapid formation in the atmosphere (given typical vapor concentrations) and significant decomposition in the APi are mutually incompatible situations (case b in Fig. 2).

4 Conclusions

In this work, we have presented the numerical results on decomposition inside an APi-TOF MS instrument of a specific class of atmospheric clusters that involve sulfuric acid and HOM molecules. A previously reported low-volatility α-pinene ozonolysis product, with the molecular formula C₁₀H₁₆O₈, was used as a representative HOM.

There are three main results from our simulations. First, decomposition of HOM clusters in the APi requires a range of cluster decomposition energies which is incompatible with efficient cluster formation in the atmosphere given sub-parts-per-billion vapor concentrations. This result has been obtained by computing the cluster survival probability in the APi as a function of its decomposition energy, and then comparing the range of energies leading to decomposition with those needed to obtain a reasonably high new-particle formation rate in atmospheric conditions. The two ranges have no overlap – the highest energy allowing decomposition differs from the lowest energy allowing atmospheric new-particle formation by roughly 10 kcal mol⁻¹. Observations of SA–HOM clusters in CLOUD experiments (Lehtipalo et al., 2018) validate the results of our model.

Our second main result is the identification of the locations of the highest-energy collisions that lead to decomposition of SA–HOM clusters in the APi. The simulation shows that they are mainly localized in the quadrupole region (IV). This is due to the large velocity of the cluster acquired in the previous region (III), caused by the strong electric field. Moreover, our results show that the alternating field of the quadrupole in region IV has only a minor effect on decomposition probability.

As a third main result, we have shown that the model displays low sensitivity to changes in cluster vibrational and rotational frequencies. This feature allows us to also obtain reliable results with low-level input data (e.g., semi-empirical methods and limited configurational sampling).

This study was performed for a small set of HOM clusters with particular functional groups (e.g., hydroperoxides, ketones), but in future similar simulations could be performed for other types of clusters, for example to assess whether they should be detectable in an APi-TOF MS or not, or whether they could possibly affect new-particle formation, despite not being directly detected.

Appendix A: ACDC simulation details

The configurational sampling of the H₂SO₄ molecule and molecular clusters $\left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right)({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}{)}_{\mathrm{0}-\mathrm{3}}$ has been performed as described by Kubečka et al. (2019) to obtain their thermodynamic (T=300 K) properties at a high level of theory (DLPNO-CCSD(T)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) (Kubečka et al., 2019; Myllys et al., 2016; Riplinger and Neese, 2013; Riplinger et al., 2013, 2016; Chai and Head-Gordon, 2008; Dunning, 1989; Kendall et al., 1992; Liakos et al., 2015)).

The formation Gibbs free energies were used as input for the Atmospheric Cluster Dynamic Code (McGrath et al., 2012), which uses them to obtain cluster stabilities (i.e., evaporation rates). Bisulfate ions ${\mathrm{HSO}}_{\mathrm{4}}^{-}$ in the atmosphere are formed by collisions of “generic” anions and sulfuric acid H₂SO₄ (i.e., as H₂SO₄ is practically the strongest acid in the atmosphere, collision of H₂SO₄ with almost any anion will tend to produce ${\mathrm{HSO}}_{\mathrm{4}}^{-}$).

In the simulation, we assume a constant concentration of bisulfate ions equal to 700 cm⁻³, which approximately corresponds to the steady-state anion concentration at a representative total ion pair formation rate due to galactic cosmic rays (GCRs) of J_ions=4 cm⁻³ s⁻¹ (Kirkby et al., 2016). We assume that formation of $\left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right)({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}{)}_{\mathrm{4}}$ (outgrowing cluster) or larger clusters leads to irreversible new-particle formation. Thus, at a specific H₂SO₄ monomer concentration, we simulated new-particle formation of sulfuric acid aerosols until a steady state had been reached. For each monomer concentration, the new-particle formation rate J_sa were noted down. (For simplicity, we neglect losses due to ion–ion recombination, wall loss, dilution or coagulation, as our focus is on comparing the relative effect of HOM.)

As described in Kubečka et al. (2019), we performed “weak” configurational sampling of molecules HOM₁₀, H₂SO₄ and molecular clusters $\left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right)\left({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}{)}_{\mathrm{0}-\mathrm{3}}\right({\mathrm{HOM}}_{\mathrm{10}}{)}_{\mathrm{0}-\mathrm{3}}$ to obtain representative thermodynamic properties of these clusters at a low level of theory (PM7, Stewart, 2013). By “weak” configurational sampling we mean that just one conformation of HOM₁₀ is included for generating of HOM-containing clusters. However, as only the rotational and vibrational energy levels from the PM7 calculations are actually used, and not the electronic energies, this approach is sufficient for our purposes. We further assume that the evaporation rate of H₂SO₄ from $\left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right)\left({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}{)}_x\right({\mathrm{HOM}}_{\mathrm{10}}{)}_{\mathrm{1}-\mathrm{3}}$ clusters is identical to the evaporation rate of H₂SO₄ from the previously calculated corresponding $\left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right)({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}{)}_x$ clusters. The evaporation rate γ of HOM₁₀ from $\left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right)\left({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}{)}_x\right({\mathrm{HOM}}_{\mathrm{10}}{)}_y$ is calculated as

$$\begin{array}{}\text{(A1)}& \mathit{\gamma }=\mathit{\beta }\cdot {\displaystyle \frac{p_{\mathrm{ref}}}{k_{\mathrm{B}}T}}\cdot e^{\mathrm{\Delta }G/k_{\mathrm{B}}T},\end{array}$$

where k_B is Boltzmann's constant, T is temperature, β is the collision rate of HOM₁₀ and $\left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right)\left({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}{)}_x\right({\mathrm{HOM}}_{\mathrm{10}}{)}_{y-\mathrm{1}}$, p_ref is the reference pressure, and ΔG is the free energy change of the reaction

$$\begin{array}{rl}& {\mathrm{HOM}}_{\mathrm{10}}+\left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right){\left({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}\right)}_x{\left({\mathrm{HOM}}_{\mathrm{10}}\right)}_{y-\mathrm{1}}\\ & ⟶\left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right){\left({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}\right)}_x{\left({\mathrm{HOM}}_{\mathrm{10}}\right)}_y.\end{array}$$

ΔG is calculated as

$$\begin{array}{}\text{(A2)}& \mathrm{\Delta }G=\mathrm{\Delta }{G}^{\mathrm{PM}\mathrm{7}}-\mathrm{\Delta }{E}_{\mathrm{el}}^{\mathrm{PM}\mathrm{7}}+E_{\mathrm{f}}=\mathrm{\Delta }{G}_{\mathrm{therm}\mathrm{\&}\mathrm{ZPE}}^{\mathrm{PM}\mathrm{7}}+E_{\mathrm{f}},\end{array}$$

where E_f is the parametrized value of the decomposition energy (varying between 15 and 30 kcal mol⁻¹ as shown), ΔG^PM7 and $\mathrm{\Delta }{E}_{\mathrm{el}}^{\mathrm{PM}\mathrm{7}}$ are the reaction Gibbs free energy and electronic energy (respectively) computed at the PM7 level of theory, and $\mathrm{\Delta }{G}_{\mathrm{therm}\mathrm{\&}\mathrm{ZPE}}^{\mathrm{PM}\mathrm{7}}$ corresponds to the thermal (and zero-point vibrational) correction to the reaction Gibbs free energy at the PM7 level.

Again, in the same way as in the decomposition model, we used the evaporation rates as input for the ACDC program assuming $\left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right)\left({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}{)}_{\mathrm{4}}\right({\mathrm{HOM}}_{\mathrm{10}}{)}_{\mathrm{0}-\mathrm{3}}$ or $\left({\mathrm{HSO}}_{\mathrm{4}}^{-}\right)\left({\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}{)}_{\mathrm{0}-\mathrm{3}}\right({\mathrm{HOM}}_{\mathrm{10}}{)}_{\mathrm{4}}$ to be the outgrowing clusters. Next, at given H₂SO₄ monomer concentration, and for a given value of E_f, we searched for the HOM₁₀ concentration at which the new-particle formation rate is increased by 1 cm⁻³ s⁻¹ (i.e., $J_{\mathrm{sa}+\mathrm{hom}}=J_{\mathrm{sa}}+\mathrm{1}$ cm⁻³ s⁻¹). The resulting graph is shown in the main article in Fig. 11b.

The above computations involved the following programs: the ABCluster program (Zhang and Dolg, 2015, 2016), the GFN2-xTB program (Grimme et al., 2017; Bannwarth et al., 2019), the Gaussian 16 program Revision A.03 (Frisch et al., 2016), the GoodVibes program (Funes-Ardois and Paton, 2016) and the Orca program version 4.0.1.2 (Neese, 2012).