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 non-interacting particles, which may significantly underesti5 mate 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 have studied the statistics of collisions of sulfuric acid molecules in vacuum by atomistic molecular dynamics 10 simulations. We have found that the effective collision cross section of the H2SO4 molecule, as described by an OPLSAll Atom force field, is significantly larger than the hardsphere diameter assigned to the molecule based on the liquid density of sulfuric acid. As a consequence, the actual col15 lision coefficient is enhanced by a factor 2.2 at 300 K, compared to 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 20 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.


25
New particle formation from condensable vapours gives 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 atmo-spheric aerosols to the planet's radiative balance are still not 30 fully understood, and currently constitute one of the largest uncertainties in climate modelling. The earliest stage of new particle formation involves collisions of individual molecules leading to the appearance of a new molecular complex. In many theoretical approaches, the statistics of such collisions 35 are simply taken from kinetic gas theory, i.e. the molecules are considered to be non-interacting 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 40 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 where k B is Boltzmann's constant, T is the temperature, µ is 45 the reduced mass and R is the radius of the spheres. It is well known that acid-base clusters, and 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). 50 Such molecules, however, are not necessarily spherical and despite being charge neutral, exhibit long-ranged 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 55 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 forma-tion 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 sulphuric 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 sulphuric acid, dimethylamine and water and concluded that an enhancement factor of 2.7 and 2.3, respectively, was needed to match experimental particle for-10 mation 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 15 been approximated by two terms: First, an anisotropic interaction between permanent dipoles, proportional to r −3 , where r is the distance between the centres of mass of the molecules. Second, an isotropic term due to the interaction between permanent dipole and induce dipole, and the inter-20 action between induced dipoles, proportional to r −6 . 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 25 regime (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 Jr. and Magee, 1967;Chesnavich et al., 1980;Su and Chesnavich, 1982;Maergoiz et al., 1996a, b). Atomistic simulations have been 30 used to study collisions of Lennard-Jones clusters and atmospherically relevant molecules, but these studies did not analyze or report thermal collision rate coefficients (Napari et al., 2004;Loukonen et al., 2014). Recently, Yang et al. (2018) have studied the condensation rate coefficients for Au and 35 Mg clusters at various gas temperatures using molecular dynamics calculations. The influence of Van der Waals forces on the collision rate has also been considered in Brownian coagulation models of ultra-fine aerosol particles (Marlow, 1980;Sceats, 1986Sceats, , 1989.

40
In the present work, we use atomistic molecular dynamics simulations to study the statistics of collisions between sulfuric acid molecules in vacuum, determine the collision rate coefficient and calculate the enhancement factor over kinetic gas theory. We are here focusing on "reactive" col-45 lisions, 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 50 large number of individual trajectories to cover a representative range of impact parameters and relative velocities makes this impossible. In the present study we are modelling the collision rate enhancement due to long-range interactions, which can be decently described by empirical force fields.

55
In section 2 we discuss technical matters of the choice of force field and the simulation setup and give a brief overview of the theoretical background for collisions of atmospheric particles. In section 3, simulation results are presented, discussed within the theoretical framework and compared to an-60 alytical and experimental results.
2 Simulation details and theoretical models

Force field benchmark
We have considered two force fields to describe the sulfuric acid molecules in the present study. The first choice was 65 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 vacuum. The second choice was the force field by Loukonen et al. (2010), who had fitted interaction parameters for sulfuric acid, bisulfate 70 and dimethylammonium according to the OPLS-All Atom procedure (Jorgensen et al., 1996). Both force fields are fitted to reproduce the C 2 geometry of the isolated H 2 SO 4 molecule in vacuum, and the atoms' partial charges create dipole moments of 3.52 and 3.07 Debye, for Ding et al. and 75 Loukonen et al., respectively, in agreement with experiments (2.7-3.0 Debye) and ab initio calculations (2.7-3.5 Debye) (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 dis-80 tance and energy parameters σ ij and ε ij , and Coulomb interactions between the partial charges q i and q j , However, in the force field by Ding et al., the geometry of 85 the individual molecule is simply constrained by harmonic potentials with force constants k ij between all pairs of atoms, while in OPLS the intramolecular interactions consist of the 90 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 torsions (dihedral angles φ) between atoms separated by three covalent bonds, To validate the force fields, we compare the structures and 5 energies of four stable configurations of the sulfuric acid dimer illustrated in Fig. 1(a-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 summarised in Tab. 1. The force field by Ding et 10 al. correctly predicts the lowest energy for dimer structure "a", and the relative energy differences ∆∆E between optimised structures are closer to those obtained in the ab initio calculation than for the OPLS force field, which assigns the lowest energy to structure "d", which is the highest en-15 ergy 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 (Hintze et al., 2003;Chackalackal and Stafford, 1966;Miller 30 et al., 2005), the force field by Ding et al. 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 35 for transferable potentials, as opposed to the Ding et al. 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. for the collision simulations.
2.2 Potential of mean force of two sulfuric acid molecules 5 We first calculated the binding free energy of two sulfuric acid molecules in vacuum as described by the force fields of Loukonen et al. and Ding et al. The potential of mean force (PMF) as a function of the sulfur-sulfur distance was calculated from a well-tempered metadynamics simulation 10 ( 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 interac- 15 tions 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 B T were deposited every 500 steps along the collective variable and a harmonic wall was used to restrict it to values below 35 Å. A bias fac-20 tor of 5 was chosen, and a Nosé-Hoover chain thermostat of Table 1. Relative energies ∆∆E (eV) and hydrogen bond distances dO···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.   (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. length 5 with a time constant of 0.1 ps was used to keep a temperature T = 300 K. The combined length of the trajectories was 120 ns for each force field. Both PMFs, shown in

Collision simulation setup
Molecular dynamics simulations were performed with the 35 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 40 constitute a closed system (in atmospheric conditions collisions with the carrier gas are rare on the time scale 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 setup was used: 45 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 50 r Figure 3. Potential of mean force between two H2SO4 molecules as a function of the sulfur-sulfur distance calculated by metadynamics simulation for the OPLS force field (Loukonen et al., 2010)  the centre of mass motion of each molecule removed separately. Then the system was evolved for 50 ps, to randomise 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 55 the x direction, v x = ±v/2, were 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 60 relative velocities follow the Maxwell-Boltzmann distribution with reduced mass µ = m/2. We sampled relative velocities between 50 and 800 ms −1 , in steps of 50 ms −1 . 99 % of the distribution lies within this range at T = 300 K. For each value of the impact parameter b and the relative veloc-65 ity v, 1000 simulations were carried out starting with different initial atomistic velocities, to ensure good sampling. The simulation setup is very similar to the one recently used by Yang et al. (2018). Sulfuric acid molecules bind to each other through the formation of one or more hydrogen bonds. 70 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 75 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 formation of one or more hydrogen bonds. Three examples 80 of simulated trajectories with relative velocity closest to the mean velocity at 300 K (350 ms −1 ) and impact parameter of 8.5 Å are illustrated in Fig. 4.
The results from the atomistic simulations will be compared to different theoretical models described in the follow-5 ing.

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 10 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 emphasise this conceptual difference between simulations and theoretical models, we use the word capture instead of collision to refer to theory- 15 based results.
The interaction between two identical polar molecules is usually written as where d 1 and d 2 are the dipole moment vectors of the 20 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 −6 . The capture rate constant for such a potential can only be calculated numerically (Maergoiz et al., 1996c). In the 25 present section we consider only the isotropic part U (r) of the interaction described by Eq. (5), the effect of 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 30 point-like colliding particles in 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), where r is the distance of the colliding body from the centre of the field, L the angular momentum. Both the total energy of the system (which equals the initial translational energy µv 2 /2 at r → ∞) and the angular momentum are conserved. 40 The centrifugal term introduces an energy barrier, and for a successful capture at the barrier (r = r max ) the translational energy µṙ 2 max /2 has to be positive. Since the angular momentum equals L = µvb, the condition for b 2 to ensure a capture is In case of a simple attractive potential (repulsive forces can be neglected, as the studied velocities are relatively low), It is preferable to consider the squared value of b, since the capture cross section is calculated as σ c = πb 2 c . It is important to note that in the Langevin model, the total energy is divided strictly to the translational and potential energy, the 55 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. In the study by Kürten et al. (2014), a model of Brownian coagulation in a field of force (Sceats, 1986(Sceats, , 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 65 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 2 SO 4 /H 2 O particles with diameters of 49-127 nm, yielding a collisions enhancement factor value 5 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).

Brownian model of aerosol coagulation
For the attractive potential described by Eq. (8), the collision enhancement factor over the kinetic gas theory rate, W B = β B /β HS , from the Brownian coagulation model in the free molecule limit can be written as (Sceats, 1989) To compare the Brownian and the Langevin models to atom- 15 istic simulation results using the OPLS force field, we have determined the attractive potential parameters and r 0 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 20 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).

25
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 30 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 probabil- 35 ity curves are shown for each velocity separately. Predominantly, the collision probability approaches to unity with lower relative velocity and shorter impact parameter. However, the collision probability decreases slightly at high relative velocities where rapid re-dissociation of the complex 40 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 ms −1 , and small values of the impact parameter, the collision probability is also reduced. This happens because in some cases the fluctuations in the inter-45 molecular 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 more detailed discussion). The dynamical collision cross section, obtained from the 50 integral over the collision probability functions, 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 55 capture cross sections σ c , the velocity dependence of the change in σ d is in very close agreement with the Langevin model solution where v 0 is a reference velocity, as shown in Fig. 6. The 5 importance of contributions from long-ranged interactions to the collision cross section is evident, as σ c is proportional to v −4/n , for interactions decaying with r −n . Furthermore, as can be seen in Fig. B1, the critical impact parameters from the Langevin model are matching rather well with the tails of 10 the simulated collision probability curves, the intersection is located without exception at P (b, v) ≈ 0.2. 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 15 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 transfered to rotational and vibrational modes, and therefore the motion over the centrifugal barrier is hindered, and the critical impact parameter effectively reduced. Additionally, steric hindrance caused by intermolecular orienta-25 tions 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 on average 0.82 for the collision of two sulfuric acid molecules.

30
The canonical collision rate coefficient can be calculated similarly as Eq. (1), but since 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 35 both the relative velocity distribution f (v) and b 2 as 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 40 collision probability function is not necessarily temperatureinvariant. However, in Appendix B it has been 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. We therefore used the collision 45 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). In case of the Langevin model (Eq. (9)) the expression for the canonical capture rate coefficient can be simplified to Since the potential of mean force is required for the Langevin model, the coefficients are calculated only at 250, 300 and 400 K (see Appendix B for further details). In addition to the coefficients obtained by different ap-55 proaches, we examine the enhancement factor W relative to the kinetic gas theory rate expressed in Eq. (1), where a hardsphere radius R = 2.77 Å was calculated from the bulk liquid density of sulfuric acid, ρ = 1830 kgm −3 , assuming a volume fraction of one. Thus, after performing the integration, 60 the enhancement factor obtained using the Langevin model can be expressed analytically as where Γ(x) denotes the Gamma function, and and r 0 are the parameters based on the potential of mean force between two 65 sulfuric acid molecules in 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 70 300 K. Both the collision rate coefficients and the enhancement factors are shown in Fig. 7. As the hard-sphere collision rate is linearly increasing as the molecular velocity is proportional to √ T , the rate coefficient based on atomistic simulations is only slightly increasing with temperature due 75 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 is decreasing with increasing temperature. This is due to the neglect of the anisotropy of the dipole-dipole interactions 80 in the Langevin model: at higher temperatures, the effect of anisotropy becomes less important and therefore the model overestimates the collision rate less, compared to lower temperatures. The same effect can be observed in the Brownian coagulation model. Indeed, as seen from Fig. 7, the enhancement factors based on the theoretical models reach bet-5 ter agreement with the factor obtained from atomistic simulations as the temperature rises. As a result of the abovementioned reasons, the enhancement factor over the hardsphere collision rate coefficient is lower at higher temperatures regardless of the chosen approach.
The enhancement factor obtained by atomistic simulations is in very good agreement with the kinetic modelling on recent experimental results of formation of atmospheric sulfuric acid dimers (Kürten et al., 2014) and small clusters of sulfuric acid, dimethylamine and water (Lehtipalo et al., 2016;15 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 experiment, as mentioned earlier, whereas according to our molecular model of long-range interaction the basic Brownian model using the 20 attractive part of the potential of mean force overestimates the rate enhancement factor by about 40 %.

Conclusions
In summary, we have benchmarked two classical force fields against experimental and ab initio data and determined that 25 the OPLS force field by Loukonen et al. 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 vacuum by molecu-30 lar dynamics simulations and compared our results against simple theoretical models. We have found that the effective collision cross section of two H 2 SO 4 molecules, as described by the OPLS force field, is significantly larger than the hardsphere diameter assigned to the molecule based on the liq-35 uid density of sulfuric acid. As a consequence, we find the collision coefficient for sulfuric acid molecules is enhanced by a factor 2.2, compared to kinetic gas theory at 300 K. This enhancement factor obtained from atomistic simulation is consistent with the discrepancy observed between exper-40 imental 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 depen-45 dence 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 a bit overestimated due to the imperfect treatment of 50 the dipole-dipole interaction, yet 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 55 presented in this work can be applied to other atmospherically relevant molecules, clusters, or ions, exhibiting dipoles of varying magnitude -and in some cases several times larger than the one of the sulfuric acid molecule -to help understand the effect of long-range interactions in cluster for-60 mation 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 ACDC (McGrath et al., 2012)), it is necessary to obtain the enhancement factors for all the relevant 65 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. Future work therefore should also be aimed at finding simple models for predicting approximate rate enhance-70 ments, 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 75 Maergoiz et al. (1996c), using classical trajectory integration, calculated the capture rate constant when two identical polar molecules interact through a potential containing anisotropic (∝ r −3 ) and isotropic (∝ r −6 ) terms, where d 1 and d 2 are the dipole moments 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 85 field is given by where the thermal capture cross section can be calculated using a fitting function κ(θ, M ) as 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Å 2 , which deviates slightly from 10 values 100.66 and 104.94 amuÅ 2 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 it is expressed as where z = a 2 + ln θ.
The reported fitting parameters are (Maergoiz et al., 1996c) Since 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-25 dipole interaction from the isotropic interaction, the constant C is written as where f is a factor denoting the relative magnitude of the anisotropic interaction between permanent dipoles with re-30 spect 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 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 40 2 % (Maergoiz et al., 1996c)).  Since we are unable to distinguish the actual dipole-dipole interaction from the total attractive potential, we have estimated the interaction using the Keesom equation (see Fig. 3): According to Eqs. (8) and (A11), about one third of the attractive potential is due to dipole-dipole interaction, and consequently the enhancement factor is lower than for the isotropic field. The estimated rate enhancement factors at 50 250, 300 and 400 K are shown in Tab. A1. Thus, by taking into account the anisotropy of the intermolecular potential, the estimated capture rate is in better agreement with the result obtained using atomistic simulations. Relative velocity (ms ) 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 simulation 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 2 = (2R) 2 ), or obtained from the Langevin capture model (Eq. (9)), are indicated by the solid black, and dashed coloured lines, respectively.
Appendix B: Temperature dependence of collision 55 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 have tested the effect of the different rotational and vibrational motion on the collision statistics in an atmospherically relevant temperature 5 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 K and 400 K, instead of 300 K. As shown in Fig. B1, such moderate change in temperature in-10 deed does not affect the collision probabilities between two sulfuric acid molecules. 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 15 was calculated at 250 K, 300 K, and 400 K, and the parameters describing the attractive intermolecular interaction (Eq. (8)) are reported in Table A1. 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 in-25 termolecular 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 30 energy can alter the translational motion and eventually lead to a definitive separation of the molecules. 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 ms −1 where a collision occurs 35 and a second one where the molecules are repelled at range. While in both cases the fluctuation of the intermolecular energy exceeds the initial translational energy of 1.27 meV, in the trajectory exhibiting a repulsion, the large positive energy fluctuations are longer lived and dominate the interaction.

40
Author contributions. RH and BR planned the simulation setup and performed benchmark calculations. BR carried out collision and metadynamics simulations and wrote the first draft of the paper. RH and EZ provided the theoretical framework. RH analyzed the simulation data. TK and HV helped plan the project and all authors 45 contributed to writing the manuscript.
Competing interests. The authors declare that they have no conflict of interest.