Rate enhancement in collisions of sulfuric acid molecules due to long-range intermolecular forces

Collisions of molecules and clusters play a key role in determining the rate of atmospheric new particle formation and growth. Traditionally the statistics of these collisions are taken from kinetic gas theory assuming spherical noninteracting particles, which may significantly underestimate the collision coefficients for most atmospherically relevant molecules. Such systematic errors in predicted new particle formation rates will also affect large-scale climate models. We studied the statistics of collisions of sulfuric acid molecules in a vacuum using atomistic molecular dynamics simulations. We found that the effective collision cross section of the H₂SO₄ molecule, as described by an optimized potentials for liquid simulation (OPLS). OPLS all-atom force field, is significantly larger than the hard-sphere diameter assigned to the molecule based on the liquid density of sulfuric acid. As a consequence, the actual collision coefficient is enhanced by a factor of 2.2 at 300 K compared with kinetic gas theory. This enhancement factor obtained from atomistic simulation is consistent with the discrepancy observed between experimental formation rates of clusters containing sulfuric acid and calculated formation rates using hard-sphere kinetics. We find reasonable agreement with an enhancement factor calculated from the Langevin model of capture, based on the attractive part of the atomistic intermolecular potential of mean force.

1 Introduction

New particle formation from condensable vapours provides an important contribution to the composition of aerosols in the atmosphere which affects air quality as well as the Earth's climate. The positive and negative contributions of atmospheric aerosols to the planet's radiative balance are still not fully understood, and currently constitute one of the largest uncertainties in climate modelling. The earliest stage of new particle formation involves the collisions of individual molecules leading to the appearance of a new molecular complex. In many theoretical approaches, the statistics of such collisions are simply taken from kinetic gas theory, i.e. the molecules are considered to be noninteracting hard spheres, and a collision occurs when the impact parameter, i.e. the perpendicular distance between the spheres' trajectories, is smaller than the sum of the hard spheres' radii. The hard-sphere collision cross section is independent of the relative velocity of the colliding bodies, and the collision rate coefficient for hard spheres of identical radii is customarily expressed as

$$\begin{array}{}\text{(1)}& {\mathit{\beta }}_{\mathrm{HS}}=\sqrt{{\displaystyle \frac{\mathrm{8}{k}_{\mathrm{B}}T}{\mathit{\pi }\mathit{\mu }}}}\mathit{\pi }(\mathrm{2}R{)}^{\mathrm{2}},\end{array}$$

where k_B is the Boltzmann constant, T is the temperature, μ is the reduced mass, and R is the radius of the spheres.

It is well known that acid–base clusters, in particular clusters containing sulfuric acid and ammonia, or amines, are very relevant in nucleation and growth of particles that can serve as cloud condensation nuclei (Almeida et al., 2013). However, such molecules are not necessarily spherical and, despite being charge neutral, exhibit long-range attraction due to interactions between permanent dipoles, permanent and induced dipoles, or induced dipoles (Israelachvili, 2011). Therefore, it is reasonable to expect that particle growth rates or cluster size distributions predicted using collision coefficients from kinetic gas theory will have a systematic error, which needs to be accounted for. In fact, systematic discrepancies have been found between experimental particle formation rates and values predicted from kinetic modelling and cluster dynamics simulations, where hard-sphere collisions are assumed. Kürten et al. (2014) measured the kinetic formation rate of sulfuric acid dimers and found that an enhancement factor of 2.3 needed to be applied to the formation rate obtained from a kinetic model. Lehtipalo et al. (2016) and Kürten et al. (2018) have studied particle formation rates in systems containing sulfuric acid, dimethylamine, and water and concluded that a respective enhancement factor of 2.7 and 2.3 was needed to match experimental particle formation rates.

The effect of long-range interactions between neutral polar molecules on the capture rate constant has been studied by classical trajectory integration (Maergoiz et al., 1996c). The interaction potential between the colliding parties has been approximated by two terms: first, an anisotropic interaction between permanent dipoles proportional to r⁻³, where r is the distance between the centres of mass of the molecules. Second, an isotropic term due to the interaction between permanent dipoles and induced dipoles, and the interaction between induced dipoles, proportional to r⁻⁶. However, such an approximation is inaccurate when the distance between the colliding particles is comparable to their size. Rate coefficients for ion–molecule capture processes have also been studied theoretically in both classic and quantum regimes (Moran and Hamill, 1963; Su and Bowers, 1973; Su et al., 1978; Chesnavich et al., 1980; Clary, 1985; Troe, 1987) or by using trajectory calculations (Dugan and Magee, 1967; Chesnavich et al., 1980; Su and Chesnavich, 1982; Maergoiz et al., 1996a, b). Atomistic simulations have been used to study collisions of Lennard-Jones clusters and atmospherically relevant molecules, but these studies did not analyse or report thermal collision rate coefficients (Napari et al., 2004; Loukonen et al., 2014). Recently, Yang et al. (2018) studied the condensation rate coefficients for Au and Mg clusters at various gas temperatures using molecular dynamics (MD) calculations. The influence of van der Waals forces on the collision rate has also been considered in Brownian coagulation models of ultrafine aerosol particles (Marlow, 1980; Sceats, 1986, 1989).

In the present work, we use atomistic molecular dynamics simulations to study the statistics of collisions between sulfuric acid molecules in a vacuum, determine the collision rate coefficient, and calculate the enhancement factor over kinetic gas theory. Here, we focus on “reactive” collisions, defined by the formation of one or more hydrogen bonds between the molecules. Detailed modelling of e.g. proton transfer processes related to hydrogen bond formation in such reactive collisions would require first-principle simulations (Loukonen et al., 2014); however, the need to simulate a large number of individual trajectories to cover a representative range of impact parameters and relative velocities makes this impossible. In the present study we model the collision rate enhancement due to long-range interactions, which can be decently described by empirical force fields.

