Clustering effects on the reactivity of alkoxy radicals: rate coefficients of ³(RO⋯OR) complexes accounting for multiple conformers

Abstract

Atmospheric peroxy (RO₂˙) and alkoxy (RO˙) radical species are crucial intermediates in the formation of secondary organic aerosol (SOA). Recent computational work has demonstrated that recombination reactions of peroxy radicals (RO₂˙ + RO₂˙) proceed through triplet complexes consisting of two alkoxy radicals (³(RO⋯OR)). To understand how peroxy recombination reactions branch into different product channels, it is thus necessary to thoroughly investigate the reactions of these triplet alkoxy complexes. Although the reactions of free alkoxy radicals have been extensively studied, the reactivity of triplet alkoxy complexes remains relatively less explored. In this study, we have therefore developed a systematic conformer sampling workflow for ³(RO⋯OR) and applied it to four typical alkoxy systems (AceO, β-ISOPO, PhCH₂O, and PhC(O)O). Rate coefficients (k) of key reactions have been calculated using multi-conformer transition state theory (MC-TST) and lowest-conformer transition state theory (LC-TST), allowing the quantitative evaluation of conformer effects. Our results demonstrate that the presence of the other RO˙ in the complex has a noticeable effect on k values. For the β-ISOPO system in particular, the predicted β-scission k values in the ³(RO⋯OR) complex are over 100 times higher than those of the free radicals.

---

1 Introduction

Climate change and air pollution have become major challenges in current times. As a pollutant, small atmospheric particles account for around 10 million global deaths and 10⁸ disability-adjusted life years annually.^1,2 On the other hand, atmospheric aerosols can cool the climate through cloud formation and scattering.³ Sustainable development actions have significantly reduced the emission of inorganic species (e.g. NO_x, SO₂, etc.), initiating a shift in aerosol compositions: before those actions, atmospheric aerosols were dominated by primary inorganic aerosols, while after those actions, secondary organic aerosol (SOA) particles have become more prominent.⁴

Volatile organic compounds (VOCs) are a major source of SOA, and come from various biogenic (e.g., isoprene molecules released by vegetation)⁵ and anthropogenic (e.g., fuel combustion, usage of organic solvents)⁶ emissions. After emission, VOCs are oxidized and converted into peroxy radicals (RO₂˙). Depending on its molecular structure, RO₂˙ can undergo both unimolecular reactions, such as H-shifts and ring closures, and bimolecular reactions with species, including NO_x, HO_x and RO₂˙.⁷ Among these channels, we are especially interested in the peroxy recombination reaction ((R1) in Fig. 1), as some subsequent channels can generate accretion products⁸ (e.g. (R2.3) and (R2.4) in Fig. 1), which are critical for SOA formation. Compared to the initial VOCs, accretion products have a higher number of carbon atoms and lower volatility,⁹ allowing them to cluster very efficiently.

Fig. 1 The known product channels of peroxy recombination reactions.

As shown in Fig. 1, the peroxy recombination reaction proceeds via a tetroxide intermediate, and generally† leads to a triplet alkoxy complex (³(RO⋯OR)) after releasing one ³O₂ molecule.^11,12 The ³(RO⋯OR) complex consists of two alkoxy radicals (RO˙) and can undergo various reactions. Therefore, the branching of peroxy recombination channels reflects different fates of ³(RO⋯OR).

So far, four reaction channels of ³(RO⋯OR) have been discovered: dissociation¹³ ((R2.1) in Fig. 1), intermolecular H-shift^14,15 ((R2.2) in Fig. 1), intersystem crossing (ISC)¹⁶ to the singlet surface ((R2.3) in Fig. 1), and β-scission of one of the alkoxy radicals^17,18 ((R2.4) in Fig. 1). Rate coefficients (k) reported for these channels vary widely with both the structure of ³(RO⋯OR) and the calculation method. Nevertheless, all four channels share a common feature: they can proceed very rapidly, with estimated k(298 K) values often exceeding 10⁹ s⁻¹.^14,15,17

In addition to β-scissions, other intramolecular reactions may also take place depending on the structure of the reactive RO˙. These intramolecular reactions include CO₂ elimination in acyloxy radicals,¹⁹ intramolecular H-shifts^20,21 and epoxide formation.^22,23 We expect that any sufficiently fast intramolecular RO˙ reaction²⁴ can also occur in ³(RO⋯OR)—but the in-complex rate is likely to differ from that of the free-radical reaction.

Among the reaction channels of ³(RO⋯OR), ISC and β-scission are of special importance, as they may lead to recombination processes that form accretion products. In general, the ester or ether products from β-scission (ROR′) are expected to be more stable and, therefore, longer-lived than the peroxide products from ISC (ROOR).¹⁷ However, note that neither ISC nor β-scission reactions guarantee the accretion product formation, as other competing processes may intervene before recombination. For example, H-shift and dissociation reactions can take place in ¹(RO⋯OR) (from ISC)²⁵ or ³(R′⋯OR⋯CH₂O) (from β-scission) complexes. Additionally, the ROR′ formation after β-scission probably requires a subsequent ISC.¹⁷

Given the importance of atmospheric ³(RO⋯OR) complexes, a natural progression is to investigate the rate coefficients of their different reactions. Numerous experimental studies have probed the final product distribution of peroxy radical recombination. For example, Frandsen et al.¹⁸ recently studied the ozonolysis of tetramethyletylene (TME), where the dominant first-generation RO₂˙ is the acetonyl peroxy radical, CH₃C(O)CH₂OO˙ (AceOO˙). Products with the elemental composition of C₆H₁₀O₄ were abundantly observed, suggesting ISC as the major channel. On the other hand, Peräkylä et al.¹⁷ detected C₁₉ esters as the primary accretion product of α-pinene (C₁₀) ozonolysis, implying that β-scission dominates the accretion product formation.

Although experiments can probe the final product distribution of the RO₂ + RO₂ reactions and even have recently captured the tetroxide intermediate in the gas phase (observed as R₂O₄(H₂O)H⁺),²⁶ current techniques are still unable to detect ³(RO⋯OR) directly, as this complex is too short-lived for spectroscopic detection or isolation. Alternatively, theoretical methods can reveal the structures of key intermediates and estimate the stepwise reaction rates.

Previous theoretical studies^14,15 have shown the difficulties of studying the in-complex reactions: (1) the ³(RO⋯OR) complex, an open-shell species with two unpaired electrons, is different from the benchmarking species used by major computational chemistry methods, up to and including state-of-the-art methods such as coupled cluster. It is noteworthy that the training sets of the xTB methods, which are commonly used for semi-empirical conformer sampling, are primarily closed-shell species.^27–29 (2) The structural variety of ³(RO⋯OR) conformers affects the accuracy of the predicted thermodynamic and kinetic parameters. By definition, each conformer corresponds to a distinct potential energy minimum.³⁰ Energetically low-lying conformers influence the rate coefficient as their Boltzmann population is non-zero.^31,32 Previous computational studies have suggested that ignoring conformers other than the global minimum may lower the accuracy of the rate coefficient prediction, especially at higher temperatures relevant to combustion.^31,33 Furthermore, compared to RO˙, a larger number of conformers can be expected for ³(RO⋯OR) due to intermolecular interactions.

Based on the discussions above, to compute the reaction rate coefficients of ³(RO⋯OR) more accurately, a systematic conformer search workflow is required. In this study, we improve the conformer sampling workflow reported by Hasan et al.^14,15 and Møller et al.,³¹ with a special focus on intermolecular H-shift reactions. Since these reactions involve both RO˙ in the ³(RO⋯OR) complex and certain computational tools (e.g. GOAT) had not been developed at that time, previous studies^14,15 encountered challenges to sample the transition state (TS) conformers of intermolecular H-shift reactions comprehensively. Then, we apply our updated workflow to four ³(RO⋯OR) systems corresponding to representative RO₂˙ + RO₂˙ pairs in the atmosphere. We determine the k values of β-scission, intermolecular H-shift and other selected reactions, such as RO˙ addition to a double bond (R3.2.2(a) and (b) in Fig. 3) and epoxy formation–decomposition (R3.3.2(1) and (2) in Fig. 4). In addition to investigating the relevance of these novel pathways, our main purpose is to assess how clustering with another RO˙ affects the reactivity of alkoxy radicals, which is achieved by comparing the free-radical and in-complex k values.

