Simple Physical Model for the Estimation of Irreversible Dissociation Rates for Bimolecular Complexes

In this article, we propose a simple method of estimating dissociation rates of bimolecular van der Waals complexes (“wells”), rooted in rigid body dynamics, requiring as input parameters only the bimolecular binding energy, together with the intermolecular equilibrium distance and moments of inertia of the complex. The classical equations of motion are solved for the intermolecular and rotational degrees of freedom in a coordinate system considering only the relative motion of the two molecules, thus bypassing the question of whether the energy of the complex is statistically distributed. Well-escaping trajectories are modeled from these equations, and the escape rate as a function of relative velocity and angular momentum is fitted to an empirical function, which is then integrated over a probability distribution of said quantities. By necessity, this approach makes crude assumptions on the shape of the potential well and neglects the impact of energy quantization, and, more crucially, the coupling between the degrees of freedom included in the equations of motion with those that are not. We quantify the error caused by the first assumption by comparing our model potential with a quantum chemical potential energy surface (PES) and show that while the model does make several compromises and may not be accurate for all classes of bimolecular complexes, it is able to produce physically consistent dissociation rate coefficients within typical atmospheric chemistry confidence intervals for triplet state alkoxyl radical complexes, for which the detailed balance approach has been shown to fail.


Centrifugal Correction to Probability Distribution
The centrifugal distortion of rotational energy levels is a well-known effect in spectroscopy, in which we may express the (first-order) distortion as a perturbative correction to the rigid rotor energy (Equation 14 in the main article): In this section we investigate if centrifugal distortion has any impact of the probability distribution of angular momentum, by deriving the centrifugal distortion constant F for a Lennard-Jones potential well as well as a a suitable cut-off value l c for the probability distribution. Bimolecular complexes are obviously non-rigid rotors, otherwise they would not be able to dissociate in the first place. Since the Lennard-Jones expression we use for the interaction energy of bimolecular complexes is admittable crude, we will first derive the first-order centrifugal distortion term using a general potential energy before plugging in our Lennard-Jones approximation into it.
First, we have to solve the minimum of the effective potential V ef f (r) = V (r) + L 2 φ 2µr 2 as a function of L φ . We will call this minima r L . Assuming r e ≈ r L at all physically reasonable values of L φ , we take the Taylor series of V ef f (r) in the vicinity of r e up to the second order: S-2 Comparison of this approximation with numerically determined values of r L for the floppier (MetO-XO) complexes shows that the approximation slightly underestimates the deviation from r e starting from L φ values corresponding the final 5-10 percentiles of ρ(L φ , T ) at room temperature. For the stronger binding complexes (α-pin-derived) the approximation is accurate for all modelled values of L φ . As higher values of L φ contribute more to dissociation rates, the third order Taylor series term might be important for some systems. However, the resulting expression for r L is a complicated quadratic equation solution, which we will not consider here. Now, let's see what corrections result from applying our Lennard-Jones potential: V LJ (r) = D r e r 12 − 2 r e r 6 = D (r 12 e − 2r 6 e r 6 ) r 12 (3) r 2 e into equation 2 results in: Next, we want to expand V ef f (r L ) = V (r L ) + : We are able to use the fact that dV dr re = 0 to our advantage. Differentiating equation 4 with L 2 φ results in: Resulting in a first-order centrifugal correction to the rotational energy of: From this, we can derive the 1st order centrifugal coefficient Θ D (for a Lennard-Jones well) in Equation 18 in the main article: As mentioned in the main text, in one includes a first order centrifugal distortion into the rotational energy, the probability distribution ρ(L φ , T ) becomes divergent, and we thus need an integration limit below infinity. We call this quantity the critical angular momentum L c , named after the critical temperature because we are deriving it using the same method as critical temperatures and pressures are derived from the Van der Waals gas law. S1 We look for an inflection point in V ef f (r) at which both the first and second derivatives are zero: This equality of these two functions should only apply at a specific value of L φ . Thus we find the critical angular momentum at the intersection of the following two functions of r: For a Lennard-Jones well: Values for all centrifugal distortion-related parameters are found in Table S6. A visualisation of the difference between ρ D (l, T ) (Equation 18 in the main article) and the centrifugally uncorrected ρ(l, T ) is shown in Figure S1, based on which it was judged that centrifugal distortion corrections will not impact the average dissociation rates significantly. Figure S1: The difference between the corrected and uncorrected probability distributions of angular momentum. All observed changes are less than 10 −3 of the peak value of the respective uncorrected distribution (0.00214 s m for (MetO) 2 .)