In Sect. 2 we discuss technical matters regarding the choice of the force field and the simulation set-up and give a brief overview of the theoretical background of collisions of atmospheric particles. In Sect. 3, simulation results are presented, discussed within the theoretical framework, and compared to analytical and experimental results.

2 Simulation details and theoretical models

2.1 Force field benchmark

We considered two force fields to describe the sulfuric acid molecules in the present study. The first choice was the force field by Ding et al. (2003), fitted specifically to reproduce DFT structures and energies of small clusters of sulfuric acid, bisulfate, and water, in a vacuum. The second choice was the force field by Loukonen et al. (2010), who fitted interaction parameters for sulfuric acid, bisulfate, and dimethylammonium according to the OPLS all-atom procedure (Jorgensen et al., 1996). Both force fields are fitted to reproduce the C₂ geometry of the isolated H₂SO₄ molecule in a vacuum, and the atoms' partial charges create dipole moments of 3.52 and 3.07 debyes, for Ding et al. (2003) and Loukonen et al. (2010), respectively, which is in agreement with experiments (2.7–3.0 debyes) and ab initio calculations (2.7–3.5 debyes) (Sedo et al., 2008). In both force fields intermolecular interactions are described by the sum of Lennard-Jones potentials between atoms i and j separated by a distance r_ij, with distance and energy parameters σ_ij and ε_ij, and Coulomb interactions between the partial charges q_i and q_j,

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{U}_{\mathrm{inter}}& =\sum _{i=\mathrm{1}}^{N_{\mathrm{1}}}\sum _{j=\mathrm{1}}^{N_{\mathrm{2}}}\mathrm{4}{\mathit{\epsilon }}_{ij}\left[{\left({\displaystyle \frac{{\mathit{\sigma }}_{ij}}{r_{ij}}}\right)}^{\mathrm{12}}-{\left({\displaystyle \frac{{\mathit{\sigma }}_{ij}}{r_{ij}}}\right)}^{\mathrm{6}}\right]\\ & +\sum _{i=\mathrm{1}}^{N_{\mathrm{1}}}\sum _{j=\mathrm{1}}^{N_{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }{\mathit{ϵ}}_{\mathrm{0}}}}{\displaystyle \frac{q_i{q}_j}{r_{ij}}}.\end{array}\end{array}$$

However, in the force field by Ding et al. (2003), the geometry of the individual molecule is simply constrained by harmonic potentials with force constants k_ij between all pairs of atoms,

$$\begin{array}{}\text{(3)}& U_{\mathrm{intra}}^{\mathrm{Ding}}=\sum _{i=\mathrm{1}}^{N_{\mathrm{1}}-\mathrm{1}}\sum _{j=i+\mathrm{1}}^{N_{\mathrm{1}}}{\displaystyle \frac{k_{ij}}{\mathrm{2}}}{\left(r_{ij}-r_{ij}^{\mathrm{0}}\right)}^{\mathrm{2}},\end{array}$$

whereas in OPLS the intramolecular interactions consist of the usual sum of two, three, and four-body potentials, i.e. harmonic bonds between covalently bonded atoms, harmonic angles θ between atoms separated by two covalent bonds, and torsion (dihedral angles ϕ) between atoms separated by three covalent bonds,

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{U}_{\mathrm{intra}}^{\mathrm{OPLS}}& =\sum _{i=\mathrm{1}}^{N_{\mathrm{bonds}}}{\displaystyle \frac{k_i^b}{\mathrm{2}}}{\left(r_i-r_i^{\mathrm{0}}\right)}^{\mathrm{2}}+\sum _{j=\mathrm{1}}^{N_{\mathrm{angles}}}{\displaystyle \frac{k_j^{\mathit{\theta }}}{\mathrm{2}}}{\left({\mathit{\theta }}_j-{\mathit{\theta }}_j^{\mathrm{0}}\right)}^{\mathrm{2}}\\ & +\sum _{k=\mathrm{1}}^{N_{\mathrm{dihedrals}}}\sum _{n=\mathrm{1}}^{\mathrm{4}}{\displaystyle \frac{V_n}{\mathrm{2}}}\left[\mathrm{1}+\cos (n{\mathit{\varphi }}^k-{\mathit{\varphi }}_n^k)\right].\end{array}\end{array}$$

To validate the force fields, we compare the structures and energies of four stable configurations of the sulfuric acid dimer illustrated in Fig. 1a–d to ab initio structures and energies at the RI-MP2/CBS//6-31+G* level of theory calculated by Temelso et al. (2012). The results of the benchmark are summarized in Table 1. The force field by Ding et al. (2003) correctly predicts the lowest energy for dimer structure “a”, and the relative energy differences ΔΔE between optimized structures are closer to those obtained in the ab initio calculation than for the OPLS force field; the OPLS force field assigns the lowest energy to structure “d”, which is the highest energy structure in the ab initio calculation. The geometries of the structures agree well with the ab initio result for both force fields, with the OPLS force field reproducing the ab initio hydrogen bond lengths slightly better than the force field by Ding et al. (2003). The binding energies at T=0 K are slightly lower for the force fields (−0.64 and −0.67 eV for OPLS and Ding et al., 2003, respectively) compared with the ab initio value of −0.72 eV. Overall, the force field by Ding et al. (2003) performs slightly better in terms of energetics.

Figure 1 Four minimum energy structures for the sulfuric acid dimer (a–d) used to benchmark the force fields by Ding et al. (2003) and Loukonen et al. (2010) against ab initio calculations by Temelso et al. (2012). Sulfur atoms are yellow, oxygen atoms are red, and hydrogen atoms are white. Hydrogen bonds are indicated by dotted lines and enumerated according to Table 1.

Ding et al. (2003)Ding et al. (2003)