This study does not include the rate coefficients of dissociation and ISC channels. For the dissociation channel, k strongly depends on the binding energy (D) of ³(RO⋯OR),¹³ and the D values are reported in Table S2 of the supplementary information (SI) (relevant notations are defined in Table S1). For the ISC channel, determining both the ISC rate and the post-ISC product distribution requires additional calculations on the highly multireference open-shell singlet surface, which is beyond the scope of this work.

The names and RO˙ molecular formulas of the four systems are listed in Table 1. Each system includes both free radicals and triplet complexes. The simplified notations for ³(RO⋯OR) are similar to those of their corresponding RO˙.

Table 1 Systems involved in this study and the names used to refer to them throughout this paper

Simplified notation in this paper	Full name	Molecular formula of RO˙	Structure of RO˙
AceO radical or AceO˙	Acetonyl alkoxy radical	C3H5O2
β-ISOPO radical or β-ISOPO˙	(S)-1-Methyl-1-vinyl-2-hydroxyethyl alkoxy radical	C5H9O2
PhCH2O radical or PhCH2O˙	Benzyl alkoxy radical	C7H7O
PhC(O)O radical or PhC(O)O˙	Benzyl acyloxy radical	C7H5O2

---

These four systems were selected based on experimental evidence and atmospheric importance: AceO˙ can be formed in tetramethylethylene (TME) ozonolysis, and several experimental studies^18,34,35 have investigated its product distribution. β-ISOPO˙, as the name implies, can be generated in the ˙OH-oxidation processes of isoprene.³⁶ Since stereochemistry is not the focus of this study, we selected only (S)-enantiomers of β-ISOPO˙ for simplicity. PhC(O)O˙ and PhCH₂O˙ were selected as representative aromatic alkoxy radicals as they are formed during the oxidation of toluene. In brief, Cl˙, ˙OH or NO₃˙ initiate a methyl H-abstraction reaction from toluene,³⁷ followed by O₂ addition to form PhCH₂OO˙. Recombination reactions of PhCH₂OO˙ then generate the corresponding alkoxy radical, PhCH₂O˙, as well as benzaldehyde (PhCHO).³⁸ Subsequently, PhCHO can be oxidized into PhC(O)O˙.³⁹

2 Computational details

2.1 Notations and computational chemistry programs applied in this study

For clarity and succinctness, Table 2 summarizes the simplified notations of the types of reactants and transition states, computational methods and energy terms used in this study. Detailed discussion of computational methods is provided in Sections 2.2–2.4. A schematic representation of the conformer sampling workflow is reported in Fig. S1.

Table 2 Simplified notations of conformers, computational methods and energy terms used in this study

Simplified notation in this paper	Detailed descriptions
Reactant monomer	An isolated RO˙ radical.
TS monomer	A transition state conformer of intramolecular RO˙ reactions.
Reactant dimer	A triplet alkoxy complex, ^3(RO⋯OR).
TS H-shift conformer	A transition state conformer of the intermolecular H-shift reaction of ^3(RO⋯OR), which involves both RO˙ in the complex.
TS dimer	A transition state conformer of other intramolecular reaction channels of ^3(RO⋯OR), which takes place within one RO˙ of the complex.
B3LYP	B3LYP/ma-def2-SVP.
ωB97X-D3 (Grid2)	ωB97X-D3/ma-def2-TZVP, with the default grid size of ORCA 6.0.1.
ωB97X-D3 (Grid3)	ωB97X-D3/ma-def2-TZVP, with the largest predefined grid size of ORCA 6.0.1.
DLPNO-CCSD(T)-F12	UHF-DLPNO-CCSD(T)-F12/cc-pVDZ-F12, with the keyword “tightPNO” and “tightSCF”, cc-pVDZ-F12-CABS as the auxiliary basis set, and aug-cc-pVDZ/C as the RI approximation basis set.
f(DFT)	Values calculated at the ωB97X-D3/ma-def2-TZVP (DefGrid3) level of theory, where f can be an energy term (E) or rate coefficient (k).
f(CC)	Values calculated at the UHF-DLPNO-CCSD(T)-F12/cc-pVDZ-F12 level of theory, where f can be E or k.
Eel	Electronic energy.
Eelzc	Zero-point corrected electronic energy.
G	Gibbs free energy.
D	Dissociation energy of a complex.
Esp	Single-point energy.
ΔbE	The energy barrier of a reaction (ΔbE = ETS − Ereactant), which reflects the kinetic feasibility of the reaction.
ΔrE	The reaction energy (ΔrE = EProduct − EReactant), which reflects the thermodynamic feasibility of the reaction.

---

The initial structures of isolated RO˙ were built using SPARTAN′24.⁴⁰ Reactant and transition state (TS) dimers were generated using the artificial bee colony algorithm for cluster global optimization (ABCluster) 3.3 program.^41,42 Transition state structures of intermolecular H-shift reactions (TS H-shift) were generated using the global optimization algorithm (GOAT)⁴³ implemented in ORCA 6.0.1.⁴⁴ Semi-empirical calculations were performed using XTB 6.7.1.^45,46 Calculations at higher level of theories, including density functional theory (DFT) and coupled cluster methods, were performed using ORCA 6.0.1.⁴⁴ Filtering and analysis of the results were done using the Jammy Key for Configurational Sampling (JKCS) 2.1 program.^47–49

2.2 Method for rate coefficient calculation and relevant energy barriers

Earlier studies reported that β-scission and intermolecular H-shift reactions of ³(RO⋯OR) proceed through a single transition state,^14,15 and can thus be treated as elementary reactions.⁵⁰ In addition, reactions taking place within a ³(RO⋯OR) complex do not involve collisions or exchanges with other molecular entities, and can therefore be treated as unimolecular elementary reactions. These reactions are further categorized as follows: if the reaction involves only one RO˙ in the ³(RO⋯OR) complex, it is a unimolecular intramolecular reaction. In contrast, if the reaction involves both RO˙ in the complex, it is a unimolecular intermolecular reaction. The first category includes β-scission, epoxy formation and decomposition, while the second category includes intermolecular H-shift and RO˙ addition to the double bond.

Therefore, reaction rate coefficients in this study are calculated using the elementary transition state theory (TST).^51,52 Although absolute k values predicted by elementary TST are often inaccurate, the relative rates are more reliable (e.g., the factor, which is further discussed in Section 3.5). In addition, according to Hasan et al.,¹⁴ variational effects are unlikely to be a major source of error even for rapid reactions in alkoxy complexes. The differences between their H-shift rate coefficients computed using elementary and variational TST are generally within a factor of two,¹⁴ suggesting that elementary TST is capable of describing such systems, and consequently, sufficient for this study.

In practice, elementary TST is applied in two forms when calculating k:³¹ multi-conformer TST (MC-TST),^31,53 which accounts for contributions from multiple conformers (within an energy cutoff) in eqn (1); or lowest-conformer TST (LC-TST),³¹ which includes only the conformer with the lowest (free) energy value, and thus with the highest Boltzmann probability, in eqn (2).

In eqn (1) and (2), κ_min is the tunneling coefficient, T is the temperature, k_B is the Boltzmann constant, h is the Planck constant, R is the gas constant, G_TS is the Gibbs free energy value of the transition state and G_R is the Gibbs free energy value of the reactant. The subscripts indicate the following information: “min” is for the lowest-energy conformer, “i” is for any TS conformer, and “j” is for any reactant conformer.

As Møller et al.³¹ reported, one κ value is often sufficient to represent the tunneling effect for all conformers of the same reaction. This coefficient corresponds to the tunneling coefficient of the lowest-energy TS conformer,³¹ denoted κ_min. In practice, the κ_min value was calculated using the Eckart tunneling method,⁵⁴ which requires the forward and backward barrier heights and imaginary frequency (ω_i) of the TS on the potential energy surface (PES). The PES was obtained using the intrinsic reaction coordinate (IRC)⁵⁵ approach in ORCA 6.0.1. The IRC endpoints were then optimized, and their zero-point corrected electronic energy (E_elzc) values were used to compute the barrier heights for the κ_min calculation.^31,32

