A cluster-of-functional-groups approach for studying organic enhanced atmospheric cluster formation

The role of organic compounds in atmospheric new particle formation is difficult to disentangle due to the myriad of potentially important oxygenated organic molecules (OOMs) present in the atmosphere. Using state-of-the-art quantum chemical methods, we here employ a novel approach, denoted the “cluster-of-functional-groups” approach, for studying the involvement of OOMs in atmospheric cluster formation. Instead of the usual “trial-and-error” approach of testing the ability of experimentally identified OOMs to form stable clusters with other nucleation precursors, we here study which, and how many, intermolecular interactions are required in a given OOM to form stable clusters. In this manner we can reverse engineer the elusive structure of OOM candidates that might be involved in organic enhanced atmospheric cluster formation.

We calculated the binding free energies of all combinations of donor and acceptor organic functional groups to investigate which functional groups most preferentially bind with each other and with other nucleation precursors such as sulfuric acid and bases (ammonia, methyl-, dimethyl- and trimethylamine). We find that multiple carboxyl groups lead to substantially more stable clusters compared to all other combinations of functional groups. Employing cluster dynamics simulations, we investigate how a hypothetically OOM composed of multiple carboxyl groups can stabilize sulfuric acid–base clusters and provide recommendations for potential atmospheric multi-carboxylic acid tracer compounds that should be explicitly studied in the future.

The presented “cluster-of-functional-groups” approach is generally applicable and can be employed in many other applications, such as ion-induced nucleation and potentially in elucidating the structural patterns in molecules that facilitate ice nucleation.

1 Introduction

Secondary organic aerosols (SOAs) constitute a major fraction of organic matter in the atmosphere (Jimenez et al., 2009) and influence both human health (Pelucchi et al., 2009) and global climate (Canadell et al., 2021). The term SOA formation usually comprises both the initial aerosol nucleation and the subsequent growth of existing aerosol particles via vapor uptake. However, while oxygenated organic molecules (OOMs) are widely accepted as being important for aerosol growth, it remain ambiguous whether OOMs are important for aerosol nucleation in the atmosphere (Kirkby et al., 2023).

The onset of aerosol nucleation is governed by the formation of stable atmospheric molecular clusters (Kulmala et al., 2013). Sulfuric acid (Sipilä et al., 2010) and atmospheric bases such as ammonia (AM) (Kirkby et al., 2011) and amines (Almeida et al., 2013; Jen et al., 2014; Glasoe et al., 2015) (methylamine (MA), dimethylamine (DMA) and trimethylamine (TMA)) are important species for nucleation both over land and sea. In addition, iodine species have been shown to contribute to nucleation at coastal and marine environments (Baccarini et al., 2020; Beck et al., 2021; He et al., 2021, 2023). The global modeling study by Dunne et al. (2016) indicated that almost all nucleation occurring in the atmosphere involves sulfuric acid coupled with either ammonia and/or organics. However, the exact role of organics in aerosol nucleation remains elusive. Pure organic ion-induced nucleation has been reported both at the CLOUD chamber (Kirkby et al., 2016) and in the field (Rose et al., 2018). While some early work has inferred that organics play a role in multicomponent cluster formation (Zhang et al., 2004; Metzger et al., 2010; Schobesberger et al., 2013; Riccobono et al., 2014), it remains unknown whether this is also true for multicomponent SA–base–OOM nucleation or whether it proceeds via two isolated pathways (one for SA–base and one for pure OOMs). Both experiments (Lehtipalo et al., 2018; Kirkby et al., 2023) and quantum chemical calculations (Elm, 2019b) allude to a decoupled mechanism, but this has not yet been definitively confirmed.