Table 1 Relative energies ΔΔE (eV) and hydrogen bond distances d_O⋯H (Å) for the sulfuric acid dimer structures (a)–(d) in Fig. 1 obtained from ab initio calculations by Temelso et al. (2012) and with the force fields by Ding et al. (2003) and Loukonen et al. (2010), following the OPLS doctrine.

Energy unit conversion: 1 eV ≈96.49 kJ mol${}^{-\mathrm{1}}\approx \mathrm{23.06}$ kcal mol${}^{-\mathrm{1}}\approx \mathrm{38.68}{k}_{\mathrm{B}}T$ at
T=300 K.

However, the vibrational spectra, calculated from the Fourier transform of the velocity autocorrelation functions of an isolated H₂SO₄ molecule in a vacuum exhibit strong differences: the force field by Loukonen et al. (2010) is able to reproduce the experimental and ab initio spectra very well (Hintze et al., 2003; Chackalackal and Stafford, 1966; Miller et al., 2005), whereas the force field by Ding et al. (2003) is not, as shown in Fig. 2. Intramolecular vibrations are relevant in the context of molecular collisions, e.g. when studying energy transfer between different internal degrees of freedom during and after the collision. Also, the OPLS all-atom procedure allows for transferable potentials, as opposed to the Ding et al. (2003) force field which cannot easily be extended to other chemical compounds in future studies. For these two reasons, we decided to use the OPLS force field by Loukonen et al. (2010) for the collision simulations.

Figure 2 Vibrational spectra of the sulfuric acid molecule obtained with OPLS (Loukonen et al., 2010) (blue area) and the force field by Ding et al. (2003) (orange area). The spectra of Ding et al. (2003) has its highest frequency mode at 6565 cm⁻¹ outside of the range of the figure. The positions of the peaks in the experimental spectra (Hintze et al., 2003; Chackalackal and Stafford, 1966; Miller et al., 2005) and ab initio calculations (Miller et al., 2005) are indicated by dashed and solid grey lines, respectively.

2.2 Potential of mean force of two sulfuric acid molecules

We first calculated the binding free energy of two sulfuric acid molecules in a vacuum as described by the force fields of Loukonen et al. (2010) and Ding et al. (2003) The potential of mean force (PMF) as a function of the sulfur–sulfur distance was calculated from a well-tempered metadynamics simulation (Barducci et al., 2008), using the plumed plug-in (Tribello et al., 2014) for lammps (Plimpton, 1995). We used a velocity Verlet integrator with a time step of 1 fs, to correctly resolve the motion of the hydrogen atoms. The Lennard-Jones interactions were cut off at 14 Å and electrostatic interactions were only evaluated in direct space, with a cut-off at 40 Å. We employed 24 random walkers, and Gaussians with a width of 0.1 Å and initial height of k_BT were deposited every 500 steps along the collective variable; a harmonic wall was used to restrict the collective variable to values below 35 Å. A bias factor of 5 was chosen, and a Nosé–Hoover chain thermostat of length 5 with a time constant of 0.1 ps was used to maintain a temperature of T=300 K. The combined length of the trajectories was 120 ns for each force field. Both PMFs, shown in Fig. 3, exhibit a minimum at r=4.1 Å, and the binding free energies are $\mathrm{\Delta }F=-\mathrm{0.29}$ and −0.27 eV for Loukonen et al. (2010), and Ding et al. (2003), respectively. This is in excellent agreement with the ab initio value of $\mathrm{\Delta }G=-\mathrm{0.30}$ eV obtained from Boltzmann-averaging over the four minimum energy dimer structures by Temelso et al. (2012) at 298.15 K, and more recent calculations at a higher level of theory, which predict slightly weaker binding (−0.23 to −0.26 eV) (Elm et al., 2016; Myllys et al., 2017).

Figure 3 Potential of mean force between two H₂SO₄ molecules as a function of the sulfur–sulfur distance calculated by metadynamics simulation for the OPLS force field (Loukonen et al., 2010) (solid black line) and the force field by Ding et al. (2003) (dashed black line) at 300 K. For comparison, the Keesom (thermally averaged permanent dipole–permanent dipole) interaction between two point dipoles (see Eq. A11) of 3.07 debyes is depicted by the dotted line. The dashed red line shows the attractive potential (Eq. 8) based on the calculated PMF curve of the OPLS force field.

2.3 Collision simulation set-up

Molecular dynamics simulations were performed with the lammps code, using a velocity Verlet integrator with a time step of 1 fs. The Lennard-Jones interactions were cut off at 14 Å and electrostatic interactions were only evaluated in direct space, with a cut-off at 120 Å. The simulations were carried out in the NVE ensemble, as the colliding molecules constitute a closed system (under atmospheric conditions collisions with the carrier gas are rare on the timescale of collisions between sulfuric acid molecules). In order to determine the molecules' collision probability as a function of impact parameter and relative velocity, the following set-up was used: first, two sulfuric acid molecules were placed in the simulation box, separated by 100 Å along x and the impact parameter b was set along the z direction, ranging from 0 to 17.5 Å in steps of 0.5 Å. Atomic velocities were randomly assigned from a Maxwell–Boltzmann distribution at T=300 K, and the centre of mass motion of each molecule was removed separately. Then the system was evolved for 50 ps, to randomize the intermolecular orientation and ensure equipartition of energy along the intramolecular degrees of freedom. At t=50 ps, each molecule received a translational velocity along the x direction, $v_x=\pm v/\mathrm{2}$, where v denotes the relative velocity, to set them on a potential collision course. The simulation was continued for another 250 ps, to ensure the possibility of a collision, even at the smallest relative velocities. For a collision of two identical molecules with molecular mass m, the relative velocities follow the Maxwell–Boltzmann distribution with reduced mass $\mathit{\mu }=m/\mathrm{2}$. We sampled relative velocities between 50 and 800 m s⁻¹, in steps of 50 m s⁻¹. A total of 99 % of the distribution lies within this range at T=300 K. For each value of the impact parameter b and the relative velocity v, 1000 simulations were carried out starting with different initial atomistic velocities, to ensure good sampling. The simulation set-up is very similar to the one recently used by Yang et al. (2018). Sulfuric acid molecules bind to each other via the formation of one or more hydrogen bonds. However, even if a collision course leads to an attachment of the two molecules, a portion of the kinetic energy will be redistributed on the degrees of freedom of the formed complex, and this excess energy can lead to a rapid dissociation in the absence of a thermalizing medium. To automate the analysis of over half a million individual trajectories, we define a collision as a trajectory during which the electrostatic energy (E_Coul) is lower than a threshold value of −0.25 eV in at least 10 consequent frames (100 fs), indicative of the formation of one or more hydrogen bonds. Three examples of simulated trajectories with a relative velocity closest to the mean velocity at 300 K (350 m s⁻¹) and an impact parameter of 8.5 Å are illustrated in Fig. 4.