2.3 Conformer sampling of free-radical reactions

For each reaction channel, a single RO˙ structure (reactant monomer) from SPARTAN was optimized at two density functional theory (DFT) levels in ORCA: first at the B3LYP^56,57/ma-def2-SVP^58,59 level and then at the ωB97X-D3^60,61/ma-def2-TZVP^58,59 level. The choice of B3LYP and ωB97X-D3 functionals was made based on previous studies on ³(RO⋯OR) complexes,^14,15,17,21 as well as benchmarking work on general chemical databases⁶² and atmospheric molecular clusters.⁶³ Tests of different basis sets on ³(AceO⋯OAce)⁶⁴ and small gas-phase clusters⁶⁵ indicate that ma-def2-SVP and ma-def2-TZVP basis sets are sufficiently accurate for the systems studied in this paper with good computational efficiency.

After optimization at the B3LYP level, the reactant monomer was optimized at the ωB97X-D3 level in two subsequent steps: first with the default grid size (no additional keywords required, or specified by the keyword “DefGrid2”)^44,66 and then with the largest predefined grid size in ORCA 6.0.1 (invoked by the keyword “DefGrid3”).^66,67 These two steps at ωB97X-D3 level will be referred to as “ωB97X-D3 (Grid2)” and “ωB97X-D3 (Grid3)” throughout this paper. Grid3 was applied to eliminate the imaginary frequencies caused by numerical noise.⁶⁸ If there were still small imaginary frequencies after optimizations with Grid3, approaches such as subtle displacement along the vector or applying tighter convergence criteria would be applied.

The TS monomer was determined by relaxed surface scans at the B3LYP level of theory. The TS structure from scans were then optimized at the B3LYP, ωB97X-D3 (Grid2) and ωB97X-D3 (Grid3) levels, similar to its corresponding reactant monomer.

Based on the optimized reactant and TS monomer structures, conformer searches were performed using GOAT with the GFN2-xTB^28,45 method in ORCA, with an energy cutoff of 15 kcal mol⁻¹.‡ Conformers generated by GOAT were then optimized at the B3LYP level, followed by filtering based on their uniqueness and relative energy values. Detailed selection criteria are:

• Any two selected conformers must differ by more than 0.01 Å in gyration radius (ΔR_g > 0.01), 0.001 Hartree in electronic energy (ΔE_el > 0.001 Hartree), and 0.1 Debye in dipole moment (Δp > 0.1 Debye).^14,49 This criterion, which is also the default setting in JKCS 2.1, will hereafter be referred to as the “uniqueness filter” or “filtered by uniqueness”.

• The Gibbs free energy value of every selected conformer must not be more than 5 kcal mol⁻¹ (ref. 14 and 31) above the global minimum (G − G_min < 5 kcal mol⁻¹). This criterion will be referred to as “G_cutoff = 5 kcal mol⁻¹” hereafter. Other energy-based filters will also be written in the form of “E_{energy type, cutoff} = x kcal mol⁻¹”.

B3LYP-optimized and filtered conformers were further optimized at the ωB97X-D3 (Grid2) level, and the optimized structures were filtered by uniqueness. Following this, the conformers were optimized at the ωB97X-D3 (Grid3) level, and then filtered by uniqueness and G_cutoff = 2 kcal mol⁻¹.

Conformers satisfying all the criteria discussed above were forwarded to the next step. To improve the accuracy of relative energy predictions, single-point energy calculations on top of the ωB97X-D3 (Grid3) optimized structures were performed using the following coupled cluster method: UHF⁶⁹-DLPNO-CCSD(T)⁷⁰-F12⁷¹/cc-pVDZ-F12,⁷² with the keywords “tightPNO”⁷³ and “tightSCF”, cc-pVDZ-F12-CABS⁷² as the auxiliary basis set, and aug⁷⁴-cc-pVDZ⁷⁵/C as the RI approximation basis set. This method was suggested by our benchmarking results and previous studies on atmospheric molecular clusters.^63,76 Coupled cluster calculations were only carried out on lowest-energy conformers (G_min), as suggested by previous studies.^14,15,31,32

Following the DLPNO-CCSD(T)-F12 calculations, energy corrections were computed as E_sp(CCSD(T)) − E_el(DFT) for the lowest-energy reactant and TS monomers, respectively, and applied to other conformers of the same type. The coupled cluster-corrected G_R and G_TS values were applied to calculate the k(CC) of free-radical reactions.

Finally, k(DFT) and k(CC) were computed using both MC-TST (eqn (1)) and LC-TST (eqn (2)) methods to evaluate the effects of coupled-cluster corrections and contributions from multiple conformers. Consequently, four k values were reported for each reaction channel.

2.4 Conformer sampling of in-complex reactions

2.4.1 Conformer sampling of reactant dimers. Based on the reactant monomers from Section 2.3, reactant dimers were generated using ABCluster 3.3, executed through JKCS2.1. Monomers were treated rigidly (i.e., with fixed internal coordinates)⁴² and placed against each other with a CHARMM36^77,78-based force field. ABCluster requires three types of input data:^41,42 monomer coordinates (in the .xyz format), atom-wise partial charges and force field parameters (ε and σ). Partial charges were calculated using the natural bond orbital (NBO) 7.0.10^79,80 program interfaced to ORCA 6.0.1. The ε and σ values were generated using the “topgen” module of ABCluster and then examined and adjusted by the authors when necessary. With the input data above, ABCluster was invoked by the following command: JKCS2_explore -pop 1280 -gen 320 -lm 100 -exploded -repeat 10. This command would guarantee that the initial number of dimers generated for each monomer combination was 1000,§ as previous studies suggested.^15,81

Reactant dimers were collected and optimized using the GFN1-xTB^27,45 method with necessary constraints. For atmospheric molecular clusters, the GFN1-xTB method was reported to have better agreement with the DFT energy values than GFN2-xTB.^81,82 In addition, GFN1-xTB showed less problems than GFN2-xTB during the geometry optimization of alkoxy dimers. GFN1-optimized structures were then filtered by uniqueness and E_el, cutoff = 15 kcal mol⁻¹.¹⁵

The filtered conformers were then optimized at the B3LYP level in ORCA. Based on our benchmarking results, applying the E_el cutoff of 5 kcal mol⁻¹ was equivalent to applying the G cutoff of 5 kcal mol⁻¹. Thus, frequency calculations were omitted at this step to reduce computational cost. After B3LYP optimization, reactant dimers were filtered by uniqueness and E_el, cutoff = 5 kcal mol⁻¹.^15,31

Subsequently, reactant dimers were optimized at the ωB97X-D3 level in ORCA. Geometry optimizations and frequency calculations were first performed with Grid2 and filtered by uniqueness. Similar calculations were then performed with Grid3, with special attention to convergence and frequencies, and the final structures were filtered by uniqueness.

The reactant dimer with the lowest G value at the ωB97X-D3 (Grid3) level was selected for the coupled cluster calculation. Its E_sp was calculated using the DLPNO-CCSD(T)-F12 method, and the energy correction (see Section 2.3) was applied to other reactant dimers. The coupled cluster-corrected G values of the reactant dimers were substituted into the G_R terms of eqn (1) and (2).

2.4.2 Conformer sampling of TS H-shift conformers. As discussed in Section 2.4.1, ABCluster 3.3 requires “isolated” monomers as input to generate the corresponding dimers. However, the TS H-shift conformer resembles more a connected moiety (two RO˙ connected through a C–H⋯O structure) rather than two separate monomers. Thus, ABCluster is not suitable for building initial TS H-shift structures. Alternatively, generation and semi-empirical optimization of TS H-shift conformers were performed using GOAT, with the GFN1-xTB method and an energy cutoff of 15 kcal mol⁻¹.

TS H-shift conformers from GOAT were then optimized using the B3LYP method. A constrained optimization without frequency calculation was first performed. The constrained optimization step fixes the distance between the two atoms strongly involved in the TS, while relaxing and optimizing the rest of the structure.⁸³ For the alkoxy systems studied, this step can reduce the computational cost and allow the following TS optimization to converge more easily. Each imaginary frequency was checked carefully to ensure it corresponded to the expected vibrational mode. Conformers that converged with correct vibrational modes at the B3LYP level were filtered by uniqueness and G_cutoff = 5 kcal mol⁻¹.