The puzzle of the role of organics in aerosol nucleation originates from the fact that there exists a myriad of OOMs in the atmosphere. Volatile organic compounds (VOCs) are rapidly transformed into less volatile species in the air due to oxidation reactions with either OH, NO₃, or O₃, making it challenging to disentangle which OOMs are important for nucleation and which are important for SOA growth. Identification of OOMs is usually performed with mass spectrometer techniques such as chemical ionization atmospheric pressure interface time of flight (CI-APi-TOF) (Jokinen et al., 2012). For instance, the seminal work by Ehn et al. (2014) measured the chemical composition of OOM monomers (C₁₀H_14-16O_7-11) and OOM covalently bound dimers (C_19-20H_28-32O_10-18) from α-pinene (C₁₀H₁₆) oxidation. While this gives important information about the chemical composition of the OOMs, it yields no information about the exact functional groups or their arrangement in the molecules. Quantum chemical (QC) calculations can be applied to yield the cluster structures and in conjugation with cluster dynamics modeling can elucidate the role of specific OOMs in nucleation. However, this requires that the given OOM structure is known. This has led to an overabundance of QC studies, including work by our group, that investigate the binding affinity of different experimentally identified OOMs to themselves and other nucleation precursors. We recently reviewed the entire QC literature on the role of organics in cluster formation, and unfortunately, not a single OOM has definitively been proven to participate in nucleation in the planetary boundary layer (Elm et al., 2023). The lack of progress could be ascribed to the fact that previous work has been looking at the wrong compounds. All studies have been performed on the organic monomers, while recent evidence from the CLOUD chamber has shown that it is in fact the covalently bound organic dimers, which have ultra low volatilities (ULVOCs) (Simon et al., 2020), that drive nucleation (Lehtipalo et al., 2018; Dada et al., 2023). We will hereon refer to these dimers as accretion products to clearly distinct them from dimer clusters. Studying large C_19-20H_28-32O_10-18 OOM accretion products such as those from α-pinene oxidation is challenging using QC methods, and to date not a single study exists. This is caused by the fact that QC methods scale steeply with system size. In addition, larger, more flexible molecules can exist in numerous conformations, which rapidly increase number of calculations required to identify the global minimum cluster structures.

Using state-of-the-art quantum chemical methods, we here employ a new approach, denoted the “cluster-of-functional-groups” approach (Elm et al., 2023), for studying the involvement of OOMs in atmospheric cluster formation. Instead of the usual “trial-and-error” approach of testing the ability of identified OOMs to form stable clusters with themselves or other nucleation precursors, we here study exactly which and how many intermolecular interactions (in the form of functional groups) are required in a given OOM to form stable clusters. This allows us to reverse engineer the potential structure of OOMs involved in organic enhanced atmospheric cluster formation. We explicitly study all possible donors (alcohol, peroxide) and acceptors (ether, epoxide, aldehyde, ketone, acid anhydride and ester), as well as carboxylic acids. Based on the functional groups that bind strongest, we extend the analysis to study the cluster dynamics of a hypothetical OOM binding to $\left(\text{SA}{)}_{\mathrm{1}-\mathrm{2}}\right(\text{base}{)}_{\mathrm{1}-\mathrm{2}}$ clusters, where the bases are AM, MA, DMA and TMA.

2 Methods

2.1 Computational details

The semi-empirical GFN1-xTB (Grimme et al., 2017) energy calculations and geometry optimization were calculated with the xtb 6.4.0 program (Bannwarth et al., 2021). Gaussian16, version B.01 (Gau, 2016), was used for the density functional theory (DFT) calculations with default convergence criteria. The ωB97X-D functional (Chai and Martin, 2008) with the 6-31++G(d,p) basis set was chosen for geometry optimization and vibrational frequency calculations based on its performance in numerous benchmarks (Elm and Mikkelsen, 2014; Myllys et al., 2016; Elm and Kristensen, 2017; Schmitz and Elm, 2020; Jensen et al., 2022). Given that some vibrational frequencies were low, the quasi-harmonic approximation, as recommended by Grimme (2012), was used to treat vibrational frequencies below 100 cm⁻¹ as free rotors. ORCA 5.0.4 (Neese, 2012, 2018; Neese et al., 2020; Neese, 2022) was used to calculate single-point energy corrections using Domain Local Pair Natural Orbital, DLPNO–CCSD(T₀) (Riplinger and Neese, 2013; Riplinger et al., 2013), with the aug-cc-pVTZ basis set, using the TightSCF convergence criteria and the NormalPNO setting (Liakos et al., 2015). The workflow and subsequent data processing was automated using the JKCS program (Kubečka et al., 2023). All data have been added to the Atmospheric Cluster DataBase (ACDB) (Elm, 2019a). The atmospheric cluster dynamics were simulated with the Atmospheric Cluster Dynamics Code (ACDC) (McGrath et al., 2012) using a modified version of the code given in the online repository by Tinja Olenius (Olenius et al., 2013; Roldin et al., 2019) available in the JKCS suite. The simulations were done at 278.15 K with a constant coagulation sink of $-\mathrm{1.6}\times {\mathrm{10}}^{-\mathrm{3}}$ s⁻¹, matching typical sink values (Dal Maso et al., 2008; Kontkanen et al., 2017). The clusters with a single sulfuric acid monomer were counted towards the total sulfuric acid concentration.

2.2 Configurational sampling

For the configurational sampling procedure, a funnel-type approach was employed (Temelso et al., 2018; Odbadrakh et al., 2020; Kubečka et al., 2019), where the level of theory is increased in each step as the number of candidate cluster structures is decreased. The workflow can be described as follows:

$$\begin{array}{rl}\text{ABCluster}& \to {\text{GFN1-xTB}}^{\text{OPT}}\to \text{ArbAlign}\\ & \to \mathit{\omega }{\text{B97X-D}}^{\text{OPT+FREQ}}\to {\text{DLPNO}}^{\text{SP}}.\end{array}$$

ABCluster (Zhang and Dolg, 2015, 2016) was used to generate a large pool of cluster structures with the CHARMM force field. Settings as recommended by Kubečka et al. (2019) with population size, SN=3000, maximum generations, g_max=200, and number of scout bees, g_limit=4 were used. We employed neutral, anionic and cationic monomers in all combinations, leading to overall neutral clusters, and saved 1000 local minima for each combination. All local minima were optimized at the GFN1-xTB level. ArbAlign (Temelso et al., 2017) was used to filter out identical clusters, based on root-mean-square deviation (RMSD) values. Based on previous studies on atmospheric clusters (Kildgaard et al., 2018b, a), an RMSD cutoff of 0.38 Å was chosen. As a large number of structures remained after ArbAlign, only the 100 cluster structures with the lowest electronic energy were selected for further optimization at the ωB97X-D/6-31++G(d,p) level of theory. The five structures with the lowest free energy were then selected for calculation of single-point energies at the DLPNO-CCSD(T₀)/aug-cc-pVTZ level.

2.3 The “cluster-of-functional-groups” approach

Small organic molecules were chosen to represent the functional groups that act as hydrogen bond donors (alcohol (CH₃OH) and peroxide (CH₃OOH)), as well as functional groups that act as hydrogen bond acceptors (ether (CH₃OCH₃), epoxide (C₂H₄O), aldehyde (CH₃CHO), ketone (CH₃COCH₃), acid anhydride $\left({\mathrm{CH}}_{\mathrm{3}}\mathrm{C}\right(=\mathrm{O}\left)\mathrm{OC}\right(=\mathrm{O}\left){\mathrm{CH}}_{\mathrm{3}}\right)$ and ester (COOCH₃)). Besides these groups, carboxylic acid, (HCOOH), which is both an acceptor and a donor, was also included. The studied functional groups were included in the funneling approach in order to identify the global minimum cluster structures. In this manner the functional groups should orient in the most favorable positions and can be viewed as ideal “contact points” for assembling the molecular backbone afterwards. Hence, the “cluster-of-functional-groups” approach should be viewed as a single proxy OOM. The strength of this approach is that we do not need to explicitly consider the origin of the OOM, i.e. anthropogenic (aromatics) or biogenic (isoprene/terpenes), or whether it is a monomer or an accretion product as the approach will inherently identify the structural patterns that are important directly based on their ability to participate in cluster formation (Elm et al., 2023). However, if the functional groups are optimally oriented far from each other, it could allude to the fact that the potential compound is an accretion product.

We note the caveat that the “cluster-of-functional-groups” approach assumes that the binding free energies of the individual groups are additive. This might not necessarily be the case for realistic atmospheric OOMs, where several effects potentially make the free energies deviate from additivity, and there can be expected some degree of cancellation of errors. For instance, the enthalpy contribution is expected to be more or less additive given that there are no intramolecular hydrogen bonds in the OOM. In addition, it is unlikely that the multiple moieties of any realistic OOM can simultaneously reach the ideal contact points without introducing some strain in the backbone, which will also lead to a higher enthalpy. The entropy contribution will most likely not be additive as a major contribution comes from the loss of translational and rotational degrees of freedom, i.e., the conversion into lower-entropy vibrational degrees of freedom. For each clustering functional group 6 high-entropy degrees of freedom are lost, while only a total of 6 high-entropy degrees of freedom are lost in the OOM. Despite these deficiencies, the “cluster-of-functional-groups” approach can still be employed for screening purposes and to yield some indication of which combinations of functional groups might potentially be important in atmospheric cluster formation.

Figure 1 The binding free energy of all organic–organic interactions between all combinations of hydrogen bond donor and acceptor, calculated at the DLPNO-CCSD(T₀)/aug-cc-pVTZ //ωB97X-D/6-31++G(d,p) level of theory with quasi-harmonic cutoff of 100 cm⁻¹, at 298.15 K and 1 atm. CA is short for carboxylic acid.

Figure 2 The binding free energy of clusters with one SA and one organic functional group, as well as with one SA fully saturated by organic functional groups, calculated at the DLPNO-CCSD(T₀)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory with quasi-harmonic cutoff of 100 cm⁻¹, at 298.15 K and 1 atm.

2.4 Atmospheric cluster dynamics