Figure 4 Three example MD trajectories with a relative velocity of 350 m s⁻¹ and an impact parameter of 8.5 Å after 50 ps of energy relaxation and the first 25 ps of collision simulation. Panel (a) shows both the distance between the two sulfur atoms (solid lines) and the centre-of-mass distance of the molecules (dotted lines) during the simulations. The electrostatic energies E_Coul, which is used to determine the possible collision event, for successful (1), unsuccessful (2), and successful but instantly evaporating (3, continuously 0.18 ps below the threshold) trajectories are shown in panel (b) with the threshold value (dashed black line).

The results from the atomistic simulations will be compared to different theoretical models described in the following.

2.4 Classical model of capture in a field of force

As the collision rate in the context of atomistic simulations is defined as the reaction rate of hydrogen bonding, the related theoretical models are often based on the assumption that if the trajectory of the colliding molecules is able to surmount a centrifugal barrier the reaction is certain. This is known as the capture approximation; to emphasize this conceptual difference between simulations and theoretical models, we use the word capture instead of collision to refer to theory-based results.

The interaction between two identical polar molecules is usually written as

$$\begin{array}{}\text{(5)}& V={\displaystyle \frac{({\mathit{d}}_{\mathrm{1}}\cdot {\mathit{d}}_{\mathrm{2}})-\mathrm{3}({\mathit{d}}_{\mathrm{1}}\cdot \mathit{n})({\mathit{d}}_{\mathrm{2}}\cdot \mathit{n})}{r^{\mathrm{3}}}}+U\left(r\right),\end{array}$$

where d₁ and d₂ are the dipole moment vectors of the molecules, n is the unit vector along the distance vector r connecting the centres of mass of the molecules, and U(r) is a spherically symmetric potential, usually proportional to r⁻⁶. The capture rate constant for such a potential can only be calculated numerically (Maergoiz et al., 1996c). In the present section we consider only the isotropic part U(r) of the interaction described by Eq. (5), the effect of the anisotropic part (first term in Eq. 5) is discussed in Appendix A. Then, the Langevin model of capture (Langevin, 1905) can be used to calculate the critical impact parameter beyond which point-like colliding particles in a vacuum will escape from each other. Here, the motion of the two colliding molecules is reduced to a one-body problem in an external central field by using an effective potential containing dispersion and centrifugal terms (Landau and Lifshitz, 1976):

$$\begin{array}{}\text{(6)}& U_{\mathrm{eff}}\left(r\right)=U\left(r\right)+{\displaystyle \frac{L^{\mathrm{2}}}{\mathrm{2}\mathit{\mu }{r}^{\mathrm{2}}}},\end{array}$$

where r is the distance of the colliding body from the centre of the field, and L is the angular momentum. Both the total energy of the system (which equals the initial translational energy μv²∕2 at r→∞) and the angular momentum are conserved. The centrifugal term introduces an energy barrier, and for a successful capture at the barrier (r=r_max) the translational energy $\mathit{\mu }{\dot{r}}_{\max }^{\mathrm{2}}/\mathrm{2}$ has to be positive. As the angular momentum equals L=μvb, the condition for b² to ensure a capture is

$$\begin{array}{}\text{(7)}& b^{\mathrm{2}}<r_{\max }^{\mathrm{2}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{2}U\left(r_{\max }\right)}{\mathit{\mu }{v}^{\mathrm{2}}}}\right).\end{array}$$

In case of a simple attractive potential (repulsive forces can be neglected, as the studied velocities are relatively low),

$$\begin{array}{}\text{(8)}& U\left(r\right)=-\mathit{ϵ}{\left({\displaystyle \frac{r_{\mathrm{0}}}{r}}\right)}^{\mathrm{6}},\end{array}$$

the square of the critical impact parameter can be written as

$$\begin{array}{}\text{(9)}& b_{\mathrm{c}}^{\mathrm{2}}={\left({\displaystyle \frac{\mathrm{27}\mathit{ϵ}}{\mathrm{2}\mathit{\mu }{v}^{\mathrm{2}}}}\right)}^{\mathrm{1}/\mathrm{3}}{r}_{\mathrm{0}}^{\mathrm{2}}.\end{array}$$

It is preferable to consider the squared value of b, as the capture cross section is calculated as ${\mathit{\sigma }}_{\mathrm{c}}=\mathit{\pi }{b}_{\mathrm{c}}^{\mathrm{2}}$. It is important to note that in the Langevin model, the total energy is divided strictly to the translational and potential energy, the internal degrees of freedom of the two bodies are considered to be completely decoupled, i.e. exchange of translational energy to rotations and vibrations that will occur in a real molecule is completely neglected.

2.5 Brownian model of aerosol coagulation