Following that, TS H-shift conformers were optimized at the ωB97X-D3 level in three steps: constrained optimization, TS optimization with Grid2, and TS optimization with Grid3. Frequencies were carefully examined, and de-noise approaches described in Section 2.3 were applied if necessary. A uniqueness filter was applied both before and after the ωB97X-D3 (Grid3) optimization.

Finally, a single point calculation of the TS H-shift conformer at the DLPNO-CCSD(T)-F12 level was performed on top of the ωB97X-D3 (Grid3) optimized structure with the lowest G value. Energy corrections were calculated as discussed in Section 2.3. The coupled cluster-corrected G values of TS H-shift conformers were then substituted into the G_TS terms of eqn (1) and (2). Together with the G_R terms from Section 2.4.1, these values were applied to calculate the k(CC) of the intermolecular H-shift reaction. In addition, the k(DFT) values were also calculated as described in Section 2.3.

2.4.3 Conformer sampling of TS dimers. For other intramolecular reaction channels of ³(RO⋯OR), their TS conformers behave more like a complex consisting of two separate RO˙ moieties, similar to the reactant dimers described in Section 2.4.1. Therefore, ABCluster can be used to generate these TS dimers with a suitable input (one reactant monomer and one TS monomer). Semi-empirical optimizations were performed using the GFN1-xTB method. After applying the uniqueness and E_el, cutoff = 15 kcal mol⁻¹ filters, TS dimers were optimized at the B3LYP level in two steps: constrained and TS optimization. Geometry convergence and frequencies were carefully examined, followed by the uniqueness and G_cutoff = 5 kcal mol⁻¹ filters.

After B3LYP optimizations and filtering, the TS dimers were optimized at the ωB97X-D3 level in three steps, as described in Section 2.4.2. A uniqueness filter was applied both before and after the ωB97X-D3 (Grid3) optimization. Frequencies were examined as mentioned in Section 2.4.2, and when required, de-noise approaches (see Section 2.3) were applied. Thereafter, the E_sp values were calculated at the DLPNO-CCSD(T)-F12 level, and energy corrections were applied as mentioned in Section 2.3. The coupled cluster-corrected G values of TS dimers were substituted into the G_TS terms of eqn (1) and (2). Together with the G_R terms from Section 2.4.1, these values were employed to compute the k(CC) of other intramolecular ³(RO⋯OR) reactions. Additionally, the k(DFT) values were also calculated as described in Section 2.3.

3 Results and discussion

This section first presents the number of conformers for each alkoxy system after sampling (Table 3) and the calculated k values using different methods (Table 4). Further discussion on each system is provided in the following subsections. Multireference and spin contamination diagnostics, conformer-wise relative energy values (including conformers from sampling and IRC endpoints) and imaginary frequencies that correspond to vibrational modes of the transition state are reported in Tables S3–S10.

Table 3 Final conformers for k calculation after ωB97X-D3 (Grid3) optimization and filtering

System	Conformer type	Number
a The total number of conformers was counted after excluding the high-energy conformers which had a small imaginary frequency that could not be eliminated (see Sections 3.3 and 3.4 for detailed discussion).
AceO	Reactant monomer	2
	TS monomer (β-scission)	2
	Reactant dimer	7
	TS dimer (β-scission)	3
	TS H-shift conformer	11
β-ISOPO	Reactant monomer	8
	TS monomer (β-scission, hydroxymethyl side)	6
	Reactant dimer	25
	TS dimer (β-scission, hydroxymethyl side)	26
PhCH2O	Reactant monomer	1
	TS monomer (epoxy formation)	1
	Epoxy radical	1
	TS monomer (epoxy decomposition)	1
	Reactant dimer	3
	TS dimer (epoxy formation)	7
	Epoxy–PhCH2O complex	12
	TS dimer (epoxy decomposition)^a	12
	TS H-shift conformer	3
PhC(O)O	Reactant monomer	1
	TS monomer (β-scission)	1
	Reactant dimer^a	9
	TS dimer (β-scission)^a	7

---

Table 4 Rate coefficients of four alkoxy systems at 298.15 K. The k values were calculated at the ωB97X-D3 (Grid3) (labeled “DFT”) and DLPNO-CCSD(T)-F12 (labeled “CC”) levels using LC-TST (labeled “LC”) and MC-TST (labeled “MC”) methods, and are reported in s⁻¹

System	Reaction type	κmin (DFT)	κmin (CC)	kLC (DFT)	kLC (CC)	kMC (DFT)	kMC (CC)
a Only the β-scission reaction from the hydroxymethyl side is reported here (Rm3.2.1 in Fig. 3, see Section 3.2 for further discussion).b The epoxy formation and decomposition reactions of PhCH2O˙ and the β-scission reaction of PhC(O)O˙ only have one conformer for each species. Thus, for these reactions kMC-TST = kLC-TST.c High-energy dimers with a small imaginary frequency were excluded from k calculations (see Sections 3.3 and 3.4 for detailed discussion).d The reactant monomer obtained from IRC showed nonphysical imaginary frequencies after DFT optimization. Instead, a reactant monomer with the same geometry (RMSD = 0.09 Å) was used to calculate κmin values.
AceO	β-Scission, free-radical	1.0	1.0	2.33 × 10^8	3.57 × 10^9	2.33 × 10^8	3.56 × 10^9
	β-Scission, in-complex	1.0	1.0	1.65 × 10^9	1.54 × 10^10	6.40 × 10^8	5.99 × 10^9
	Intermolecular H-shift	1.7	1.7	8.91 × 10^6	6.22 × 10^5	1.29 × 10^7	8.98 × 10^5
β-ISOPO^a	β-Scission, free-radical	1.2	1.2	2.91 × 10^9	7.44 × 10^8	3.70 × 10^9	9.47 × 10^8
	β-Scission, in-complex	1.1	1.1	7.49 × 10^11	6.46 × 10^11	3.73 × 10^11	3.22 × 10^11
PhCH2O	Epoxy formation, free-radical^b	1.3	1.3	1.19 × 10^3	1.59 × 10^4	1.19 × 10^3	1.59 × 10^4
	Epoxy formation, in-complex	1.2	0.9	2.93 × 10^3	2.33 × 10^4	3.13 × 10^3	2.49 × 10^4
	Epoxy decomposition, free-radical^b	1.7	1.7	1.21 × 10^6	1.29 × 10^7	1.21 × 10^6	1.29 × 10^7
	Epoxy decomposition, in-complex^c	1.7	1.7	7.19 × 10^5	9.99 × 10^6	5.40 × 10^5	7.50 × 10^6
	Intermolecular H-shift	3.5	4.0	6.79 × 10^6	2.96 × 10^6	6.30 × 10^6	2.75 × 10^6
PhC(O)O	β-Scission, free-radical^b	1.1^d	1.1^d	1.97 × 10^7	2.49 × 10^8	1.97 × 10^7	2.49 × 10^8
	β-Scission, in-complex^c	1.1	1.0	1.05 × 10^7	1.32 × 10^8	6.40 × 10^6	8.02 × 10^7

---

3.1 The AceO system

Structures of key species during the β-scission and intermolecular H-shift reactions are shown in Fig. 2. Given the relatively simple structure of AceO˙, only a small number of conformers were found for the AceO system after the sampling workflow (see Table 3). Rate coefficients calculated at different levels of theory are listed in Table 4.

Fig. 2 Reaction scheme of ³(AceO⋯OAce). All 3D structures presented here are conformers with the lowest G after the ωB97X-D3 (Grid3) optimization and filtering. Color coding: gray for C, white for H, and red for O.

Table 4 and previous theoretical studies indicate that the in-complex β-scission is among the fastest of the four reaction channels of ³(AceO⋯OAce): its rate coefficient (10⁹–10¹⁰ s⁻¹) exceeds that of the intermolecular H-shift (10⁵–10⁷ s⁻¹, this study) and ISC (4 × 10⁸ s⁻¹ (ref. 14 and 16)) reactions, and is comparable to that of dissociation (1.5 × 10⁹ s⁻¹ (ref. 13)). Moreover, for all four methods employed to compute k values, the presence of another AceO˙ in ³(AceO⋯OAce) consistently increases the β-scission rate (see the k values in the same column of Table 4).