To study the time development of molecular cluster distributions, the free energy is used to solve the birth–death equations with ACDC (McGrath et al., 2012). The birth–death equations describe how clusters are destroyed and created by condensation and evaporation and are given as the change in concentration of cluster i,

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{\text{d}{c}_i}{\text{d}t}}& =\sum _{j<i}{\mathit{\beta }}_{j,(i-j)}{c}_j{c}_{(i-j)}+\sum _j{\mathit{\gamma }}_{(i+j)\to i}{c}_{(i+j)}\\ & -\sum _j{\mathit{\beta }}_{i,j}{c}_i{c}_j-\sum _{j<i}{\mathit{\gamma }}_{i\to j}{c}_i+Q_i-S_i.\end{array}\end{array}$$

Here, j is another cluster in the system, β_i,j is the collision coefficient between cluster i and j, and γ_i→j is the evaporation coefficient of cluster i into a smaller cluster, of which one is cluster j. Q_i covers outside sources of i and S_i other possible loss mechanisms for i. The collision coefficient is calculated using kinetic gas theory and the evaporation coefficient via mass balance based on calculated free energy,

$$\begin{array}{}\text{(2)}& {\mathit{\gamma }}_{(i+j)\to i}={\mathit{\beta }}_{i,j}{c}_{\text{ref}}\text{exp}\left({\displaystyle \frac{\mathrm{\Delta }{G}_{i+j}-\mathrm{\Delta }{G}_i-\mathrm{\Delta }{G}_j}{k_{\mathrm{b}}T}}\right),\end{array}$$

where c_ref is the monomer concentration of the vapor used to calculate the free energies. The birth–death equations are generated by checking all possible cluster configurations and examining which evaporations and collisions are able to destroy or create a given cluster (McGrath et al., 2012). As it is not possible to explicitly simulate the entire cluster size range, a user-specific size limit has to be chosen where the clusters are considered stable against evaporation and are counted towards the formation rate. We will refer to these clusters as “outgrowing”. For our systems, we chose the clusters with one additional acid compared to the cluster sizes we have data for [(SA)₃(base)₂, (SA)₂(OOM)₂(base)₂ and (SA)₃(OOM)₁(base)₂]. We exclude outgrowing collisions with bases as clusters composed of more bases than acids are usually not stable. The outgrowing clusters containing OOMs have one more molecule present compared to the pathway without OOM. This will slightly underestimate the contribution of OOM to the relative simulated J_potential rates. These outgrowing cluster sizes are still quite small and therefore artificially stabilize the systems as the critical cluster size is not necessarily captured well. To distinguish the calculated rates from actual nucleation rates, we will refer to the rates as “cluster formation potentials” (J_potential). These rates illustrate the potential of the cluster to grow to larger sizes and correspond to an upper-bound on the formation rate. We refer to the “Clusteromics I” paper (Elm, 2021a) for additional information on “cluster formation potentials”.

Figure 3 The binding free energy of the (SA)₁(base)₁(organic)₁ clusters, calculated at the DLPNO-CCSD(T₀)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory with quasi-harmonic cutoff of 100 cm⁻¹, at 298.15 K and 1 atm. Carboxylic acid, the base and sulfuric acid have been omitted in the labels for readability. REF is the reference (SA)₁(base)₁ cluster without organics present.

3 Results and discussion

3.1 Pure organic clusters

Initially, pure organic clusters were studied containing every possible donor–acceptor combination. This led to a total of 21 dimer cluster structures. The calculated free energies, at the DLPNO-CCSD(T₀)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory, are given in Fig. 1. The calculations are performed at 298.15 K and 1 atm.

All the dimer clusters have a positive change in free energy, except for the dimer cluster with two carboxylic acid (CA) groups that has a $\mathrm{\Delta }G=-\mathrm{3.23}$ kcal mol⁻¹. This value is in excellent agreement with the experimental value of −3.46 kcal mol⁻¹ for the formic acid dimer cluster given by Vander Auwera et al. (2007), at 296 K. However, this value corresponds to a high evaporation rate of the formic acid dimer. Thus, cluster formation involving only organic molecules is unlikely at 298.15 K and 1 atm, unless the molecule contains multiple carboxylic acid groups. In all cases the inclusion of a carboxyl group lowers the binding free energy considerably compared to the other functional groups. Hence, it is likely that other compounds besides organics are needed for forming the initial stable clusters.

Figure 4 The binding free energy of the (SA)₁(base)₁(CA)₁(organic)₁ clusters, calculated at the DLPNO-CCSD(T₀)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory with quasi-harmonic cutoff of 100 cm⁻¹, at 298.15 K and 1 atm. CA is short for carboxylic acid. The base and sulfuric acid have been omitted in the labels for readability. REF is the reference (SA)₁(base)₁ cluster without organics present.