In the study by Kürten et al. (2014), a model of Brownian coagulation in a field of force (Sceats, 1986, 1989) was used to estimate the collision enhancement factor for neutral cluster formation involving sulfuric acid in a free molecule regime (Chan and Mozurkewich, 2001). The model is based on solving the Fokker–Planck equation for a pair of Brownian particles whose motion is determined by a thermal random force (Sceats, 1986). In the paper by Chan and Mozurkewich (2001), the Hamaker constant describing the strength of the van der Waals interaction was fitted to experiments with uncharged H₂SO₄∕H₂O particles with diameters of 49–127 nm, yielding a collision enhancement factor value of 2.3 at 300 K. Although the Hamaker constant is usually considered to be size-independent, there may be enhanced interaction for very small particle sizes, with radii of the order of 1 nm (Pinchuk, 2012).

For the attractive potential described by Eq. (8), the collision enhancement factor over the kinetic gas theory rate, $W_{\mathrm{B}}={\mathit{\beta }}_{\mathrm{B}}/{\mathit{\beta }}_{\mathrm{HS}}$, from the Brownian coagulation model in the free molecule limit can be written as (Sceats, 1989)

$$\begin{array}{}\text{(10)}& W_{\mathrm{B}}={\left({\displaystyle \frac{\mathrm{3}\mathit{ϵ}}{k_{\mathrm{B}}T}}\right)}^{\mathrm{1}/\mathrm{3}}{\left({\displaystyle \frac{r_{\mathrm{0}}}{\mathrm{2}R}}\right)}^{\mathrm{2}}\exp \left({\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\right).\end{array}$$

To compare the Brownian and the Langevin models to atomistic simulation results using the OPLS force field, we have determined the attractive potential parameters ϵ and r₀ to be equal to −2ΔF and the sulfur-sulfur distance at which the PMF reaches its minimum, respectively.

It should be noted that the model of Brownian coagulation does not describe the correct transport physics of collisions of molecules in the gas phase. For a discussion on the transition from the free molecular (ballistic) regime to the continuum (diffusive) regime, see e.g. Ouyang et al. (2012).

3 Results and discussion

The statistics of the collision probabilities as a function of the impact parameter and relative velocity, P(b,v), obtained from the atomistic simulations are shown as a heat map in Fig. 5 where white indicates a certain collision event (defined by the formation of one or more hydrogen bonds) and black indicates zero collision probability. The critical impact parameter (above which a collision unlikely occurs) has a negative exponential dependency on the relative velocity, which corresponds well with the Langevin model. A more detailed plot is provided in Fig. B1, where the sigmoidal probability curves are shown for each velocity separately. Predominantly, the collision probability approaches unity with a lower relative velocity and a shorter impact parameter. However, the collision probability decreases slightly at high relative velocities where rapid re-dissociation of the complex can be caused by high kinetic energy and slow redistribution of the energy to vibrational modes of the formed cluster. At the slowest velocity of 50 m s⁻¹, and small values of the impact parameter, the collision probability is also reduced. This happens because in some cases the fluctuations in the intermolecular energy, just as they come within interaction range, are sufficient to exceed the very small initial translational energies of the colliding molecules, effectively repelling them (see Appendix C for a more detailed discussion).

Figure 5 Heat map of the collision probability of sulfuric acid molecules plotted as a function of impact parameter b and relative velocity v obtained from molecular dynamics simulation. The impact parameter equivalent to the hard-sphere collision area (b=2R) and the critical impact parameter obtained from the Langevin capture model (Eq. 9) are indicated by the dashed red and the solid white lines, respectively.

The dynamical collision cross section, obtained from the integral over the collision probability functions,

$$\begin{array}{}\text{(11)}& {\mathit{\sigma }}_{\mathrm{d}}\left(v\right)=\mathit{\pi }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}{b}^{\mathrm{2}}P(b,v),\end{array}$$

is consistently decreasing with relative velocity v. Even though it can be seen in Fig. 5, and especially in Fig. B1, that values of σ_d are smaller than the corresponding Langevin capture cross sections σ_c, the velocity dependence of the change in σ_d is in very close agreement with the Langevin model solution

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{c}}\left(v\right)}{{\mathit{\sigma }}_{\mathrm{c}}\left(v_{\mathrm{0}}\right)}}={\left({\displaystyle \frac{v}{v_{\mathrm{0}}}}\right)}^{-\mathrm{2}/\mathrm{3}},\end{array}$$

where v₀ is a reference velocity, as shown in Fig. 6. The importance of contributions from long-range interactions to the collision cross section is evident, as σ_c is proportional to $v^{-\mathrm{4}/n}$, for interactions decaying with r⁻ⁿ. Furthermore, as can be seen in Fig. B1, the critical impact parameters from the Langevin model match rather well with the tails of the simulated collision probability curves, the intersection is located without exception at $P(b,v)\approx \mathrm{0.2}$.

Figure 6 Ratio between the collision cross section σ and a reference value σ₀ as a function of relative velocity. Here the reference value σ₀ is the cross section corresponding to a relative velocity v₀=350 m s⁻¹. The red crosses show the ratio of dynamical cross sections obtained from MD simulations, and the solid black line shows the relation predicted by the Langevin model in Eq. (12).

The discrepancy between σ_d and σ_c is the result of the assumptions made in the Langevin model, where the capture is considered to be orientation-independent and the particles do not have any internal structure. If the anisotropy of the dipole–dipole potential is taken into account, as in Eq. (5), the capture cross section will be reduced (this has been estimated in Appendix A using a numerical approach provided by Maergoiz et al., 1996c). However, if two molecules are able to move rather close to each other, translational energy can be transferred to rotational and vibrational modes; therefore, the motion over the centrifugal barrier is hindered, and the critical impact parameter is effectively reduced. Additionally, steric hindrance caused by intermolecular orientations incompatible with the formation of hydrogen bonds will also lower the collision probability. Due to coupling, steric hindrance, and other dynamical effects, the ratio between the cross sections σ_d and σ_c is 0.82 on average for the collision of two sulfuric acid molecules.