However, in contrast to the high theoretical k values, accretion products corresponding to the β-scission reaction (acetonyl acetate, C₅H₈O₃) were only observed at concentrations close to the detection limit in a recent experimental study.¹⁸ Meanwhile, peroxide products (C₆H₁₀O₄) corresponding to ISC were more abundant. The discrepancy between theoretical and experimental results implies that competing processes, rather than direct recombination on the triplet surface, may take place after the β-scission reaction. Our preliminary calculations on the ³(CH₃C(O)⋯OAce) complex¶ indicate that although recombination on the triplet surface is thermodynamically feasible, it is prevented by a substantial barrier (see Table S11). Thus, the formation of ROR′-type products likely requires an ISC after β-scission.

For the intermolecular H-shift reaction, our k(DFT) values agree well with the reported value (8.85 × 10⁶ s⁻¹ (ref. 14)), while our k(CC) values are two orders of magnitude higher than the reported value (1.62 × 10³ s⁻¹ (ref. 14)). The lowest-energy TS H-shift conformer found by our conformer sampling workflow has lower E_elzc and G than that reported by Hasan et al.¹⁴ (see Table S12), suggesting that part of this variation may arise from different geometries used for CCSD(T) calculations. In addition, subtle differences in the CCSD(T) methods employed by the two studies may also contribute.

3.2 The β-ISOPO system

The increased structural complexity of the β-ISOPO system leads to several possible pathways: β-scissions from the hydroxymethyl (R_m3.2.1 in Fig. 3), methyl (R_m3.2.2 in Fig. 3) and ethenyl (R_m3.2.3 in Fig. 3) sides; and radical-double bond addition reactions at the β position (R3.2.2 (a) in Fig. 3) and the γ position (R3.2.2 (b) in Fig. 3). Intermolecular H-shift reactions are unlikely to occur, as no “aldehyde + alcohol” products can be formed for ³(β-ISOPO⋯β-ISOPO).

Fig. 3 Reaction scheme of the β-ISOPO system. All 3D structures presented are the lowest-G conformers, with the same level of theory and color coding as in Fig. 2.

Given the expected large number of reaction channels and conformers of the β-ISOPO system, we first calculated the energy barriers (Δ_bE_elzc and Δ_bG, as defined in Table 2) of the following reactions, prior to comprehensive conformer sampling: three β-scission channels of isolated β-ISOPO˙ (R_m3.2.1–3.2.3) and two radical-double bond addition channels of ³(β-ISOPO⋯β-ISOPO) (R3.2.2 (a) and (b)). This step aimed to determine the energetically most favorable pathway, which was then selected for systematic conformer sampling. The barrier heights are listed in Table S13.

As shown in Table S13, R_m3.2.1 has significantly lower energy barriers than the other channels. The barrier heights of the three β-scission channels align with the trend of unimolecular RO˙ reactions predicted by the structure–activity relationship (SAR) models¹⁹ (Δ_bE_elzc(–CH₂OH) < Δ_bE_elzc(–CH₃) < Δ_bE_elzc(–CH=CH₂)). Furthermore, differences among Δ_bE_elzc and Δ_bG values are large enough that the presence of another β-ISOPO˙ is unlikely to change this conclusion. Thus, the systematic conformer sampling workflow was applied only to β-scission reactions from the hydroxymethyl side (R_m3.2.1 and R3.2.1 in Fig. 3). The number of conformers and calculated k values are presented in Tables 3 and 4, respectively.

For R_m3.2.1, our k(CC) values (2.33 × 10⁸ s⁻¹) are around 10 times larger than the reported values (10⁷ s⁻¹ at 300 K and 760 Torr).^84,85 This discrepancy may result from different computational methods and kinetic models between earlier studies and this study (e.g. B3LYP^84,85 and ωB97X-D3, CCSD(T)^84,85 and DLPNO-CCSD(T)-F12, and RRKM⁸⁵ and elementary TST, where in each pair the latter method is used in this study).

According to the first two sections in Table 4, the presence of another RO˙ in ³(RO⋯OR) has a stronger impact on β-scission rates for the β-ISOPO system than for AceO. For β-ISOPO, the k(in-complex) values are over two orders of magnitude higher than k(free-radical) values. We note that the Δ_bG of R3.2.1 is only 1.30 kcal mol⁻¹ at the ωB97X-D3 (Grid3) level, implying that elementary TST may not provide accurate k values.⁸⁶ Nevertheless, the notable difference in barrier heights (Δ_bG(free-radical, DFT) − Δ_bG(in-complex, DFT) = 3.34 kcal mol⁻¹) supports our qualitative conclusion: the in-complex β-scission reaction of β-ISOPO proceeds much faster than its free-radical counterpart.

3.3 The PhCH₂O system

The β-scission reaction of isolated PhCH₂O˙ has large energy barriers (Δ_bE_elzc = 25.98 kcal mol⁻¹ and Δ_bG = 25.35 kcal mol⁻¹) at the ωB97X-D3 (Grid3) level. Since the difference in Δ_bE values between in-complex and free-radical reactions is typically less than 5 kcal mol⁻¹ (see Section S5 in the SI), the β-scission channel is unlikely to be significant for ³(PhCH₂O⋯OCH₂Ph) as well. We then examined other possible free-radical reactions of PhCH₂O˙ and proposed an epoxy formation–decomposition pathway (R3.3.2(1) and (2) in Fig. 4). Similar epoxy formation–decomposition processes have been reported for unsaturated aliphatic RO˙.^22,87 The conformer sampling workflow mentioned in Section 2 was applied to the species involved in epoxy formation–decomposition and intermolecular H-shift reactions.

Fig. 4 Reaction scheme of the PhCH₂O system. All 3D structures presented are the lowest-G conformers, with the same level of theory and color coding as in Fig. 2.

Each monomer of the PhCH₂O system (see Table 3) has only one conformer because of the rigid benzene ring. Thus, k(MC-TST) equals k(LC-TST) for free-radical epoxy formation and decomposition reactions. Besides, one high-energy TS dimer of the epoxy decomposition reaction showed a second small imaginary frequency (−8.39 cm⁻¹), which could not be eliminated by the Grid3 optimization or manual displacements as described in Section 2.3. This conformer was therefore excluded from the k calculation.

Table 4 summarizes the k values for the free-radical and in-complex reactions of the PhCH₂O system. Compared to earlier studies,^22,87 the epoxy formation reaction shows a larger energy barrier in isolated PhCH₂O˙ (Δ_bE_elzc(CC, this study) = 12.86 kcal mol⁻¹) than in aliphatic RO˙ (Δ_bE(coupled cluster, earlier studies) = 4.60²²–9.5⁸⁷ kcal mol⁻¹). This gap is significant regardless of the different coupled cluster methods applied. We further examined their reaction energies (Δ_rE, defined in Table 2). The epoxy formation reaction in aliphatic RO˙ is isoergic or slightly endoergic (Δ_rE(coupled cluster, earlier studies) = −0.66–2.23 kcal mol⁻¹ (ref. 22)), while in PhCH₂O˙ it is much more endoergic (Δ_rE_elzc(CC) = 10.27 kcal mol⁻¹) and, therefore, thermodynamically less favorable.

Although the epoxy decomposition channel has k values comparable to those of intermolecular H-shift, its contribution is not expected to exceed 1% of the final products. This is because its prerequisite step, epoxy formation, proceeds around two orders of magnitude more slowly than the intermolecular H-shift (see Table 4).

Unlike the two systems mentioned in Sections 3.1 and 3.2, the presence of another PhCH₂O˙ does not exhibit a consistent effect on the reactions of ³(PhCH₂O⋯OCH₂Ph): it slightly increases the k values of epoxy formation, but decreases those of epoxy decomposition (see Table 4).

3.4 The PhC(O)O system

Due to the structural constraints of PhC(O)O˙, only the β-scission channel (see Fig. 5) is feasible for this system, and the number of conformers is relatively small (see Table 3). Similar to the PhCH₂O system, one high-energy reactant dimer and two high-energy TS dimers of the β-scission in ³(PhC(O)O⋯O(O)Ph) showed small imaginary frequencies that could not be eliminated after using the techniques mentioned in Section 2.3. Thus, these three conformers were not included when calculating k.

