Atmospheric new particle formation enhanced by tricarboxylic acids

Organic molecules contribute substantially to the formation of aerosols in the atmosphere, forming what is known as secondary organic aerosols (SOAs). The organic molecules are emitted as volatile organic compounds (VOCs) and undergo a number of reactions in the atmosphere. Due to the variety of both VOCs and reaction pathways, it has been difficult to elucidate the exact structure of an organic molecule that is able to drive new particle formation (NPF). Using quantum chemistry methods, we have studied the NPF ability of three different oxygenated organic molecules (OOMs): 3-methyl-1,2,3-butanecarboxylic acid (MBTCA), carboxyheptanoic acid (CHA), and pinyl diaterpenylic ester (PDPE). These all contain three carboxylic acids, which, as previous works suggest, are good candidates for driving NPF and have been observed in the atmosphere, as well as in lab experiments. Using computational methods, we studied the (OOM)_1–2(SA)_0–2(base)_0–2 clusters, where SA is sulfuric acid, and the base is [ammonia (AM), methylamine (MA), dimethylamine (DMA), and trimethylamine (TMA)]. Geometry optimization and thermochemical parameters are calculated at the ωB97X-D/6-31++G(d,p) level of theory, and single-point energies are calculated at the DLPNO-CCSD(T₀)/aug-cc-pVTZ level of theory. We found that PDPE was able to produce the most stable clusters, presumably due to its higher flexibility compared to MBTCA and CHA. Cluster formation potentials are simulated using the Atmospheric Cluster Dynamics Code (ACDC). We found that all three OOMs were able to enhance cluster formation for the (OOM)(SA)(base) systems by 2–3 orders of magnitude for most systems. In particular, the (OOM)(SA)(DMA) system has a high cluster formation potential, with similar trends in the enhancement across all three OOMs.

1 Introduction

Atmospheric aerosols can act as cloud condensation nuclei (CCN), increasing the global albedo (Boucher and Lohmann, 1995), or can scatter light themselves (Canadell et al., 2021). This is known as the indirect and direct aerosol effect, respectively, both of which have a net cooling effect on the global climate (Canadell et al., 2021). Aerosols are still responsible for the largest uncertainty in modern radiative forcing models according to the Intergovernmental Panel on Climate Change (IPCC) (Canadell et al., 2021). Atmospheric aerosols are either directly emitted into the atmosphere as primary aerosols or emerge through the clustering of molecules initially emitted into the atmosphere as gaseous species. The latter process is known as new particle formation (NPF) and is the mechanism responsible for secondary aerosols (Kulmala et al., 2013). Usually, the gaseous species have to undergo oxidation, lowering their volatility, before NPF is favorable. On average, it is estimated that, globally, about 50 % of CCN arise from NPF (Boucher and Lohmann, 1995; Merikanto et al., 2009), with a higher fraction over the eastern United States and a lower fraction over mainland Europe (Zhao et al., 2024).