Figure 5 The (SA)₁(DMA)₁(acid anhydride)₁(CA)₁ and (SA)₁(DMA)₁(CA)₂ cluster geometries lowest in free energy. Calculated at the DLPNO-CCSD(T₀)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory with quasi-harmonic cutoff of 100 cm⁻¹, at 298.15 K and 1 atm. White – hydrogen, brown – carbon, red – oxygen, yellow – sulfur, blue – nitrogen. CA is short for carboxylic acid.

3.2 SA–organic clusters

As sulfuric acid is known to be essential for new particle formation, dimer clusters involving all the different organic functional groups and SA are studied. We studied the (SA)₁(organic)₁ dimer clusters where only a single organic functional group is present, as well as the (SA)₁(CA)₁(organic)₁ clusters that have been fully saturated by organic functional groups, meaning all donor and acceptor groups in the SA molecule are hydrogen-bonded. This led to a total of 46 clusters studied. The fully saturated SA clusters all include at least one carboxylic acid group as well as two other identical organic groups, except for the case where there are two carboxylic acid groups in total as this fully saturates the SA molecule. The calculated free energies (at 298.15 K and 1 atm) are shown in Fig. 2.

As expected, including SA decreases the free energy for all the dimer clusters when compared to the pure organic dimer clusters. Generally, saturating the SA molecule with organic functional groups also results in a cluster lower in free energy. However, in cases where there are added two hydrogen bond donors to the (SA)₁(CA)₁ cluster, the addition of the second functional group leads to an increase in the binding free energy. For instance, this is the case for the clusters containing either two esters and two epoxides.

For the dimer clusters, (SA)₁(CA)₁ is the most stable with $\mathrm{\Delta }G=-\mathrm{5.06}$ kcal mol⁻¹, while for the clusters with a fully saturated SA, it is (SA)₁(CA)₂, with $\mathrm{\Delta }G=-\mathrm{9.51}$ kcal mol⁻¹, indicating a trend of carboxylic acid being the most stabilizing functional group. This finding is consistent with the previous study by Elm et al. (2017a). It should be noted that the (SA)₁(CA)₁(alcohol)₂ cluster has a binding free energy very close to that of the (SA)₁(CA)₂ cluster, with a value of $\mathrm{\Delta }G=-\mathrm{8.09}$ kcal mol⁻¹. This suggests that multiple alcohol groups could also be important for cluster formation.

Figure 6 The (SA)₁(DMA)₁(alcohol/aldehyde/peroxide)₁(CA)₁ cluster geometries lowest in free energy. Calculated at the DLPNO-CCSD(T₀)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory with quasi-harmonic cutoff of 100 cm⁻¹, at 298.15 K and 1 atm. White – hydrogen, brown – carbon, red – oxygen, yellow – sulfur, blue – nitrogen. CA is short for carboxylic acid.

3.3 (SA)₁(base)₁–organic clusters

Next, bases are added to the clusters, to gauge its effect on the cluster stability. Initially, simple trimer clusters involving one organic functional group, one SA and one base (either AM, MA, DMA or TMA) are studied. The calculated binding free energies of these clusters, at 298.15 K and 1 atm, are presented in Fig. 3.

In all cases the (SA)₁(base)₁(organic)₁ clusters are lower in free energy compared to the (SA)₁(base)₁ clusters. Hence, the addition of organics is thermodynamically favorable. For instance, the addition of a carboxylic acid group to the (SA)₁(base)₁ leads to a lowering of the free energy by −4.09, −7.80, −8.09 and −4.95 kcal mol⁻¹ for AM, MA, DMA and TMA, respectively.

Generally, the cluster binding free energies follow the basicity of the base with the following trend: DMA ≃ TMA < MA < AM. This is consistent with previous work, where it has been shown both theoretically (Kurtén et al., 2008) and experimentally (Jen et al., 2014; Glasoe et al., 2015; Almeida et al., 2013) that amines are much more effective in driving new particle formation with SA than AM is, even when accounting for the large difference in atmospheric mixing ratios. The ordering of DMA and TMA depends slightly on the system, with DMA dominantly leading to the most stable clusters. In particular, the (SA)₁(DMA)₁(CA)₁, (SA)₁(TMA)₁(CA)₁ and (SA)₁(DMA)₁(acid anhydride)₁ clusters have the lowest free energies, with ΔG of −19.42, −17.77 and −18.46 kcal mol⁻¹, respectively. Carboxylic acid groups are again found to form the most favorable interactions, with the cluster with acid anhydride only being about 1 kcal mol⁻¹ less stable.

To further explore which functional groups that could yield stable clusters, we added one more organic functional group to the (SA)₁(base)₁ clusters. As carboxylic acid has been prevalent in forming the most stable clusters in the previous sections, we added the constraint that one of the organics should be a carboxylic acid forming the (SA)₁(base)₁(CA)₁(organic)₁ clusters. The calculated binding free energies (at 298.15 K and 1 atm) are given in Fig. 4.