Fig. 5 Reaction scheme of the PhC(O)O system. All 3D structures presented are the lowest-G conformers, with the same level of theory and color coding as in Fig. 2.

As shown in Table 4, free-radical and in-complex β-scission reactions are relatively fast for the PhC(O)O system. The presence of another PhC(O)O˙ in the complex increases k(CC) while decreasing k(DFT). However, based on previous studies on CH₃C(O)O˙,^88,89 the CO₂ ejection transition state likely has substantial multireference characters, suggesting that neither k(DFT) nor k(CC) is particularly accurate.

3.5 General findings of the four systems

We plotted the range of k for β-scission reactions in AceO, β-ISOPO and PhC(O)O systems, as shown in Fig. 6. The β-scission reaction is kinetically unfeasible for the PhCH₂O system and is thus not included in the figure.

Fig. 6 The range of free-radical and in-complex β-scission reactions of different systems. The origin of the x-axis is set to 10⁶ s⁻¹. Acronyms: MC for MC-TST, LC for LC-TST, DFT for ωB97X-D3 (Grid3), and CC for DLPNO-CCSD(T)-F12.

From Fig. 6, the in-complex β-scission reaction roughly follows the trend k(β-ISOPO) > k(AceO) > k(PhC(O)O), which is consistent with the SAR-predicted order of unimolecular RO˙ reactions.^19,87 While Fig. 6 shows the overall magnitude of k for β-scission reactions, we also evaluate factors that influence the calculated k values of all studied reaction channels. These factors include: (1) the reaction takes place in ³(RO⋯OR) or in isolated RO˙. (2) k was calculated using the MC-TST or LC-TST method. (3) The Δ_bG values were calculated at the ωB97X-D3 (Grid3) or DLPNO-CCSD(T)-F12 level of theory.

To quantify the effects of those factors on k, we defined three parameters: p(d/m), p(MC/LC) and p(CC/DFT). For computing p(d/m) (eqn (3)) and p(CC/DFT) (eqn (5)), only k(LC-TST) values are used to eliminate the effect of multiple conformers. Since intermolecular H-shift reactions involve two RO˙ in the ³(RO⋯OR) complex, p(d/m) does not apply to them. p(MC/LC) reflects the effects of multiple conformers and is computed using only k(in-complex) values. For some free-radical reactions, there is only one conformer for each species (see Table 3), and consequently, their p(MC/LC) always equals 1. The calculated p values are listed in Table 5.

Table 5 Results of p(d/m), p(MC/LC) and p(CC/DFT). Parameters were calculated using eqn (3), (4) and (5), respectively. Definitions of k(DFT) and k(CC) are explained in Table 2

System and reaction	p(d/m)	p(MC/LC)	p(CC/DFT)
AceO, β-scission	4.33	0.39	9.35
AceO, intermolecular H-shift	—	1.44	0.07
β-ISOPO, β-scission	868.81	0.50	0.86
PhCH2O, epoxy formation	1.46	1.07	7.96
PhCH2O, epoxy decomposition	0.77	0.75	13.89
PhCH2O, intermolecular H-shift	—	0.93	0.44
PhC(O)O, β-scission	0.53	0.61	12.52

---

According to Tables 4 and 5, the formation of ³(RO⋯OR) significantly increases the β-scission rates of aliphatic systems. In particular, the presence of the other β-ISOPO˙ raises the k(β-scission) by a factor of 868 at the DLPNO-CCSD(T)-F12 level. This increase possibly results from an intermolecular hydrogen bond that stabilizes the TS dimer (O₃–H₁₂⋯O₁ in Fig. S2b). Although the TS monomer also has a hydrogen bond (O₁–H₃⋯O₂ in Fig. S2a), it is much weaker compared to its counterpart in the TS dimer: the O⋯H distance is longer (2.4 Å vs. 1.7 Å in the TS dimer), and the O–H⋯O angle is less favorable (103° vs. 177° in the TS dimer).

In addition, the binding energy of ³(RO⋯OR) may also affect the p(d/m) values. For example, ³(β-ISOPO⋯β-ISOPO) is among the strongest-bound complexes (D_elzc(DFT) = 7.96 kcal mol⁻¹ and D_elzc(CC) = 6.52 kcal mol⁻¹; see Table S2), and its p(d/m) is exceptionally large. However, the other three complexes have similar binding energies (within 6–8 kcal mol⁻¹, see Table S2), while their p(d/m) values are between 0.5 and 5. This observation suggests that while the binding energy may contribute to differences between k(in-complex) and k(free-radical), it is unlikely to be the determining factor.

p(MC/LC) depends on the relative number of low-energy TS and reactant conformers (see Table 3). For example, the intermolecular H-shift reaction of ³(AceO⋯OAce) has 11 low-energy TS conformers, leading to p(MC/LC) > 1. However, all p(MC/LC) values are between 0.3 and 1.5, indicating that multiple conformers change the k value by less than a factor of three. Thus, the multi-conformer effects are a relatively minor source of error, as long as the global minimum conformers (with the lowest G) have been correctly identified by the sampling workflow.

p(CC/DFT) is an indicator of consistency between k(DFT) and k(CC) values. If the differences between the DFT and CCSD(T) results were consistent and systematic among the four systems, the four p(CC/DFT) values could be expected to be similar. However, p(CC/DFT) varies from 0.07 to 13.89 (see Table 5). Notably, the two H-shift reactions have p(CC/DFT) < 1, while the β-scission reactions have p(CC/DFT) values close to or above one. Considering the typical uncertainties when applying DFT and CCSD(T) methods to complex open-shell systems, the disagreement shown in Table 5 is within an acceptable range.

4 Conclusions

Based on previous studies,^14,15,31 we developed a systematic conformer sampling workflow for the key intermediate in RO₂˙ recombination reactions, ³(RO⋯OR). This workflow was applied to four representative alkoxy systems in the atmosphere: AceO, β-ISOPO, PhCH₂O and PhC(O)O. The rate coefficients for typical reactions were computed at the ωB97X-D3/ma-def2-TZVP (Grid3) and UHF-DLPNO-CCSD(T)-F12/cc-pVDZ-F12 levels of theory.

We examined the established channels of these systems, including β-scission and intermolecular H-shift reactions. Our workflow improved the conformer sampling, especially for intermolecular H-shift reactions. For example, we identified a lower-energy TS dimer for the β-scission of ³(AceO⋯OAce), compared to the best TS conformer reported previously.¹⁴ Furthermore, we investigated two novel pathways: RO˙ addition to the double bond of the other RO˙ in the complex (R3.2.2(a) and (b) in Fig. 3) and epoxy formation–decomposition of PhCH₂O˙, leading to a phenoxymethyl radical (R3.3.2(1) and (2) in Fig. 4). However, neither novel pathway was competitive compared to in-complex β-scission and intermolecular H-shift reactions (see Table 4).

Our conformer sampling workflow was tested and found to be robust for both aliphatic and aromatic systems. Including multiple conformers in k calculations has a modest effect (less than a factor of three, see Table 5), provided that the lowest-energy conformers have been correctly identified. On the other hand, CCSD(T) energy corrections influence k values by less than a factor of 15 (see Table 5). Interestingly, our CCSD(T) corrections tend to slow down intermolecular H-shift reactions, while accelerating β-scission reactions (see Table 5).

For β-scission reactions in the other three alkoxy systems, the k(in-complex) and k(free-radical) values are comparable. However, in the β-ISOPO system, the k(in-complex) is two orders of magnitude higher than k(free-radical). This rate enhancement is likely due to the intermolecular hydrogen bond that stabilizes the transition state (see Fig. S2).

Our observations on k values tentatively suggest that, for β-scission reactions, k(free-radical) can constrain the lower limit of k(in-complex), as well as estimate the order-of-magnitude for k(in-complex) in relatively simple or weakly-bound systems. However, for aerosol-relevant systems with H-bonding functional groups, k(in-complex) may be significantly enhanced. This may help explain the substantial ester formation detected in α-pinene ozonolysis.¹⁷

Author contributions