The canonical collision rate coefficient can be calculated in a similar fashion to Eq. (1), but as the collision probabilities P(b,v) obtained from the atomistic simulations depend on both the velocity and the impact parameter, the MD-based collision rate coefficient is calculated by integrating over both the relative velocity distribution f(v) and b² as

$$\begin{array}{}\text{(13)}& {\mathit{\beta }}_{\mathrm{MD}}=\mathit{\pi }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}v\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}{b}^{\mathrm{2}}vf\left(v\right)P(b,v).\end{array}$$

While the thermal velocity distribution f(v) of the colliding molecules can be altered freely to correspond with an arbitrary temperature, the effect of the internal motion to the collision probability function is not necessarily temperature-invariant. However, in Appendix B it is shown that a moderate change (simulations carried out at 250 and 400 K) in the internal kinetic energy does not effect the collision probabilities significantly. Therefore, we used the collision probability distributions calculated at 300 K to compute the collision rate coefficients for the atmospherically relevant temperature range T=225–425 K (see Fig. 7).

Figure 7 Collision rate coefficient β (a) and the enhancement factor W (b) as a function of temperature calculated for the hard-sphere, MD, Langevin, and Brownian approaches.

For the Langevin model (Eq. 9) the expression for the canonical capture rate coefficient can be simplified to

$$\begin{array}{}\text{(14)}& {\mathit{\beta }}_{\mathrm{L}}=\mathit{\pi }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}vvf\left(v\right)b_{\mathrm{c}}^{\mathrm{2}}\left(v\right).\end{array}$$

As the potential of mean force is required for the Langevin model, the coefficients are only calculated at 250, 300, and 400 K (see Appendix B for further details).

In addition to the coefficients obtained by different approaches, we examine the enhancement factor W relative to the kinetic gas theory rate expressed in Eq. (1), where a hard-sphere radius R=2.77 Å was calculated from the bulk liquid density of sulfuric acid, ρ=1830 kg m⁻³, assuming a volume fraction of one. Thus, after performing the integration, the enhancement factor obtained using the Langevin model can be expressed analytically as

$$\begin{array}{}\text{(15)}& W_{\mathrm{L}}={\displaystyle \frac{{\mathit{\beta }}_{\mathrm{L}}}{{\mathit{\beta }}_{\mathrm{HS}}}}=\mathrm{\Gamma }\left({\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}\right){\left({\displaystyle \frac{\mathrm{2}\mathit{ϵ}}{k_{\mathrm{B}}T}}\right)}^{\mathrm{1}/\mathrm{3}}{\left({\displaystyle \frac{r_{\mathrm{0}}}{\mathrm{2}R}}\right)}^{\mathrm{2}},\end{array}$$

where Γ(x) denotes the Gamma function, and ϵ and r₀ are the parameters based on the potential of mean force between two sulfuric acid molecules in a vacuum, using the functional form in Eq. (8). The enhancement factor from the Brownian model of aerosol coagulation, W_B, was calculated analytically from Eq. (10).

We find that W_MD≈2.20, W_L≈2.59, and W_B≈3.06 at 300 K. Both the collision rate coefficients and the enhancement factors are shown in Fig. 7. As the hard-sphere collision rate increases linearly as the molecular velocity is proportional to $\sqrt{T}$, the rate coefficient based on atomistic simulations only slightly increases with temperature due to the exponential narrowing of the collision cross section as the relative velocity increases. However, the dependency is different in case of the Langevin model as the rate coefficient decreases with increasing temperature. This is due to the neglect of the anisotropy of the dipole–dipole interactions in the Langevin model: at higher temperatures, the effect of anisotropy becomes less important; therefore, the model overestimates the collision rate less, compared with lower temperatures. The same effect can be observed in the Brownian coagulation model. Indeed, as seen in Fig. 7, the enhancement factors based on the theoretical models reach better agreement with the factor obtained from atomistic simulations as the temperature rises. As a result of the above-mentioned reasons, the enhancement factor over the hard-sphere collision rate coefficient is lower at higher temperatures regardless of the approach chosen.

The enhancement factor obtained by atomistic simulations is in very good agreement with the kinetic modelling on the recent experimental results of the formation of atmospheric sulfuric acid dimers (Kürten et al., 2014) and small clusters of sulfuric acid, dimethylamine, and water (Lehtipalo et al., 2016; Kürten et al., 2018). In these studies, the enhancement factor was estimated to be 2.3–2.7 using the Brownian coagulation model and van der Waals interactions fitted to experimental data, as mentioned earlier, whereas according to our molecular model of long-range interaction the basic Brownian model using the attractive part of the potential of mean force overestimates the rate enhancement factor by about 40 %.

4 Conclusions

In summary, we have benchmarked two classical force fields against experimental and ab initio data and determined that the OPLS force field by Loukonen et al. (2010) was able to describe the geometry and vibrational spectra of the isolated sulfuric acid molecule, as well as the geometry and binding free energy of the sulfuric acid dimer. We studied the statistics of collisions of sulfuric acid molecules in a vacuum using molecular dynamics simulations and compared our results against simple theoretical models. We found that the effective collision cross section of two H₂SO₄ molecules, as described by the OPLS force field, is significantly larger than the hard-sphere diameter assigned to the molecule based on the liquid density of sulfuric acid. As a consequence, we find that the collision coefficient for sulfuric acid molecules is enhanced by a factor of 2.2, compared with kinetic gas theory at 300 K. This enhancement factor obtained from atomistic simulation is consistent with the discrepancy observed between experimental formation rates of clusters containing sulfuric acid and rates calculated using hard-sphere kinetics. At a temperature range from 250 to 400 K, the rate enhancement factor is monotonously decreasing with increasing temperature, however the drop is less than 20 %. The velocity dependence of the simulated dynamical collision cross section is in good agreement with the Langevin model solution. We also note that the enhancement factor obtained from the Langevin model using the attractive part of the intermolecular potential is slightly overestimated due to the imperfect treatment of the dipole–dipole interaction; nevertheless, in the atmospherically relevant temperature range the factor is within 30 % of the result from the atomistic simulation, at a fraction of the computational cost.