The addition of one more organic functional group has a stabilizing effect on all the clusters. The cluster with two carboxylic acids is found to be lowest in free energy, with $\mathrm{\Delta }G=-\mathrm{23.96}$ kcal mol⁻¹ for the (SA)₁(DMA)₁(CA)₂ cluster. However, the cluster that contains an acid anhydride is also very close in free energy, with $\mathrm{\Delta }G=-\mathrm{23.79}$ kcal mol⁻¹. The cluster structures can be seen in Fig. 5, calculated at the DLPNO-CCSD(T₀)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory. In both cases the organics can be seen to interact both with SA and the DMA molecules.

Furthermore, the (SA)₁(DMA)₁(alcohol)₁(CA)₁, (SA)₁(DMA)₁(aldehyde)₁(CA)₁ and (SA)₁(DMA)₁(peroxide)₁(CA)₁ clusters also have low free energy values with ΔG of −19.24, −18.01 and −20.86 kcal mol⁻¹, respectively. Figure 6 presents the cluster structures. Again the organics are interacting with both SA and the DMA. Interestingly, the aldehyde in the (SA)₁(DMA)₁(aldehyde)₁(CA)₁ cluster preferentially binds to the SA and DMA compounds via several weak interactions, instead of binding to the vacant S–OH group in SA. We found that the cluster where the aldehyde resides at the S–OH group in SA is 0.62 kcal mol⁻¹ higher in free than the presented structure in Fig. 5 (see Supplement). Thus, the ability to form clusters with organic molecules including multiple alcohols, aldehydes or peroxides as functional groups should be further studied in the future.

3.4 (SA)₁₋₂(base)₁₋₂(CA)₃ clusters

As the carboxylic acid groups consistently lead to the most stable clusters, the (SA)₁₋₂(base)₁₋₂(CA)₃ clusters have also been studied. The calculated binding free energies, at 298.15 K and 1 atm, are given in Fig. 7.

Figure 7 The binding free energy of the (SA)_1-2(base)₁(CA)₃ clusters, calculated at the DLPNO-CCSD(T₀)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory with quasi-harmonic cutoff of 100 cm⁻¹, at 298.15 K and 1 atm. Carboxylic acid has been omitted in the labels for readability.

Figure 8 The (SA)₁₋₂(DMA)₁₋₂(CA)₃ cluster geometries lowest in free energy. Calculated at the DLPNO-CCSD(T₀)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory with quasi-harmonic cutoff of 100 cm⁻¹, at 298.15 K and 1 atm. White – hydrogen, brown – carbon, red – oxygen, yellow – sulfur, blue – nitrogen. CA is short for carboxylic acid.

Logically, the inclusion of three carboxylic acid groups has a stabilizing effect, just as the inclusion of more SA and base does. Following the trend from earlier of DMA producing the most stable clusters, (SA)₂(DMA)₂(CA)₃ is the most stable cluster with $\mathrm{\Delta }G=-\mathrm{56.27}$ kcal mol⁻¹. Furthermore, DMA is found as a component in the three clusters that are lowest in free energy, with (SA)₂(MA)₁(DMA)₁(CA)₃ having $\mathrm{\Delta }G=-\mathrm{55.16}$ kcal mol⁻¹ and (SA)₂(DMA)₁(TMA)₁(CA)₃ having $\mathrm{\Delta }G=-\mathrm{52.89}$ kcal mol⁻¹. AM, on the other hand, is again the base that produces the clusters with the highest free energy, both alone and in combination with other bases. The two clusters highest in free energy are (SA)₁(AM)₂(CA)₃ with $\mathrm{\Delta }G=-\mathrm{21.40}$ kcal mol⁻¹ and (SA)₁(AM)₁(CA)₃ with $\mathrm{\Delta }G=-\mathrm{20.07}$ kcal mol⁻¹. Two SA molecules lead to more stable clusters than two DMA molecules, as (SA)₂(DMA)₁(CA)₃ has a slightly lower free energy compared to (SA)₁(DMA)₂(CA)₃ with a ΔG of −36.92 and −35.59 kcal mol⁻¹, respectively. While this is consistent with previous studies (Olenius et al., 2013; Elm, 2017), we see that the inclusion of the organic acids decreases the difference between the two paths. The clusters lowest in free energy all have two SA molecules and two bases, as the inclusion of more SA and base gives rise to more acid–base interactions, as is illustrated in Fig. 8.

Figure 9 Simulated cluster formation potentials (in clusters cm⁻³ s⁻¹) as a function of oxygenated organic molecule mixing ratio, in the lower concentration limit (a) and upper concentration limit (b). The simulations are performed at 278.15 K and 1 atm.