It is well known that inorganic acids and bases are important for the initial cluster formation over land (Kulmala et al., 2013). Some confirmed aerosol precursors are sulfuric acid (SA) (Sipilä et al., 2010), water (Loukonen et al., 2010; Temelso et al., 2012a, b; Henschel et al., 2014; Rasmussen et al., 2020; Neefjes et al., 2026), and bases such as ammonia (AM) (Weber et al., 1996; Kirkby et al., 2011; Elm, 2021b; Dunne et al., 2016), as well as amines such as methylamine (MA), dimethylamine (DMA), and trimethylamine (TMA) (Almeida et al., 2013; Elm et al., 2020; Engsvang et al., 2023; Kurtén et al., 2008; Loukonen et al., 2010; Nadykto et al., 2011, 2015, 2014; Jen et al., 2014; Glasoe et al., 2015; Elm, 2021b, b; Kurtén et al., 2008; Nadykto et al., 2011; Jen et al., 2014; Nadykto et al., 2015; Glasoe et al., 2015). Over the oceans, it is the oxidation products of dimethylsulfide, such as methanesulfonic acid (Barnes et al., 2006), as well as various iodine species, that are believed to drive NPF together with SA (O'Dowd et al., 2002; Sipilä et al., 2016; He et al., 2021, 2023). Organic molecules make up a large part of the total atmospheric aerosol load (Jimenez et al., 2009; Hallquist et al., 2009), with secondary organic aerosols (SOAs) making up 15 %–80 % by mass of global PM_2.5 (Hallquist et al., 2009). Zhang et al. (2004) showed that aromatic acids enhance nucleation of SA due to the formation of a complex between SA and the aromatic acid causing a lowered nucleation barrier. In addition, multiple chamber studies have proved that this is also true for oxidation products of biogenic terpenes, like α-pinene and Δ-3-carene (Hoffmann et al., 1997; Pathak et al., 2007; Laj et al., 2009; Zhang et al., 2009; Metzger et al., 2010; Schobesberger et al., 2013; Riccobono et al., 2014; Quéléver et al., 2019; Simon et al., 2020; Thomsen et al., 2021, 2022, 2024). SOAs are especially important in rural sites, where the concentration of organics is high compared to the concentration of SA (Fang et al., 2020; Ehn et al., 2014). Fang et al. (2020) suggested that dicarboxylic acids take part in NPF, given that the observed nucleation rates in their field study could not be explained by SA–AM or SA–DMA clustering alone.

Clustering of organic molecules has been studied for almost 20 years using quantum chemical methods, starting with Nadykto and Yu (2007), who studied the interaction between AM or water and simple, common carboxylic acids such as formic and acetic acid. More recently, larger oxygenated organic molecules, including multiple alcohol and/or carboxyl groups, were studied computationally, showing binding Gibbs free energies as low as −14.37 kcal mol⁻¹ for dimer clusters of dicarboxylic acids (Tan et al., 2022).

Given the vast number of different volatile organic compounds (VOCs) emitted into the atmosphere and the complexity of the possible reaction pathways, the area of organic nucleation has proven to be difficult to study. Elm et al. (2023) concluded that, despite SOAs having been studied extensively over many years, both experimentally and computationally, the exact structure of an organic molecule that is able to drive nucleation has not yet been identified. With the limited knowledge available, it is not yet clear whether organics are important for nucleation or only enter the particle after it has been formed, contributing to particle growth (Kulmala et al., 2013). Carboxylic acids are believed to be the most promising organic functional group to be involved in new particle formation (Elm et al., 2017). In our previous work using the cluster-of-functional-groups approach (Pedersen et al., 2024), we showed that three different formic acid molecules, representing a tricarboxylic acid, should form stable clusters with SA and bases, where the bases are AM, MA, DMA, and TMA. In particular, DMA showed increased cluster formation rates when the concentration of the oxygenated organic molecule (OOM) was increased due to its high base strength.

Tan et al. (2022) studied eight different organic–SA systems and found that the carboxyheptanoic acid–SA systems, as well as the phtalic acid–SA systems, had high cluster formation potentials, suggesting that organic molecules can indeed play a role in nucleation. Large accretion products, also called dimers, are believed to be able to drive organic nucleation (Lehtipalo et al., 2018; Dada et al., 2023). These are believed to be formed by gas-phase cross reaction of two peroxy radicals, forming an accretion product (Peräkylä et al., 2023). However, Kenseth et al. (2023) studied the structures and formation mechanisms of four different dimeric accretion products found in atmospheric particles and concluded that accretion products, specifically those containing cis-pinic acid subunits, are formed through particle-phase nucleophilic addition reactions.

Figure 1 Molecular structure of MBTCA, CHA, and PDPE.

In this article, we study the nucleating ability of three different organic molecules, specifically 3-methyl-1,2,3-butanecarboxylic acid (MBTCA), carboxyheptanoic acid (CHA), and pinyl diaterpenylic ester (PDPE) (see Fig. 1). These have been chosen as they all contain three carboxylic acid moieties and have been observed experimentally. MBTCA is formed from OH-initiated oxidation of α-pinene and has been observed in chamber experiments and in field campaigns in the Sierra Nevada Mountains in California (Müller et al., 2012; Kristensen et al., 2013). PDPE, a large accretion product, was found to be formed from α-pinene oxidation using O₃ (Kristensen et al., 2016) and was also found in the same field study in the Sierra Nevada mountains (Kristensen et al., 2013), as well as in a field study in Hyytiälä, Finland (Kristensen et al., 2016). CHA was identified in a chamber experiment as a product of the photooxidation of a d-limonene–NO_x–air mixture and has also been detected in field studies in Hungary (Jaoui et al., 2006; Yasmeen et al., 2011). The clustering ability of both MBTCA and CHA has already been studied using quantum chemical methods (Tan et al., 2022), and it was found that trimers made up of both SA and the organic molecule were stable enough, leading to negligible evaporation.

The clustering ability of MBTCA, CHA, and PDPE is studied, both in purely organic clusters and in (SA)_1–2(base)_1–2 clusters, where the bases are ammonia (AM), methylamine (MA), dimethylamine (DMA), and trimethylamine (TMA). These have been shown to enhance nucleation in clusters containing SA (Kulmala et al., 2013; Kirkby et al., 2011; Schobesberger et al., 2013; Weber et al., 1996; Dunne et al., 2016; Almeida et al., 2013; Elm et al., 2020; Engsvang et al., 2023; Kurtén et al., 2008; Loukonen et al., 2010; Nadykto et al., 2011, 2015, 2014; Jen et al., 2014; Glasoe et al., 2015; Elm, 2021b). The clustering ability was investigated by calculating the binding Gibbs free energy at the DLPNO-CCSD(T₀)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory using an extensive configurational sampling workflow to find the global free energy minimum structure. The thermodynamic data were used as input in the Atmospheric Cluster Dynamics Code (ACDC) to calculate cluster formation rates for the (SA)_1–2(base)_1–2(OOM)_1–2 systems.

2 Methods

2.1 Computational details

ABCluster (Zhang and Dolg, 2015) was used for an initial configurational search with the CHARMM force field using a population size SN=3000, maximum generations g_max=200, and a number of scout bees g_limit=4, as recommended by Kubečka et al. (2019). A total of 1000 local minima were saved for each protonation state. Semi-empirical GFN1-xTB energy calculations and geometry optimizations were executed using the XTB 6.4.0 program (Grimme et al., 2017). A further search of the potential energy surface (PES) was performed using CREST 2.12 (Pracht et al., 2020) at the GFN1-xTB level within the iterative metadynamics genetic structure crossing (iMTD-GC) workflow. An energy window of 30 kcal mol⁻¹ was used, accounting for the differences in the potential energy surface (PES) at the semi-empirical level and higher levels of theory, to sample the larger configurational space. DFT calculations with the ωB97X-D functional (Chai and Martin, 2008) and 6-31++G(d,p) basis set were performed using Gaussian16, version B.01 (Gaussian 16, 2016), using default convergence criteria. DLPNO–CCSD(T₀) (Riplinger and Neese, 2013; Riplinger et al., 2013) single-point energies were calculated in ORCA 5.0.4 (Neese, 2012, 2022), with the aug-cc-pVTZ basis set, using the TightSCF convergence criteria and the NormalPNO setting (Liakos et al., 2015). Both the workflow and data processing were automated with JKCS (Kubečka et al., 2023a).

2.2 Configurational sampling workflow

For the configurational search, a funnel-type procedure was employed (Temelso et al., 2018; Odbadrakh et al., 2020; Kubečka et al., 2023a), where the level of theory is increased as the number of cluster candidates is decreased:

$$\begin{array}{rl}& \mathrm{ABCluster}\to \text{GFN1-xTB}\to \mathrm{CREST}\to \mathit{\omega }\text{B97X-D}\\ & \to \mathrm{DLPNO}.\end{array}$$

The initial cluster structures were generated using ABCluster. These were optimized using GFN1-xTB. The structure with the lowest electronic energy at this level was used as input for additional configurational sampling using CREST. The molecular dynamics and metadynamics within CREST allow for thorough exploration of the PES, including the internal rotations of bonds for flexible molecules, which were not accommodated for in ABCluster. The 100 cluster structures with the lowest electronic energy after the CREST calculation were selected to be optimized at the DFT-level. Single-point energies using DLPNO-CCSD(T₀)/aug-cc-pVTZ were calculated for the five clusters with the lowest electronic energy at the DFT-level, as recommended by Jensen et al. (2022). The binding Gibbs free energy ΔG_binding is calculated as

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }{G}_{\mathrm{binding}}=G_{\mathrm{cluster}}-\sum _i{G}_{\mathrm{monomer},i}.\end{array}$$