In the future, the atomistic collision modelling approach presented in this work can be applied to other atmospherically relevant molecules, clusters, or ions, exhibiting dipoles of varying magnitudes – and in some cases several times larger than that of the sulfuric acid molecule – to help understand the effect of long-range interactions in cluster formation rates. However, before we can quantitatively assess the influence of collision rate enhancement on atmospherical new particle formation rates obtained from cluster dynamics models (for example, Atmospheric Cluster Dynamics Code; McGrath et al., 2012), it is necessary to obtain the enhancement factors for all of the relevant collisions between clusters of different sizes and composition, as the pathway for growth may change – a formidable task, even if only the simplest acid–base clusters were considered. Thus, future work should also be aimed at finding simple models for predicting approximate rate enhancements, based on just a few physico-chemical properties, such as molecular structures, dipole moments, or charge distributions, of the interacting molecules and/or clusters.

Appendix A: Effect of anisotropy on the dipole–dipole capture rate

Maergoiz et al. (1996c), using classical trajectory integration, calculated the capture rate constant when two identical polar molecules interact via a potential containing anisotropic ($\propto r^{-\mathrm{3}}$) and isotropic ($\propto r^{-\mathrm{6}}$) terms:

$$\begin{array}{}\text{(A1)}& V={\displaystyle \frac{({\mathit{d}}_{\mathrm{1}}\cdot {\mathit{d}}_{\mathrm{2}})-\mathrm{3}({\mathit{d}}_{\mathrm{1}}\cdot \mathit{n})({\mathit{d}}_{\mathrm{2}}\cdot \mathit{n})}{r^{\mathrm{3}}}}-{\displaystyle \frac{C}{r^{\mathrm{6}}}},\end{array}$$

where d₁ and d₂ are the dipole moment vectors of the molecules, n is the unit vector along the distance vector r connecting the centres of mass of the molecules, and C is the isotropic interaction constant.

As in Eq. (1), the capture rate coefficient in an anisotropic field is given by

$$\begin{array}{}\text{(A2)}& {\mathit{\beta }}_{\mathrm{aniso}}=\sqrt{{\displaystyle \frac{\mathrm{8}{k}_{\mathrm{B}}T}{\mathit{\pi }\mathit{\mu }}}}{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{aniso}},\end{array}$$

where the thermal capture cross section can be calculated using a fitting function κ(θ,M) as

$$\begin{array}{}\text{(A3)}& {\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{aniso}}=\mathit{\pi }{\left({\displaystyle \frac{d^{\mathrm{2}}}{k_{B}T}}\right)}^{\mathrm{2}/\mathrm{3}}{\mathit{\theta }}^{\mathrm{1}/\mathrm{6}}\mathit{\kappa }(\mathit{\theta },M),\end{array}$$

where d is the molecular dipole moment (for the OPLS model of sulfuric acid, $d=\left|{\mathit{d}}_{\mathrm{1}}\right|=\left|{\mathit{d}}_{\mathrm{2}}\right|=\mathrm{3.07}$ debyes). Maergoiz et al. (1996c) use two dimensionless parameters in their model:

$$\begin{array}{}\text{(A4)}& \mathit{\theta }={\displaystyle \frac{C{k}_{B}T}{d^{\mathrm{4}}}}\end{array}$$

and

$$\begin{array}{}\text{(A5)}& M={\displaystyle \frac{\mathit{\mu }{d}^{\mathrm{4}/\mathrm{3}}}{\mathrm{2}I(k_{B}T{)}^{\mathrm{2}/\mathrm{3}}}}.\end{array}$$

Based on our MD simulations using the OPLS force field, the average moment of inertia I of a vibrating sulfuric acid molecule is 100.04 amuŲ, which deviates slightly from the values of 100.66 and 104.94 amuŲ obtained experimentally (Kuczkowski et al., 1981) and from quantum chemical calculations (Zapadinsky et al., 2019), respectively. The fitting function is obtained from classical trajectory calculations and is expressed as