Variation of Cutoff Distance
Trajectories were simulated at the l = 0 level of theory for a smaller sample of complexes to test the sensitivity of the results to the chosen cutoff distance. The results are seen in Table S1. The dependence of the dissociation rate on the cutoff distance r c should eventually converges to k d (T ) ∝ r −1 c one r c is far enough, but this seems to not be quite true yet at cutoff distance 20-30Å. A noticeable trend in these results is a slower convergence to the stated dependence for the complexes with further equilibrium distances r e , which is consistent with S-5 the fact that the potential gradient is higher closer to equilibrium ( Figure S2). The trends seen here are the reason why the cutoff distance is set to depend on r e in the trajectory simulation code, as this means that the ratio V (r c )/D will always be the same (≈ −0.0027) at the cutoff for each trajectory. (In)accuracy of Lennard-Jones Potential As covered in the main article, the Lennard-Jones potential (Equation 8) was used to model the total interaction potential for purposes of simplicity rather than accuracy. The main assumption behind this choice was that the depth of the potential well has a much larger impact on the overall dissociation rate than the curvature of the interaction potential at distances r > r e . In this section we will interrogate that assumption. Many popular force fields in Molecular Dynamics make use of the electric field multipole expansion for intermolecular interactions, resulting in a function with the general shape 9. We therefore find it noteworthy that this model is accurate for long-range interactions only, at which the intermolecular distance significantly exceed the intramolecular distances. S21 In other words, this is not a good model for close-range intermolecular interactions. Unfortunately, as seen in the typical interaction potential curve sketched in Figure S2, the close-range interactions 1 The cited source discusses this in terms of systems composed of multiple classical charges, but in our case the extent of orbital overlap is a better measure for when the long-range regime starts being a good model for intermolecular interactions. This discussed briefly in source S3 S-6 of the molecules are also the most important for determining the relative velocity of the molecules along the dissociative trajectory, as that is where the ∇V (r) is the largest. This means that explicit quantum chemical calculations varying the intermolecular distance r are required to determine V (r) with reasonable accuracy. Figure S2: A visualisation of why long-range approximations of the potential energy are not accurate for calculating dissociation.
The accuracy of the Lennard-Jones function by performing quantum chemical potential energy scans on three reasonably small bimolecular complexes with different intermolecular bonding properties: (MetO) 2 , (AceO) 2 and MetO-ProOHO. The scans were performed in Gaussian 16 S4 as a series of constrained geometry optimizations with the two radicals at different intermolecular distances generated from the global minimum geometry provided in the original source. The distance between the two centers of mass was calculated for all the optimized geometries. The electronic structure method chosen for the geometry optimization was ωB97X-D/jul-cc-pVDZ, S5,S6 which is one 'month' and one basis function step down S-7 from the basis set used in reference S11. The resulting (electronic) energies were compared to those calculated for the optimized bimolecular complex and the isolated radicals using the same level of theory, and the results are presented in Figure S3. As seen in the figure, The attractive interaction is clearly stronger than r −6 for (MetO) 2 and (AceO) 2 , whereas for MetO-ProOHO it is somewhat weaker. This is somewhat counterintuitive, as the latter is a H-bonded complex, but a closer look at the data points above the LJ curve showed that the ProOHO molecule prefers prefers an intramolecular H-bonding structure at these distances, explaining the sharper decline in binding energy as a function of r.
For each of the three complexes, a simple trial function that approximates the V (r) data better than the Lennard-Jones potential was chosen. Dissociation trajectories were simulated on levels of theory 1 & 3 and dissociation rates were calculated. The used trial functions are presented in Equation 10, and compared to scan data in Figure S4. The dissociation rates calculated based on the trial functions are presented in Table S2. As seen in the table, the differences are well within the margin of error, justifying the use of the LJ potential for all complexes.   Table S2: Comparison of dissociation rates calculated using the Lennard-Jones potential and the trial potentials.
Firstly, we assume the orientation-dependent interaction potential V (⃗ r, O) is separable into an isotropic component V (r) and a dimensionless 'anisotropy factor' f (O) 2 .
First we must address the most common method of modelling long-range potential surface anisotropy in physical chemistry, namely using the dipole and quadrupole moments of the molecules to quantify the orientation dependence of V (r, O), because it is unfortunately inadequate for our purpose. There are two reasons for this: Firstly, the usage of multipole moments determined for independent molecules at equilibrium geometry implicitly assumes that geometry stays rigid. In van der Waals complexes, the intermolecular interactions perturb both molecular structures from their respective independent equilibrium geometries, meaning that the multipole moments determined for that geometry is no longer accurate.
Secondly, the concept of a dipole moment is only physically well defined at distances noticeably longer than the dipole's characteristic length. S2 This is most certainly not the case in van der Waals complexes. In other words, one must apply a model that fully accounts for the fact that the molecules are polarizable diffuse charge distributions.
A better close-range anisotropy model would be to consider consider the three intermolec-2 This assumption may not be applicable if one of the molecules has significant orientation dependence in polarizability, resulting in V (r) decaying differently depending on the orientation.
S-10 ular DOF which couple into free rotational modes at high separation, that is, the torsional vibrations. Assuming these are independent, we can most simply model these as three rigid hindered rotors. With that model, the potential energy at equilibrium distance is: S7 Where U i are the hindered rotor barrier heights. As these are intermolecular torsional modes, σ i is determined by the rotational symmetry of the molecules. Note also that one of the three modes must correspond to the 'altitude coordinate', spanning onlt 0 ≤ θ 3 ≤ π rather than 0 ≤ θ 1,2 ≤ 2π. S8 According to our simplified Lennard-Jones model (Equation We may thus reconcile these two models by plugging Equation 11 into D: Were we have placed the U 1 + U 2 + U 3 expression outside the sum to emphasize that it is a constant. The radial equation of motion changes to: And, if we neglect the quantization of energy levels, we can formulate classical equations of motion for the torsional DOFs as well: S-11 By assuming that the initial conditions are Boltzmann-distributed, we may derive a canonical probability distribution of (initial values of) the torsional positions θ i : For the first two degrees of freedom, the normalization constant is a modified Bessel function of 0th order: 2kT cos(σ 1,2 θ 1,2 ) σ 1,2 dθ 1,2 For the third torsional coordinate the normalization constant is a hyperbolic cosecant function: From this point onward, one samples different initial values for the three torsional coordinates similarly as we have done for v and l. As seen from the equations of motion, the torsional coordinates are interdependent with the distance coordinate r, which suggests that the anisotropic effects of the dissociation are nonlinear and hard to quantify. The timestep of the simulation might have to be shortened to ensure that energy is conserved during the simulation. Another added difficulty arising from this approach is that the torsional degrees of freedom have to be characterized and analysed in the first place, and most Hindered Rotor characterization codes are better equipped for finding rotations of covalent bonds.

S-12
∆G for a complex with zero binding energy For a reaction with net production of gas molecules, the enthalpy change is ∆H = ∆E + RT ∆n in constant temperature. As such the dimerization Gibbs free energy for a hypothetical complex with zero binding energy is: In general E is a sum of the zero-point-corrected potential energy E 0 and the thermal energy E T . However, with zero binding energy ∆E 0 is zero, and ∆E comes exclusively from the thermal component. Here the partition function-based definition of thermal energy and entropy is useful: With no binding energy, the six intermolecular modes are all unbounded, and are most accurately described as either translations or rotations. Thus the components of the partition functions are described by: where RR refers to the rigid rotor approximation. For the translational modes we use note that we can use the formula for reduced mass to our advantage, as Using the reduced mass of a (MetO) 2 complex, Q −1 t,µ = 4.165 · 10 −7 when p = 1 atm and T = 298.15 K. For the rotational partition function we make two admittedly contradictory assumptions for simplicity: We assume that the rotational symmetry number is σ = 1 and that both molecules are uniformly dense (ρ(r, θ, φ) = ρ 0 ) spheres of radius re 2 and total mass m = ρ 0 V , and as such all three moments of inertia are: The 'bispherical complex' is this case is a fully symmetric prolate rotor, for which the lower moment of inertia is I 12,A = 3mr 2 e 20 , and the two higher are I 12,B = I 12,C = µr 2 e = mr 2 e 2 . The quotient of rotational partition functions is: Which, using the mass of a MetO radical and the equilibrium distance of a (MetO) 2 complex is Q r12 Q r1 Q r2 = 1.278 · 10 −4 . The total Gibbs free energy of dimerization is: Where we assumed ω 12 ω 1 ω 2 = 3 4 as with the triplet state alkoxyl complexes. This 'zeroenergy limit' is consistent with typical binding entropies for atmospheric (strongly bound) binding entropies, which are between -29.3 and -39.8 cal mol The detailed balance dissociation rate is then: cm 3 molecule s · 2.46 · 10 19 molecule cm 3 e 22 = 8 · 10 18 s −1 S-14

Complex parameter data
The main physical parameters for the model complexes are presented in Table S3. Physical parameters related to angular momentum are presented in Table S4, and model-specific derived parameters are presented in Table S5.   is the primary rotational constant expressed in units of l, presented to give an idea of the accuracy of the assumption that the probability distribution is continuous. As touched on in the main article, the variance of the minimal escape velocity on l was determined numerically and fit to the following empirical second-order polynomial, where the parameters α and β are presented in Table S5: S-16

v,l Probability distributions visualised
This section includes visualized probability distributions for all bimolecular complexes. All non-dissociative trajectories, that is, v, l-values for which v < v c (l), are cut out from the figures, just like in the code.

xyz geometries used for parameter determination
The alkoxy radical complex geometries are from sources. S11-S14 Citations are presented in the table captions. r e -parameters were determined using the formula for centre of mass: In the equation i represents a molecule and j represents an atom. The geometries are presented below. Horizontal lines are used to show where one molecule ends and the other begins. S-34

Non-Ergodicity of Bimolecular Complexes
The possible non-ergodicity of bimolecular complexes impacts the accuracy of the typical models we might use to model the dissociation. In this section we will discuss the applicability of the detailed balance and ergodicity assumptions in terms of the triplet state (MetO) 2 dimer, the simplest bimolecular complex we are interested in. The detailed balance assumption of association and dissociation holds if the equilibrium energy distribution (Boltzmann distribution) of the dissociating system and the steady-state energy distribution do not differ significantly: S15 The cited paper formulates the same argument more rigorously in terms of energy relaxation eigenvectors of the master equation matrix, but the substance of the argument stays the same: Detailed balance breaks when the steady-state energy distribution differs significantly from the Boltzmann distribution, and the important physical parameters for determining this is the ratio between the chemically significant eigenvalue (the dissociation rate) and the lowest energy relaxation eigenvalues (the slowest rates). If the former sits firmly outside the numerical continuum of the latter group, then dissociation rates calculated using Equation 2 are reliable.
If the dissociation outspeeds the intermolecular energy relaxation, than the complex is non-ergodic, in which case kinetic methodologies derived using statistical ensembles (The Eyring equation, RRKM, etc.) are all inaccurate, at least if the the complex as a whole is treated as a system in equilibrium. It has long been well-known that dissociation reactions with weakly bound complexes have this property, complexes with noble gas bonds being a particularly strong example. S16 Marcus himself made a distinction between weakly and S-56 strongly bound complexes when discussing the applicability of the RRKM model for dissociation reactions. S17 All of this depends on the magnitude of the anharmonic couplings. Energy flows from one oscillator to another due to the existence of non-linear resonances, which depends on the magnitude of the anharmonic couplings. We investigate this by performing an automated VPT2 S18 anharmonicity analysis of the (MetO) 2 vibrational modes using ωB97XD/aug-cc-pVTZ, S5,S19 the same level of quantum theory that the geometry was optimized with in Source. S11 This calculation was performed in the Gaussian 16 program suite. S4 The total anharmonic energy, excluding the rovibrational coupling is thus described by the series expression (in mass-weighted atomic units): where K are derivatives of the potential energy.   Tables S32-S36. The precise dynamics of these coupled systems is highly complex, S20 and exact quantum mechanical solutions exist only for model systems. S21,S22 One useful model system for our case, however, is a recent study by Zhang et.al, S23 focusing on modelling a 'ergodicity phase diagram' for a system composed of two fragments interacting via a Bose Statistics Triangle Rule potential, S24 in which the cubic anharmonic couplings are expressed as 1 Here V rx 'is on the order of a typical anharmonicity constantχ ij ', which in our case seems to be on average 77 cm −1 for the intermolecular couplings (Table S36). The a parameter here describes the ratio of higherorder coupling terms compared to the lower-order terms. It it typically between 0.1 and 0.3 depending on the 'softness' of the vibrational modes, coupled methyl rotations being on the softer end of the scale. In Fig. 4 of Source S23, a 'ergodicity phase diagram' is presented S-59 for various coupling strengths. If the coupling potential V rx is taken as the strength of the intermolecular modes in our case, we are located in the region Vrx Ω = 77 cm −1 1390 cm −1 ≤ 0.06 of the diagram, in which even partial ergodicity requires a very 'soft' coupling of a ≥ 0.25.
Somewhat problematically, the modelled diagrams all have higher excitation energies than our complexes would realistically have, as the 'low' excitation energy of 2000 cm −1 exceeds the dissociation energy of (MetO) 2 by over 70 %. Thus we may conclude that the white 'non-ergodic' area of the diagram will be further shifted to the right, requiring exceptionally strong coupling for ergodicity to apply. As seen Tables S32-S36, the coupling between the torsional modes of the −CH 3 groups and the stretching of the CH bonds is quite strong, and these may or may not provide the necessary energy flow for partial ergodicity. If this is the case, it should be noted that this is due to the proximity of the two methyl groups, a property the (MetO) 2 complex doesn't share with many other alkoxyl radical complexes.
Thus we may conclude that dissociation of the complex at the very least poses a serious competition for internal energy relaxation, and that equilibrium models for calculating the dissociation rate like Detailed Balance or RRKM may overrestimate the amount of available energy the dissociative modes have at their disposal.