The quasi-harmonic approximation (Grimme, 2012) is implemented, with a threshold value of 100 cm⁻¹, to mitigate unphysical stabilizing effects from low-lying vibrational frequencies.

2.3 Atmospheric cluster dynamics code

The ACDC workflow (McGrath et al., 2012) generates and solves the birth–death equations, giving the change in concentration of cluster i as a function of condensation and evaporation,

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{c}_i}{\mathrm{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}$$

where j is another cluster or monomer 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 smaller clusters, one of which is cluster j. Q_i covers outside sources of i, and S_i other possible loss mechanisms for i. The collision coefficient is given by kinetic gas theory, and the evaporation coefficient is given by mass balance based on the calculated Gibbs free energy,

$$\begin{array}{}\text{(3)}& {\mathit{\gamma }}_{(i+j)\to i}={\mathit{\beta }}_{i,j}{c}_{\mathrm{ref}}\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 reference pressure (1 atm).

Due to the computational cost, it is not possible to simulate the cluster formation all the way from gas species to fully stable particles. Therefore, a limit has been chosen where the clusters are assumed to be stable against evaporation and are counted towards the cluster formation rate. These clusters will be referred to as “outgrowing”. In this study, we explicitly calculated thermodynamic data for the SA–OOM–base clusters up to size $\mathrm{2}\times \mathrm{2}\times \mathrm{2}$, while the data for the SA–base clusters up to size 3×3 were obtained from Kubečka et al. (2023b). The outgrowing clusters are therefore chosen to be one monomer larger than we have data for: (SA)₄(base)₃, (SA)₃(OOM)₂(base)₂, (SA)₂(OOM)₃(base)₂, and (SA)₃(OOM)₁(base)₃. Given that these outgrowing clusters are fairly small, the systems are artificially stable, and the critical cluster size may not yet have been reached, leading to the calculated formation rates that are likely overestimated. Furthermore, given that these clusters do not have an equal number of SA and base monomers, this might artificially overestimate the enhancement factor of the OOM. To test this, we performed a simulation for smaller clusters with (SA)₂(PDPE)₃(DMA)₂ and (SA)₃(PDPE)₀(DMA)₂ at 278.15 K, 1 ppt DMA, an SA concentration of 10⁶ molec. cm⁻³, and PDPE at 0 and 10 ppt. We find an enhancement factor of 137 compared to 396 with the other box sizes. The enhancement is still of the same order of magnitude, illustrating that the enhancement from the OOM is not driven by the boundary conditions. In order to distinguish the calculated rates from actual nucleation rates, we denote them as cluster formation potentials (J_potential) as they represent the potential of clusters to grow to larger sizes and should be interpreted as an upper limit of the nucleation rate. A more in-depth discussion of cluster formation potentials is given in Clusteromics I (Elm, 2021b). Finally, we defined the main outgrowing paths as the pathway following the highest flux backwards from the outgrowing clusters towards the monomers.

3 Results and discussion

3.1 Cluster thermochemistry

The limiting step in atmospheric cluster formation is often the formation of the initial dimer cluster (Olenius et al., 2013; Elm et al., 2017). Usually, organics bind too weakly to themselves and SA to initiate the initial dimer cluster formation. The lowest calculated binding Gibbs free energies of the three purely organic dimer clusters, (MBTCA)₂, (CHA)₂, and (PDPE)₂, are given in Fig. 2. The (SA)₁(DMA)₁ cluster is included as a reference because this system has been studied extensively and has been shown to form stable clusters (Kurtén et al., 2008; Kubečka, 2021).

Figure 2 The (OOM)₂ cluster geometries lowest in binding Gibbs free energy at the DLPNO-CCSD(T₀)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory with a quasi-harmonic threshold of 100 cm⁻¹ at 298.15 K and 1 atm. Orange is carbon, red is oxygen, and white is hydrogen.

The MBTCA–MBTCA and CHA–CHA interactions are both slightly weaker than the SA–DMA interaction, but the difference is very small. The PDPE–PDPE interaction is the strongest, even stronger than the SA–DMA interaction, suggesting a high probability for PDPE to drive NPF. Organic clusters have been extensively studied in the literature, with much emphasis on di- and/or tri-carboxylic acids. For instance, Elm et al. (2019) studied 13 dicarboxylic acid dimer clusters at the same level of theory as used in this study. They found that the glutaric acid (five-carbon backbone) dimer has a binding Gibbs free energy of −11.6 kcal mol⁻¹, the suberic acid (six-carbon backbone) dimer cluster has a binding Gibbs free energy of −9.9 kcal mol⁻¹, and the pimelic acid (seven-carbon backbone) dimer has a binding Gibbs free energy of −13.0 kcal mol⁻¹, while the rest of the dimer clusters they studied had significantly higher binding Gibbs free energies.

Other dicarboxylic acids originating from biogenic sources have been studied. For the pinic acid dimer cluster (Elm et al., 2014), a binding Gibbs free energy of −7.05 kcal mol⁻¹ was found, calculated at the M06-2X/6-311++G(3df,3pd) level of theory. Malonic acid (MOA)₂ yields $\mathrm{\Delta }G=-\mathrm{8.29}$ kcal mol⁻¹ with the RI-MP2/cc-pVTZ//PW91PW91/6-311++G(2d,2p) level of theory (Wang et al., 2021) or $\mathrm{\Delta }G=-\mathrm{4.35}$ kcal mol⁻¹ at the M06-2X/6-311++G(3df,3pd) level of theory (Zhang et al., 2018). This illustrates that the applied level of theory has a major impact on the modeled values, and one has to be cautious, especially with the applied single-point energy. While some of the previously studied dicarboxylic acids rival the binding Gibbs free energy values of the tricarboxylic acids in Fig. 2, they also lead to “dead-end” closed structures where there is no place for more molecules to attach. In order to facilitate cluster formation the organics need to bind strongly enough to the initial cluster, as well as enhance the attachment of additional molecules to the cluster.

In Fig. 2, the cluster geometries of the three (OOM)₂ clusters are shown. Here it can be seen that the short carbon backbone in both MBTCA and CHA restricts their ability to form three carboxylic acid–carboxylic acid bonds in the dimer clusters. (CHA)₂ is more stable than (MBTCA)₂ due to its less branched – and therefore longer – carbon backbone. This is also consistent with the study of dicarboxylic acids by Elm et al. (2019), where suberic acid dimers (seven-carbon backbone) were more stable than glutaric acid dimers (five-carbon backbone). In (MBTCA)₂, the two carboxylic acid pairs are placed closer together, lowering its stability. (PDPE)₂, on the other hand, has three carboxylic acid pairs, fully utilizing its hydrogen-bonding potential. This does mean that there is no obvious places for growth for the PDPE cluster, and expanding the cluster would therefore involve some restructuring. The (MBTCA)₂ and (CHA)₂ clusters will be able to grow more easily, without major change in their geometry, due to the two unbound carboxylic acids in each cluster. However, introducing other nucleation precursors such as SA and bases into the clusters might change their structure and lead to additional hydrogen-bond moieties that support further growth.

Given that the acid–base interaction between SA and AM (Kirkby et al., 2011) or amines (Almeida et al., 2013; Elm, 2021b) is known to produce clusters, all combinations of (SA)_1–2(base)_1–2(OOM)_1–2 were also studied. The calculated binding Gibbs free energies of the (SA)_1–2(base)_1–2(OOM)₁ clusters are given in Fig. 3.

Figure 3 The binding Gibbs free energy of the (SA)_1–2(base)_1–2(MBTCA)₁, (SA)_1–2(base)_1–2(CHA)₁, and (SA)_1–2(base)_1–2(PDPE)₁ clusters, calculated at the DLPNO-CCSD(T₀)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory with a quasi-harmonic threshold of 100 cm⁻¹ at 298.15 K and 1 atm. The left panel shows clusters with 1 SA, and the right panel shows clusters with 2 SAs.

It is seen that all the (SA)_1–2(base)_1–2(OOM)₁ clusters have very similar free energies for all three OOMs within each cluster type. Not surprisingly, the overall trend shows that the additions of more SA and base both act to lower the binding Gibbs free energy of the clusters. This is in accordance with what was previously seen when the OOMs were represented by three formic acid molecules (Pedersen et al., 2024). Similarly to the pure organic dimers, PDPE was again able to produce the most stable cluster. Overall, PDPE is only the strongest binding OOM in 7 out of the 16 clusters, with CHA being the strongest binding OOM in the remaining clusters. This is in contrast to Kähärä et al. (2025), who studied a total of 143 dimers consisting of different OOMs, including accretion products. They found that inflexible molecules formed more stable molecules than flexible molecules, even if they had a low bulk saturation vapor pressure, due to their decreased likeliness to form internal hydrogen bonds as monomers. Out of the three OOMs studied here, only PDPE shows a single internal hydrogen bond. The monomer structures are included in Fig. S1 in the Supplement. This is not the definitive parameter that determines the stability of the clusters, and the flexibility of PDPE is thus favorable. Within all cluster sizes and across all OOMs, DMA is the strongest binding base. For the smallest clusters, PDPE is generally unfavorable. This could be due to its large size and high flexibility, which make it able to position the carboxylic acid groups in the most favorable positions but introduce more steric hindrance for the smaller clusters. The (SA)₁(DMA)₁(OOM)₁ and (SA)₂(DMA)₂(OOM)₁ clusters are shown in Fig. 4. As can be seen in Fig. 4, both MBTCA and CHA have an available carboxylic acid group in the (SA)₁(DMA)₁(OOM)₁ clusters, while the larger PDPE molecule is able to make use of all three carboxylic acid groups when interacting with the SA and bases. However, as stated above, this cluster is 1 kcal mol⁻¹ less stable than the (SA)₁(DMA)₁(CHA)₁ cluster due to steric hindrance introduced by the large PDPE molecule. In the (SA)₂(DMA)₂(OOM)₁ clusters, all three OOMs have three binding carboxyl groups. This means that PDPE is able to stretch out, minimizing the steric hindrance while maximizing the number of favorable interactions, yielding a binding Gibbs free energy of −61.7 kcal mol⁻¹. For the (SA)₂(DMA)₂(MBTCA)₁ (−57.4 kcal mol⁻¹) and (SA)₂(DMA)₂(CHA)₁ (−60.1 kcal mol⁻¹) clusters, the smaller OOM means the inorganic acids and bases are closer together to facilitate hydrogen bonding, which increases steric hindrance compared to the (SA)₂(DMA)₂(PDPE)₁ cluster.

Figure 4 The (SA)₁(DMA)₁(OOM)₁ and (SA)₂(DMA)₂(OOM)₁ cluster geometries lowest in binding Gibbs free energy at the DLPNO-CCSD(T₀)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory with a quasi-harmonic threshold of 100 cm⁻¹ at 298.15 K and 1 atm. Orange is carbon, red is oxygen, blue is nitrogen, yellow is sulfur, and white is hydrogen.

Figure 5 The binding Gibbs free energy of the (SA)_1–2(base)_1–2(MBTCA)₂, (SA)_1–2(base)_1–2(CHA)₂, and (SA)_1–2(base)_1–2(PDPE)₂ 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. The left panel shows clusters with 1 SA, and the right panel shows clusters with 2 SAs.

The calculated binding Gibbs free energies of the (SA)_1–2(base)_1–2(OOM)₂ clusters are given in Fig. 5. With the addition of one more OOM, the trends remain roughly the same, with a ∼ 15 kcal mol⁻¹ decrease in free energy across all clusters. However, DMA is no longer the strongest binding base across the board. For example, the (SA)₁(MA)₂(PDPE)₂ cluster has a binding Gibbs free energy of −50.6 kcal mol⁻¹, while the (SA)₁(DMA)₂(PDPE)₂ cluster has a binding Gibbs free energy of −50.0 kcal mol⁻¹. However, it should be noted that the 0.6 kcal mol⁻¹ difference is within the uncertainty of the DLPNO-CCSD(T₀) energies. Generally, PDPE has been able to take advantage of how strongly it binds to itself as PDPE is the strongest-binding OOM in 11 out of the 16 clusters with two OOMs. For the remaining clusters, CHA is the strongest-binding OOM, while MBTCA is the weakest-binding OOM for all clusters, as was the case for the clusters with one OOM.

Figure 6 The (SA)₁(DMA)₁(OOM)₂ and (SA)₂(DMA)₂(OOM)₂ cluster geometries lowest in binding Gibbs free energy at the DLPNO-CCSD(T₀)/aug-cc-pVTZ//ωB97X-D/6-31++G(d,p) level of theory with a quasi-harmonic threshold of 100 cm⁻¹ at 298.15 K and 1 atm. Orange is carbon, red is oxygen, blue is nitrogen, yellow is sulfur, and white is hydrogen.

The (SA)₁(DMA)₁(OOM)₂ and (SA)₂(DMA)₂(OOM)₂ clusters are shown in Fig. 6 as DMA was the most favorable base for these cluster sizes. In both cases, PDPE is the most favorable OOM. As seen in Fig. 6, the two PDPE molecules enclose a core consisting of the inorganic acids and bases. Although these look close to being spherical and particle-like, these cluster structures are very flat, missing the depth to be fully considered a particle. For the two clusters with CHA, it is evident that the carbon backbone on CHA is not long enough to reproduce this shell structure; instead, the inorganic acids and bases are pushed to one side of the cluster. This is especially noticeable in the (SA)₂(DMA)₂(CHA)₂ cluster and can also be seen in the (SA)₂(DMA)₂(MBTCA)₂ cluster. However, both the (SA)₂(DMA)₂(CHA)₁ and (SA)₂(DMA)₂(CHA)₂ clusters have binding free energies that are very similar to the (SA)₂(DMA)₂(PDPE)₁ and (SA)₂(DMA)₂(PDPE)₂ clusters, respectively, especially compared to the two clusters containing MBTCA. For the (SA)₁(DMA)₁(MBTCA)₂ cluster, the inorganic acid and base is sandwiched between the two MBTCA molecules, which means that two carboxylic acid groups remain unbound. This results in a cluster that is noticeably less stable (−35.3 kcal mol⁻¹) than the (SA)₁(DMA)₁(CHA)₂ (−42.9 kcal mol⁻¹) and (SA)₁(DMA)₁(PDPE)₂ (−45.8 kcal mol⁻¹) clusters. The diameter of the (SA)₂(DMA)₂(PDPE)₂ cluster is above 1 nm, and the diameters of both the (SA)₂(DMA)₂(CHA)₂ and (SA)₁(DMA)₁(PDPE)₂ are close to 1 nm in diameter, making them near the lower limit of experimental measurement techniques (Vanhanen et al., 2011). The cluster radii of the clusters in Figs. 4 and 6 are given in Table S1.

3.2 Cluster formation potentials

From the thermochemical data in Sect. 3.1, we can calculate how the three OOMs, CHA, MBTCA, and PDPE, are likely to be able to stabilize an SA–base cluster. The simulated formation potentials (J_potential) of the (SA)_1–2(base)_1–2(OOM)_1–2 systems, with the base being AM, MA, DMA, and TMA and OOM being CHA, MBTCA, and PDPE, are given in Figs. 7–9. To allow for direct comparison, the vapor concentrations used are the same as those in the Clusteromics I–V series of papers (Elm, 2021b, a, 2022; Knattrup and Elm, 2022; Ayoubi et al., 2023; Pedersen et al., 2024). The sulfuric acid concentration was fixed at 10⁶ molec. cm⁻³, and the concentrations of the bases were studied for two extremes with a “lower limit” and an “upper limit”. These were set as follows: AM (10 ppt, 10 ppb), MA (1 ppt, 100 ppt), and DMA/TMA (1 ppt, 10 ppt), with the low concentration limit likely being the best representation of the actual concentrations observed in the ambient atmosphere. The OOM concentration was varied from 0 to 10 ppt, where 10 ppt should be considered to be an upper-bound estimate, only possible with close to 1 ppb concentration of the precursor, about 10 % yield of the tricarboxylic acid, and a condensation sink below 0.001 s⁻¹. The simulations were performed at 278.15 K and 1 atm, reflecting the conditions of springtime boreal forest areas, using the settings described in Sect. 2.3.

Figure 7 Simulated cluster formation potentials (in clusters cm⁻³ s⁻¹) as a function of CHA mixing ratio at the lower-concentration limit (left) and upper-concentration limit (right). The simulations are performed at 278.15 K and 1 atm.

Figure 8 Simulated cluster formation potentials (in clusters cm⁻³ s⁻¹) as a function of MBTCA mixing ratio at the lower-concentration limit (left) and upper-concentration limit (right). The simulations are performed at 278.15 K and 1 atm.

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

Given that, for the OOMs, only the density of MBTCA was available in the literature, the effect of changing density was tested on the system, yielding the highest nucleation rates (10 ppt PDPE, [SA] = 10⁶ molec. cm⁻³ and 10 ppt DMA) and the lowest nucleation rates (10 ppt MBTCHA, [SA] = 10⁶ molec. cm⁻³ and 10 ppt AM). When the density was set to that of agaric acid, a tricarboxylic acid with a relatively low density of 1.115 g cm⁻³, the cluster formation potential was 260.71 and $\mathrm{3.99}\times {\mathrm{10}}^{-\mathrm{7}}$ cm⁻³ s⁻¹ for lowest and highest nucleation rates, respectively. When the density was set to that of citric acid, another tricarboxylic acid, but with a relatively high density of 1.665 g cm⁻³, the cluster formation potential was 206.54 and $\mathrm{3.22}\times {\mathrm{10}}^{-\mathrm{7}}$ cm⁻³ s⁻¹, respectively. This is a factor of 1.26 and 1.24 for lowest and highest nucleation rates, respectively. If we assume that the error in the cluster formation potential is directly given by the error in the evaporation rate, the factor by which the rate changes is given as $f=\exp (\mathrm{\Delta }{G}_{\mathrm{error}}/\mathrm{RT})$. This means that, at room temperature, where RT≈0.59 kcal mol⁻¹, an ln (1.26)×0.59 kcal mol⁻¹ ≈ 0.14 kcal mol⁻¹ error would yield the same change. The exponential growing error in the binding Gibbs free energies is therefore expected to be more severe than the error in the assumed density, and we therefore set the density for all three molecules to that of MBTCA, 1.430 g cm⁻³ (Kostenidou et al., 2018).

Across all three OOMs, the trends are very similar, showing an increase in cluster formation potential as the OOM concentration increases, indicating that OOM enhances the nucleation rate. This is as expected since the OOM lowers the binding energy sufficiently to enhance the cluster formation potential, and introducing more OOM will give more potential for binding, leading to more particles. However, the extent of this enhancement is largely dependent on the specific base in the SA–base–OOM system and, to some extent, is also dependent on the base concentration, as was also evident in our previous work (Pedersen et al., 2024). The largest enhancement is seen in the SA–AM–OOM systems, especially those in the lower concentration limit. We hypothesize that this is due to the weak interaction between SA and AM, meaning that the OOM has a relatively large stabilizing effect on the cluster as it is otherwise weakly bound. However, since the increase in the formation potential is of the magnitude of 10⁻¹⁴ to 10⁻⁷, 10⁻³, and 10⁻⁴ for the SA–AM–MBTCA, SA–AM–CHA, and SA–AM–PDPE, respectively, the formation potentials for these systems are still negligible in all cases.

In agreement with the thermochemical data, the systems with the largest formation potentials are the SA–DMA–OOM systems. Interestingly, the specific OOM does not have a noticeable impact on the cluster formation potentials of these systems, with the J_potential value increasing by 2 orders of magnitude from 1 to 10² in the upper-concentration regime and from 10⁻¹ to 10¹ in the lower-concentration regime. The lack of noticeable differences in the enhancement between these systems highlights the importance of the explicit functional groups rather than that of the specific molecule itself. This could also indicate that it might be possible to lump the compounds into groups based on their constituent functional groups. However, for the remaining bases, the trends are not as identical across the three OOMs as they are for the SA–DMA–OOM system. While the majority of the systems follow the formation potential trend of DMA > TMA > MA > AM, the SA–base–CHA system at the upper concentration limit has the formation potential of the SA–AM–CHA system, overtaking that of SA–MA–CHA at a CHA concentration of 4.0 ppt. Similarly, the formation potential of the SA–AM–PDPE system catches up to that of SA–MA–PDPE at a PDPE concentration of 10 ppt.

Examining the fluxes reveals that OOMs are present in all outgrowing clusters. For all systems where the base is either AM, MA, or DMA, the OOM contributes to over 90 % of the outgrowing clusters. This was also the case for the SA–TMA–CHA and SA–TMA–PDPE systems in the lower concentration regime. However, at 0.5 ppt in the high-base-concentration regime, CHA contributes 51.67 % to the outgrowing clusters, while PDPE contributes 44.33 %, and MBTCA only contributes 38.44 %. MBTCA also had a low contribution of only 45.46 % in relation to the outgrowing clusters at 0.5 ppt in the low-base-concentration regime. This low contribution of OOM to the SA–TMA–OOM clusters is caused by the strong SA–TMA interaction, as well as the bulky TMA molecule hindering the binding of the OOM to the cluster due to steric hindrance. MBTCA is the most branched – and therefore also the most rigid – of the three OOMs studied, amplifying the effect of the steric hindrance.

4 Conclusions

Based on previous results using the cluster-of-functional-groups approach, we were able to identify three different oxygenated organic molecules (OOMs), specifically 3-methyl-1,2,3-butanecarboxylic acid (MBTCA), carboxyheptanoic acid (CHA), and pinyl diaterpenylic ester (PDPE), that we believe are able to form thermodynamically stable clusters in the atmosphere with sulfuric acid (SA) and nitrogen-containing bases due to each of them containing three carboxyl groups. Using quantum chemical calculations, we studied the intermolecular interactions between the OOMs and SA–base clusters and used the thermochemical data to study the cluster formation potentials of SA–base–OOM clusters. While the thermochemistry reveals distinct differences among PDPE, CHA, and MBTCA, the cluster formation potentials show similar trends cross all three OOMs, especially in the SA–DMA–OOM systems. We find that all three OOMs are present in the most important outgrowing clusters. The insensitivity to the specific OOM suggests that the functional groups, rather than the organic molecule itself, are most important for cluster formation and growth.

The purely organic dimers were found to be relatively stable. The (PDPE)₂ cluster was 0.7 kcal mol⁻¹ more stable than the (SA)₁(DMA)₁ cluster, while the (CHA)₂ cluster was only 1.2 kcal mol⁻¹ less stable than (SA)₁(DMA)₁. The introduction of SA and bases lowered the binding Gibbs free energy considerably for all three OOMs, with a decrease of roughly 10−25 kcal mol⁻¹. Given that these mixed organic and inorganic clusters are larger than the purely organic clusters studied, extending the study up to (OOM)₃₋₄ may reveal that organics can form stable clusters without the need for inorganic acids and bases.

It would also be interesting to study mixed OOM systems further as these could introduce additional geometric flexibility and potentially expose additional COOH functional groups that facilitate further growth.