While the clusters become lower in free energy as they increase in size, it is not a guarantee that they are stable against evaporation in the atmosphere. For instance, the subsequent addition of more carboxylic acid groups to the (SA)₁(DMA)₁ cluster leads to a decrease of −8.04, −4.54 and −3.13 kcal mol⁻¹ for the first, second and third addition, respectively. Hence, if these are treated as individual molecules it is unlikely that they form stable clusters at realistic atmospheric conditions. However, if all three functional groups are treated as single OOM, the free energy for adding a hypothetical idealized tricarboxylic to the (SA)₁(DMA)₁ cluster is $\mathrm{\Delta }G=-\mathrm{15.71}$ kcal mol⁻¹. It should be noted that this approach neglects the loss of high-entropy translation-and-rotation degrees of freedom (see discussion in Sect. 2.3). A pragmatic approach to remedy this effect would be to model the single OOM as a cluster of the three carboxylic acid groups. Essentially, we want to compare the following two reactions:

$$\begin{array}{}\text{(3)}& {\displaystyle }\left(\text{SA}{)}_{\mathrm{1}}\right(\text{DMA}{)}_{\mathrm{1}}+\mathrm{3}\times \text{HCOOH}⟶\left(\text{SA}{)}_{\mathrm{1}}\right(\text{DMA}{)}_{\mathrm{1}}(\text{HCOOH}{)}_{\mathrm{3}}\text{(4)}& {\displaystyle }\left(\text{SA}{)}_{\mathrm{1}}\right(\text{DMA}{)}_{\mathrm{1}}+{\text{(HCOOH)}}_{\mathrm{3}}⟶\left(\text{SA}{)}_{\mathrm{1}}\right(\text{DMA}{)}_{\mathrm{1}}(\text{HCOOH}{)}_{\mathrm{3}}.\end{array}$$

Calculating the binding free energy difference between Reactions (3) and (4) we obtain

$$\begin{array}{}\text{(5)}& \mathrm{\Delta }\mathrm{\Delta }G=\mathrm{\Delta }G\left[{\text{(HCOOH)}}_{\mathrm{3}}\right]-\mathrm{\Delta }G[\mathrm{3}\times \text{HCOOH}].\end{array}$$

Computing this value at 298.15 K and 1 atm leads to a perfect cancellation of errors with a ΔΔG value of −0.06 kcal mol⁻¹. However, as the enthalpy and entropy terms are different for Reaction (3) and (4), we can expect a different temperature dependence. Recalculating ΔΔG at 278.15 K and 1 atm leads to a +1.5 kcal mol⁻¹ difference. Hence, there is very little difference between the two methods of calculating the binding free energies, and we will in the following stick to the simple method given by Reaction (3). Nevertheless, this is a very strong binding and would correspond to a quite stable cluster. Hence, by employing the “cluster-of-functional-groups” approach we have identified that tricarboxylic acids are likely candidates for forming stable clusters with SA and bases. This conclusion was already hypothesized but not explicitly proven in our earlier work (Elm et al., 2017a).

3.5 Cluster formation potentials

Based on the calculated thermochemistry in the previous sections, we can study how a hypothetical OOM composed of three carboxyl groups can potentially stabilize sulfuric acid–base clusters. Letting the three carboxyl groups represent a single OOM, we simulated the cluster formation potential (J_potential) of the (SA)₁₋₂(base)₁₋₂(OOM)₁ systems, with base = AM, MA, DMA and TMA. We calculated the collision coefficient of the OOM in the ACDC simulations based on the liquid density for formic acid and 3 times the mass of formic acid. To allow direct comparability, we conformed with our previous studies and used the same vapor concentrations as in the Clusteromics series of papers (Elm, 2021a, b, 2022; Knattrup and Elm, 2022; Ayoubi et al., 2023). Sulfuric acid was fixed at 10⁶ molecules cm⁻³, and the concentration of the bases was studied giving two extremes with a “lower limit” and an “upper limit”. These were set as follows: AM (10, 10 ppb), MA (1, 100 ppt) and DMA/TMA (1, 10 ppt). It should be noted that the low concentration limit likely best represents the actual concentrations observed in the ambient atmosphere.

The simulations were performed at 278.15 K and 1 atm using the settings described in Sect. 2.4. Figure 9 presents the simulated cluster formation potentials (J_potential) for the four bases as a function of OOM concentration (0–10 ppt).