$$\begin{array}{}\text{(A6)}& \ln \mathit{\kappa }(\mathit{\theta },M)=a_{\mathrm{0}}+{\left({\displaystyle \frac{a_{\mathrm{1}}^{\mathrm{2}}z}{\sinh \left(z\right)}}+{\displaystyle \frac{z^{\mathrm{2}}}{\mathrm{36}}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where

$$\begin{array}{}\text{(A7)}& z=a_{\mathrm{2}}+\ln \mathit{\theta }.\end{array}$$

The reported fitting parameters according to Maergoiz et al. (1996c) are

$$\begin{array}{}\text{(A8a)}& {\displaystyle }{a}_{\mathrm{0}}=\mathrm{0.2406}-\mathrm{0.1596}{\left(\mathrm{1}+\mathrm{1.9192}{M}^{\mathrm{0.9935}}\right)}^{-\mathrm{1}},\text{(A8b)}& {\displaystyle }{a}_{\mathrm{1}}=\mathrm{0.253}-\mathrm{0.04573}{\left(\mathrm{1}+\mathrm{1.1645}{M}^{\mathrm{0.6422}}\right)}^{-\mathrm{1}},\text{(A8c)}& {\displaystyle }{a}_{\mathrm{2}}=\mathrm{1.7617}+\mathrm{0.9577}{\left(\mathrm{1}+\mathrm{1.9192}{M}^{\mathrm{0.9935}}\right)}^{-\mathrm{1}}.\end{array}$$

As all different long-range interactions are included in the attractive part of the potential of mean force between two sulfuric acid molecules (Eq. 8), to exclude the dipole–dipole interaction from the isotropic interaction, the constant C is written as

$$\begin{array}{}\text{(A9)}& C=(\mathrm{1}-f)\mathit{ϵ}{r}_{\mathrm{0}}^{\mathrm{6}},\end{array}$$

where f is a factor denoting the relative magnitude of the anisotropic interaction between permanent dipoles with respect to the total interaction. Thus, the rate coefficient is controlled by the relative dipole–dipole interaction and the enhancement factor over the kinetic gas theory can be written as

$$\begin{array}{}\text{(A10)}& W_{\mathrm{aniso}}\left(f\right)={\displaystyle \frac{{\mathit{\beta }}_{\mathrm{aniso}}\left(f\right)}{{\mathit{\beta }}_{\mathrm{HS}}}}={\displaystyle \frac{{\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_{\mathrm{aniso}}\left(f\right)}{\mathit{\pi }(\mathrm{2}R{)}^{\mathrm{2}}}}.\end{array}$$

Figure A1 shows the rate enhancement as a function of the interaction factor f at 300 K. As the anisotropic part does not contribute, i.e. f=0, the enhancement factor is less than 4 % higher than the value obtained from the Langevin model (the statistical error of the thermal capture cross section is about 2 %; Maergoiz et al., 1996c).

Figure A1 Enhancement factor W_aniso over the kinetic gas theory calculated from the anisotropic dipole–dipole collision cross section as a function of the interaction factor f at 300 K (black solid line). The red cross denotes the case where the dipole–dipole contribution of the total interaction equals the Keesom equation Eq. (A11) (W_aniso=2.41). The enhancement factors obtained from the Langevin model and the MD simulations are shown as the dotted and dashed grey lines, respectively.

As we were unable to distinguish the actual dipole–dipole interaction from the total attractive potential, we estimated the interaction using the Keesom equation (see Fig. 3):

$$\begin{array}{}\text{(A11)}& U_{\mathrm{K}}=-{\displaystyle \frac{d^{\mathrm{4}}}{\mathrm{24}{\mathit{\pi }}^{\mathrm{2}}{\mathit{ϵ}}_{\mathrm{0}}^{\mathrm{2}}{k}_{\mathrm{B}}T{r}^{\mathrm{6}}}}.\end{array}$$

According to Eqs. (8) and (A11), about one-third of the attractive potential is due to dipole–dipole interaction; consequently, the enhancement factor is lower than for the isotropic field. The estimated rate enhancement factors at 250, 300, and 400 K are shown in Table A1. Thus, by taking the anisotropy of the intermolecular potential into account, the estimated capture rate is in better agreement with the result obtained using atomistic simulations.

Table A1 The attractive potential parameters ϵ and r₀ for H₂SO₄–H₂SO₄ interaction based on the PMF calculations with the estimated anisotropic interaction factor $f=U_{\mathrm{K}}/U$ and the corresponding enhancement factors calculated by the Langevin model W_L, the anisotropic approach W_aniso(f), and the atomistic simulations W_MD.

Appendix B: The temperature dependence of collision probabilities and interaction parameters

The canonical collision rate coefficient can be calculated from the collision probabilities obtained from atomistic simulation at arbitrary temperatures by shifting the Maxwell–Boltzmann relative velocity distribution, provided that changes in the internal motion of the molecules do not affect the collision probabilities. We tested the effect of the different rotational and vibrational motion on the collision statistics in an atmospherically relevant temperature range by carrying out MD collision simulations for a subset of impact parameters b where the molecules' atomistic velocities were drawn from Maxwell–Boltzmann distributions corresponding to 250 and 400 K, instead of 300 K. As shown in Fig. B1, such a moderate change in temperature does not affect the collision probabilities between two sulfuric acid molecules.

Figure B1 Collision probabilities of sulfuric acid molecules, as a function of the impact parameter squared, for different values of the relative velocity, obtained from molecular dynamics simulations at 300 K (solid coloured lines), at 250 K (coloured dots), and at 400 K (coloured crosses). The step-like collision probabilities for a hard-sphere model (b²=(2R)²), or those obtained from the Langevin capture model (Eq. 9), are indicated by the solid black, and dashed coloured lines, respectively.

In order to vary temperature in calculating the thermal collision rate coefficient using the Langevin approach, the potential of mean force between two sulfuric acid molecules was calculated at 250, 300, and 400 K, and the parameters describing the attractive intermolecular interaction (Eq. 8) are reported in Table A1.

Appendix C: Intermolecular repulsion at low velocities in MD simulations

As shown in Figs. 5 and B1, for small values of the impact parameter and the initial relative velocity between two colliding molecules in the atomistic simulations, the collision probability can be considerably smaller than unity, which seems counter-intuitive at first. This is due to the fact that the intermolecular interaction is anisotropic and the molecules are rotating, which can lead to instantaneous repulsion even at distances where the intermolecular potential of mean force is slightly attractive. If the initial translational kinetic energy is low enough, the temporary fluctuations in intermolecular energy can alter the translational motion and eventually lead to a definitive separation of the molecules.

Figure C1 Time evolution of the intermolecular energy (red line) and the sulfur–sulfur distance (solid black line) in one MD trajectory with relative velocity v=50 m s⁻¹ and b=0 Å where a collision occurs (a) and one where the molecules are repelled at range (b). The initial translational kinetic energy of the molecules (1.27 meV) is indicated by the dashed black line.

This process is illustrated in Fig. C1, which shows the evolution of the intermolecular energy and distance in one trajectory with b=0 Å and v=50 m s⁻¹ where a collision occurs and a second one where the molecules are repelled at range. While the fluctuation of the intermolecular energy exceeds the initial translational energy of 1.27 meV in both cases, in the trajectory exhibiting a repulsion, the large positive energy fluctuations are longer lived and dominate the interaction.