Hongye Fraise Zhao: investigation (main), methodology (main), formal analysis, and writing – original draft. Lauri Franzon: investigation (support), methodology (support), supervision (support), and writing – review & editing (support). Severi Juttula: software (support) and writing – review & editing (support). Robert Skog: methodology (support) and writing – review & editing (support). Theo Kurtén: project administration, funding acquisition, resources, supervision (main), and writing – review & editing (main). Nanna Myllys: funding acquisition, supervision (support), and writing – review & editing (support).

Conflicts of interest

There are no conflicts to declare.

Data availability

Supplementary information (SI) contains data supporting the findings in this article. Supplementary information: detailed ORCA output files (in .out format) and ωB97X-D3 (Grid3)-optimized structures (in the .xyz format) of all conformers, which contain necessary information (e.g. frequencies and energy values) for k calculations, are published in the Zenodo archive [DOI: https://doi.org/10.5281/zenodo.19210535]. See DOI: https://doi.org/10.1039/d5cp04635a.

Acknowledgements

We gratefully acknowledge the Oskar Öflunds Stiftelse for awarding the 2025 scientific research grant that supported this project. We also acknowledge funding from the Research Council of Finland (Center of Excellence Virtual Laboratory for Molecular Level Atmospheric Transformations (VILMA) (grant no. 346369) and the Academy Research Fellowship (grant no. 347775)) and the European Union (ERC project AtmosCOC no. 101161989). In addition, we acknowledge the CSC-IT Center for Science Ltd (CSC) of Finland for providing the high-performance computing (HPC) environment and technical assistance. We also acknowledge the support from the Doctoral Program in Chemistry and Molecular Sciences (CHEMS-DP) at the University of Helsinki.

Notes and references

1. K. Vohra, A. Vodonos, J. Schwartz, E. A. Marais, M. P. Sulprizio and L. J. Mickley, Environ. Res., 2021, 195, 110754.
2. T. Muenzel, M. Sørensen, O. Hahad, M. Nieuwenhuijsen and A. Daiber, Nat. Rev. Cardiol., 2023, 20, 651–669.
3. U. Pöschl, Angew. Chem., Int. Ed., 2005, 44, 7520–7540.
4. A. P. Tsimpidi, S. Scholz, A. Milousis, N. Mihalopoulos and V. A. Karydis, Atmos. Chem. Phys., 2025, 25, 10183–10213.
5. H. Wang, X. Liu, C. Wu and G. Lin, Atmos. Chem. Phys., 2024, 24, 3309–3328.
6. Z. He, X. Wang, Z. Ling, J. Zhao, H. Guo, M. Shao and Z. Wang, Atmos. Chem. Phys., 2019, 19, 8801–8816.
7. J. J. Orlando and G. S. Tyndall, Chem. Soc. Rev., 2012, 41, 6294–6317.
8. K. C. Barsanti and J. F. Pankow, Atmos. Environ., 2004, 38, 4371–4382.
9. J. H. Kroll and J. H. Seinfeld, Atmos. Environ., 2008, 42, 3593–3624.
10. S. E. Murphy, J. D. Crounse, K. H. Møller, S. P. Rezgui, N. J. Hafeman, J. Park, H. G. Kjaergaard, B. M. Stoltz and P. O. Wennberg, Environ. Sci.: Atmos., 2023, 3, 882–893.
11. V.-T. Salo, R. Valiev, S. Lehtola and T. Kurtén, J. Phys. Chem. A, 2022, 126, 4046–4056.
12. G. Ghigo, A. Maranzana and G. Tonachini, J. Chem. Phys., 2003, 118, 10575–10583.
13. L. Franzon, J. Phys. Chem. A, 2023, 127, 5956–5966.
14. G. Hasan, V.-T. Salo, R. R. Valiev, J. Kubečka and T. Kurtén, J. Phys. Chem. A, 2020, 124, 8305–8320.
15. G. Hasan, V.-T. Salo, T. Golin Almeida, R. R. Valiev and T. Kurtén, J. Phys. Chem. A, 2023, 127, 1686–1696.
16. R. R. Valiev, G. Hasan, V.-T. Salo, J. Kubečka and T. Kurtén, J. Phys. Chem. A, 2019, 123, 6596–6604.
17. O. Peräkylä, T. Berndt, L. Franzon, G. Hasan, M. Meder, R. R. Valiev, C. D. Daub, J. G. Varelas, F. M. Geiger and R. J. Thomson, et al., J. Am. Chem. Soc., 2023, 145, 7780–7790.
18. B. N. Frandsen, L. Franzon, M. Meder, D. Pasik, E. Ahongshangbam, N. Vinkvist, N. Myllys, S. Iyer, M. P. Rissanen and M. Ehn, et al., ACS Earth Space Chem., 2025, 9(6), 1322–1337.
19. L. Vereecken and J. Peeters, Phys. Chem. Chem. Phys., 2009, 11, 9062–9074.
20. L. Vereecken and J. Peeters, Phys. Chem. Chem. Phys., 2010, 12, 12608–12620.
21. L. Franzon, M. Camredon, R. Valorso, B. Aumont and T. Kurtén, Atmos. Chem. Phys., 2024, 24, 11679–11699.
22. V. S. Nguyen and J. Peeters, J. Phys. Chem. A, 2015, 119, 7270–7276.
23. L. Vereecken, P. Carlsson, A. Novelli, F. Bernard, S. Brown, C. Cho, J. Crowley, H. Fuchs, W. Mellouki and D. Reimer, et al., Phys. Chem. Chem. Phys., 2021, 23, 5496–5515.
24. J. J. Orlando, G. S. Tyndall and T. J. Wallington, Chem. Rev., 2003, 103, 4657–4690.
25. I. Mandal, C. D. Daub, R. Valiev, T. Kurtén and R. B. Gerber, Phys. Chem. Chem. Phys., 2025, 27, 2395–2404.
26. B. Nozière and R. Patrick, Sci. Adv., 2026, 12, eaeb6495.
27. S. Grimme, C. Bannwarth and P. Shushkov, J. Chem. Theory Comput., 2017, 13, 1989–2009.
28. C. Bannwarth, S. Ehlert and S. Grimme, J. Chem. Theory Comput., 2019, 15, 1652–1671.
29. S. Spicher and S. Grimme, Angew. Chem., Int. Ed., 2020, 59, 15665–15673.
30. International Union of Pure and Applied Chemistry (IUPAC), conformer, 2025 DOI:10.1351/goldbook.C01262.
31. K. H. Møller, R. V. Otkjær, N. Hyttinen, T. Kurtén and H. G. Kjaergaard, J. Phys. Chem. A, 2016, 120, 10072–10087.
32. Q. Zhao, K. H. Møller, J. Chen and H. G. Kjaergaard, J. Phys. Chem. A, 2022, 126, 6483–6494.
33. D. Ferro-Costas, E. Martinez-Nunez, J. Rodriguez-Otero, E. Cabaleiro-Lago, C. M. Estévez, B. Fernandez, A. Fernandez-Ramos and S. A. Vazquez, J. Phys. Chem. A, 2018, 122, 4790–4800.
34. T. Berndt, W. Scholz, B. Mentler, L. Fischer, H. Herrmann, M. Kulmala and A. Hansel, Angew. Chem., Int. Ed., 2018, 57, 3820–3824.
35. K. Zuraski, F. J. Grieman, A. O. Hui, J. Cowen, F. A. Winiberg, C. J. Percival, M. Okumura and S. P. Sander, J. Phys. Chem. A, 2023, 127, 7772–7792.
36. P. O. Wennberg, K. H. Bates, J. D. Crounse, L. G. Dodson, R. C. McVay, L. A. Mertens, T. B. Nguyen, E. Praske, R. H. Schwantes and M. D. Smarte, et al., Chem. Rev., 2018, 118, 3337–3390.
37. R. Wu, S. Pan, Y. Li and L. Wang, J. Phys. Chem. A, 2014, 118, 4533–4547.
38. B. Feng, C. Sun and S. Zhang, Atmos. Environ., 2019, 201, 18–27.
39. Y. Ren, L. Zhou, A. Mellouki, V. Daële, M. Idir, S. S. Brown, B. Ruscic, R. S. Paton, M. R. McGillen and A. R. Ravishankara, Atmos. Chem. Phys., 2021, 21, 13537–13551.
40. Wavefunction Inc., SPARTAN'24, https://www.wavefun.com/spartan, 2024, Accessed August 2025.
41. J. Zhang and M. Dolg, Phys. Chem. Chem. Phys., 2015, 17, 24173–24181.
42. J. Zhang and M. Dolg, Phys. Chem. Chem. Phys., 2016, 18, 3003–3010.
43. B. de Souza, Angew. Chem., Int. Ed., 2025, 64, e202500393.
44. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2025, 15, e70019.
45. C. Bannwarth, E. Caldeweyher, S. Ehlert, A. Hansen, P. Pracht, J. Seibert, S. Spicher and S. Grimme, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2021, 11, e1493.
46. S. Grimme et al., XTB version 6.7.1, GitHub repository, 2024, https://github.com/grimme-lab/xtb/releases/tag/v6.7.1, Accessed August 2025.
47. J. Kubečka, V. Besel, T. Kurtén, N. Myllys and H. Vehkamäki, J. Phys. Chem. A, 2019, 123, 6022–6033.
48. J. Kubečka, A. S. Christensen, F. R. Rasmussen and J. Elm, Environ. Sci. Technol. Lett., 2022, 9, 239–244.
49. J. Kubečka, JKCS2.1: Jammy Key for Configurational Sampling, https://github.com/kubeckaj/JKCS2.1, 2024, Accessed August 2025.
50. International Union of Pure and Applied Chemistry (IUPAC), elementary reaction, https://doi.org/10.1351/goldbook.E02035, 2025.
51. H. Eyring, J. Chem. Phys., 1935, 3, 107–115.
52. K. J. Laidler and M. C. King, J. Phys. Chem., 1983, 87, 2657–2664.
53. L. Vereecken and J. Peeters, J. Chem. Phys., 2003, 119, 5159–5170.
54. C. Eckart, Phys. Rev., 1930, 35, 1303.
55. K. Ishida, K. Morokuma and A. Komornicki, J. Chem. Phys., 1977, 66, 2153–2156.
56. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
57. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785.
58. J. Zheng, X. Xu and D. G. Truhlar, Theor. Chem. Acc., 2011, 128, 295–305.
59. F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys., 2005, 7, 3297–3305.
60. J.-D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys., 2008, 10, 6615–6620.
61. Y.-S. Lin, G.-D. Li, S.-P. Mao and J.-D. Chai, J. Chem. Theory Comput., 2013, 9, 263–272.
62. L. Goerigk, A. Hansen, C. Bauer, S. Ehrlich, A. Najibi and S. Grimme, Phys. Chem. Chem. Phys., 2017, 19, 32184–32215.
63. J. Elm, J. Kubečka, V. Besel, M. J. Jääskeläinen, R. Halonen, T. Kurtén and H. Vehkamäki, J. Aerosol Sci., 2020, 149, 105621.
64. H. Zhao, Rapid alkoxy radical β-scission in bimolecular complexes, 2024, https://hdl.handle.net/10138/589458.
65. N. Myllys, J. Elm and T. Kurtén, Comput. Theor. Chem., 2016, 1098, 1–12.
66. B. Helmich-Paris, B. de Souza, F. Neese and R. Izsák, J. Chem. Phys., 2021, 155, 104109.
67. F. Neese et al., ORCA 6.0 Manual Section 7.3. Details on the numerical integration grids, 2024, https://www.faccts.de/docs/orca/6.0/manual/contents/detailed/numint.html, Accessed November 2025.
68. S. P. Sitkiewicz, R. Zaleśny, E. Ramos-Cordoba, J. M. Luis and E. Matito, J. Phys. Chem. Lett., 2022, 13, 5963–5968.
69. J. A. Pople and R. K. Nesbet, J. Chem. Phys., 1954, 22, 571–572.
70. C. Riplinger, B. Sandhoefer, A. Hansen and F. Neese, J. Chem. Phys., 2013, 139, 134101.
71. F. Pavošević, F. Neese and E. F. Valeev, J. Chem. Phys., 2014, 141, 054106.
72. K. A. Peterson, T. B. Adler and H.-J. Werner, J. Chem. Phys., 2008, 128, 084102.
73. D. G. Liakos, M. Sparta, M. K. Kesharwani, J. M. Martin and F. Neese, J. Chem. Theory Comput., 2015, 11, 1525–1539.
74. R. A. Kendall, T. H. Dunning Jr and R. J. Harrison, J. Chem. Phys., 1992, 96, 6796–6806.
75. T. H. Dunning Jr, J. Chem. Phys., 1989, 90, 1007–1023.
76. G. Schmitz and J. Elm, ACS Omega, 2020, 5, 7601–7612.
77. R. B. Best, X. Zhu, J. Shim, P. E. Lopes, J. Mittal, M. Feig and A. D. MacKerell Jr, J. Chem. Theory Comput., 2012, 8, 3257–3273.
78. A. D. M. Jr. et al., CHARMM Force Fields, https://mackerell.umaryland.edu/charmm_ff.shtml, 2025, Accessed August 2025.
79. E. D. Glendening, J. K. Badenhoop, A. E. Reed, J. E. Carpenter, J. A. Bohmann, C. M. Morales, P. Karafiloglou, C. R. Landis and F. Weinhold, NBO 7.0, Version 7.0.10, Theoretical Chemistry Institute, University of Wisconsin, Madison, WI, 2018.
80. E. D. Glendening, C. R. Landis and F. Weinhold, J. Comput. Chem., 2019, 40, 2234–2241.
81. H. Wu, M. Engsvang, Y. Knattrup, J. Kubečka and J. Elm, ACS Omega, 2023, 8, 45065–45077.
82. Y. Knattrup, J. Kubečka, H. Wu, F. Jensen and J. Elm, RSC Adv., 2024, 14, 20048–20055.
83. F. Neese et al., ORCA 6.0 Manual Section 6.3.6: Constrained Optimizations, 2024, https://www.faccts.de/docs/orca/6.0/manual/contents/typical/optimizations.html#constrained-optimizations, Accessed August 2025.
84. W. Lei and R. Zhang, J. Phys. Chem. A, 2001, 105, 3808–3815.
85. J. Park, J. C. Stephens, R. Zhang and S. W. North, J. Phys. Chem. A, 2003, 107, 6408–6414.
86. J. L. Bao and D. G. Truhlar, Chem. Soc. Rev., 2017, 46, 7548–7596.
87. A. Novelli, C. Cho, H. Fuchs, A. Hofzumahaus, F. Rohrer, R. Tillmann, A. Kiendler-Scharr, A. Wahner and L. Vereecken, Phys. Chem. Chem. Phys., 2021, 23, 5474–5495.
88. C. M. Rosado-Reyes and J. S. Francisco, J. Phys. Chem. A, 2006, 110, 4419–4433.
89. Y. Z. Zhou, S. Li, Q. S. Li and S. W. Zhang, J. Mol. Struct. THEOCHEM, 2008, 854, 40–45.

---

Footnotes

† According to the experimental results of Murphy et al.,¹⁰ for certain peroxy radicals such as HOCH₂CH₂OO˙, the recombination reaction might also proceed through a ¹(RO⋯O₂⋯OR) complex and produce HO₂˙. The reaction can be summarized as follows: RO₂˙ + RO₂˙ → ¹(RO⋯O₂⋯OR) → R_–H=O + HO₂˙ + RO˙.

‡ The conformer search will stop going uphill if E(new conformer) − E(initial conformer) > 15 + 3N_atom kcal mol⁻¹, where N_atom is the number of atoms in the conformer. This is invoked by the keyword “MaxEn 15.0”.⁴³

§ For example, if there are two conformers of AceO˙, the number of possible ³(AceO⋯OAce) combinations will be three, and the ABCluster step will generate 3000 initial structures of ³(AceO⋯OAce).

¶ This complex corresponds to the product complex of the β-scission reaction after ejecting the formaldehyde. The processes can be summarized as follows: ³(AceO⋯OAce) → ³(CH₃C(O)⋯CH₂O⋯OAce) → ³(CH₃C(O)⋯OAce) + CH₂O.

---