It is seen that the influence of the added OOM is highly dependent on the base in the SA–base–OOM systems. In the lower concentration limit the J_potential value increases around 1 order of magnitude with 10 ppt OOM present for MA and DMA. This is roughly 4 orders of magnitude for the SA–AM–OOM systems. However, the increase is from $\mathrm{1.49}\times {\mathrm{10}}^{-\mathrm{10}}$ to $\mathrm{9.59}\times {\mathrm{10}}^{-\mathrm{6}}$ cm⁻³ s⁻¹, implying that the absolute cluster formation potential is still negligible. The cluster formation potential of the SA–TMA–OOM system is found to be more or less unaffected by the presence of the OOM and thereby does not increase much with increased OOM concentration. This is caused by the three bulky methyl groups in TMA, which impede the attachment of the OOM due to steric hindrance (See Supplement). The formation potentials present the following trend at 10 ppt of OOM: DMA > TMA > MA > AM. In the upper concentration limit we see very similar trends, except that the SA–AM–OOM system begins to overtake the SA–MA–OOM system, due to the very high mixing ratio of 10 ppb of AM.

Inspecting the fluxes out of the systems, in the low concentration regime, shows that the OOM (at 1 ppt) contributes to 97.96 %, 98.86 %, 72.44 % and 7.52 % of the outgrowing fluxes for the AM, MA, DMA and TMA systems, respectively. The fluxes for each system at various concentrations can be seen in the Supplement. Hence, the more weakly bound the SA–base system is without the OOM present, the more the OOM will contribute to the cluster formation potential. At 10 ppt of OOM 100 % of the outgrowing clusters contain an OOM for AM, MA and DMA, both in the lower and upper concentration regime of the bases. At 10 ppt of OOM the outgrowing clusters containing an OOM for TMA are 78.97 % and 87.22 % in the lower and upper concentration limits. Again, it should be noted that the simulated J_potential values show the potential to grow to larger sizes and not the actual nucleation rate. Hence, whether the OOM will evaporate from larger clusters or further contribute to the cluster growth remains unknown.

Overall, we see that the hypothetical OOM substantially contributes to the cluster formation potential, but as this is an idealized compound, explicit SA–base–tricarboxylic acid clusters should be further studied. Cluster formation with 3-methyl-1,2,3-butanecarboxylic acid (MBTCA) (Müller et al., 2012) has previously been carried out for SA–MBTCA clusters (Ortega et al., 2016; Elm et al., 2017b) and SA–AM–MBTCA clusters (Myllys et al., 2017), but it remains ambiguous whether MBTCA is important for cluster formation. Recent evidence has suggested that SA–MBTCA clusters primarily grow along the MBTCA coordinate, i.e., without the participation of SA (Elm, 2019b). Other potential atmospheric relevant tricarboxylic acids could be the 3-carboxyheptanoic acid (CHA) compound identified from limonene oxidation (Jaoui et al., 2006; Yasmeen et al., 2011). Recent work indicates that CHA might be involved in SA–CHA cluster formation (Tan et al., 2022), and it would be worth to explicitly study SA–base–CHA clusters further. Finally, there has been identified a large accretion product, known as pinyl diaterpenylic ester (PDPE), from α-pinene oxidation (Kristensen et al., 2013). Such a large compound has never been studied in atmospheric cluster formation and would be worth investigating in the future. This is also aligned with experimental observations that deem large accretion products as the most likely OOMs to drive nucleation (Lehtipalo et al., 2018; Dada et al., 2023).

4 Conclusions

Using quantum chemical calculations we have employed a new method, denoted the “cluster-of-functional-groups” approach, for studying the role of oxygenated organic molecules (OOMs) in organic enhanced atmospheric cluster formation. Studying all combinations of organic donor and acceptor functional groups and their intermolecular interactions with sulfuric acid–base clusters, it is found that carboxyl groups lead to the most stable clusters. Based on our findings, we study the cluster formation potential of a hypothetical tricarboxylic acid molecule, represented as a cluster of three carboxyl groups. We find that the proxy OOMs composed of three carboxylic acid groups are present almost in all the relevant outgrowing clusters. Thereby the OOM is directly participating in SA–base–OOM cluster formation. Hence, the role of explicit tricarboxylic acids should be further studied, and especially large flexible tricarboxylic acid accretion products might be important for cluster formation.

We hypothesize that a “cluster-of-functionals-groups” approach composed of multiple carbonyls, alcohols and hydroperoxide groups might rival the binding strength of three carboxylic acids in SA–base clustering and should be further studied in the future. The approach should also be extended to study ion-induced nucleation, which might substantially enhance pure organic nucleation. In addition, the effect of temperature should be investigated, as all the functional groups will bind more strongly at lower temperatures.

Overall, our presented “cluster-of-functional-groups” approach is general and applicable to other topics where multifunctional compounds are important and where the exact chemical structure remains unknown. Hence, we speculate that the work could also be extended to study the potential molecular structures that can enhance ice nucleation in the form of identifying potential functional groups that facilitate the formation of hexagonal ice I_h structures.