Quantum Dynamics Studies of Anomalous Isotope Effects

1. Quantum Dynamics Studies of Anomalous Isotope Effects Dmitri Babikov Marquette University, Chemistry Department Milwaukee, Wisconsin, USA.

2. Outline What we learned about the MIF from studies of the anomalous isotope effect in O3? What remains a challenge for theory and experiment on ozone? What has already been done on the S-containing species? What can be done in the near future on the S-containing species?

3. Traced to the three-body recombination reaction which forms ozone: O + O2 + M  O3 + M Rates can differ by more then 50%, remarkably large isotope effect taking into account a small change of mass! (Mauersberger and co., J. Geophys. Res.95, D1, 901, 1990). Discovery of Enrichments in O3 In the atmosphere: “anomalous” In the laboratory studies: equal in 17O and 18O “mass independent” (Heidenreich & Thiemens, JCP78, 892, 1983) (Mauersberger, Geophys. Res. Lett.8, 935, 1981)

4. This is not a complete list… First Round of Theoretical Research on Anomalous Isotope Effect

5. State of the Problem before the Y 2000 “Currently, in spite of 15 years of intensive experimental and theoretical investigation,.. the mechanism that is responsible remains unidentified.” “Despite the progress that has been made during the past 10 years, a convincing physical explanation of the process that results in enrichment is still missing.” M. H. Thiemens, Science 283, 341 (1999) K. Mauersberger, Science 283, 370 (1999) • Theory was not particularly successful. • - Experiments results were not complete and were not summarized in the form useful to theoreticians. • - Situation changed dramatically in the XXI century!

6. Very strong isotope effects. No any obvious correlation with O3 masses. Relative Rates (Exp.) (Mauersberger and co., PCCP3, 4718, 2001)

7. 18→18 17→17 17→18 18→17 17→16 16→16 16→18 18→16 16→17 (Mauersberger and co., PCCP3, 4718, 2001) Experimental Rates vs.DZPE xO + yOzO → (xOyOzO)* →xOyO + zO + M → O3 Their explanation of the correlation: ZPEs of O2 on the reactant and product sides are different. The atom exchange reaction is slightly exothermic or endothermic. This affects rate of the reaction, and lifetime of the intermediates. This paper stimulated new round of theoretical work.

8. Second Round of Theoretical Research on Anomalous Isotope Effect

9. Statistical Theory of Anomalous Isotope Effect (Marcus) Nb† Na† RRKM • Ya,b – partitioning factor ( due to DZPE ); • h – symmetry factor ( in either w or r ). 1618+18 16+1818 O3 - Densities of states are computed from models (free rotor; hindered rotor). - Simple models for deactivation are assumed (strong collision; step-ladder; exponential). h =1.18 DE =210 cm-1 - Parameters DEand h are tuned to fit expe- rimental results for 16+1818 and 18+1616. Findings: - Effects of the Ya,b cancel in the “scrambled” conditions and have no effect on enrichments. - The h-effect is essential for enrichments. h =1 Marcus and co., JCP 117, 1536 (2002).

10. The First Quantum Dynamics Treatment (Clary) ● Vibrational motion of O3 is described by the wave function: ● A reduced dimensionality approximation is employed in which the bending angle in Jackoby coordinates is fixed: This wave function is represented numerically by a 2D-grid of points (140x140) ● Bound states and scattering resonances (metastable states) are computed for J = 0 by solving the TISE using the stabilization method: ~ 90 bound states resonances O3* ~ 60 resonances Charlo and Clary, JCP 117, 1660 (2002). O3 bound states

11. The Coupled-Channel Approach The collision of O3+ M is treated by introducing spherical polar coordinates for M: The scattering of M is described by a multi-dimensional wave function in the form: Sudden collision approximation is used (O3 cant rotate during the collision): • dependence is parametric (essentially 1D). Separate calculations are performed for a large number of fixed values of angles (10x24). Final results are obtained by averaging over the angles. Next approximation neglects the ro-vibrational couplings for O3 states (IOS). same as in O3

12. The Coupled-Channel Approach Substitution into the TISE leads to a system of coupled-channel equations: - wave vector in channel v - potential coupling matrix ~ 100 open, ~ 20 closed channels. Solved on a grid of points in RM (~150): Result of the calculations is a scattering matrix for state-to-state transitions: Done for ~ 30 total energies E, separately for ℓ = 0 to 150:

13. The First Quantum Results (Clary) ● Results of this very nice work were not particularly encouraging: only one reaction showed large isotope effect in the right direction. In all other cases the isotope effect was either small, or in wrong direction. ● No any correlation was observed. No clear mechanism proposed. ● Due to a number of approximations used it was somewhat hard to figure out the problem. Note: This method was successful in the treatment of HCN + Ar energy transfer. What possibly could be a problem in O3 case? My guess is: Unlucky combination of the PES features + reduced dimensionality implemented using Jackobi coordinates.

14. Role of Quantum Scattering Resonances (Babikov) ● Focus on energies and lifetimes of O3* states. ● Vibrational wave functions are full dimensional (3D): - the adiabatically-adjusting principal-axes hyper-spherical coordinates (APH) are employed. ● Accurate ab initio global PES was used. ● Only the J = 0 case was considered. ● Coupled-Channel treatment of O + O2 scattering is employed. O + O2 O3* O3*+ M  O3 + M

15. 18O18O 16O18O 16O18O18O J = 0 j =7 j =5 E O2 OO Threshold j =4 ZPEO2 . 0 Energy (eV) j =3 j =5 PES j =2 j =1 j =0 O3 j =3 DZPE j =1 Lifetime (ps) The Lifetime Spectrum Babikov et al, J. Chem. Phys. 118, 6298 (2003).

16. 18O18O 16O18O 16O18O 16O16O 16O16O18O 16O18O18O j =7 j =5 j =5 j =6 j =4 j =5 j =3 Energy (eV) Energy (eV) j =3 j =5 j =2 j =4 j =1 j =1 j =0 j =3 j =3 DZPE j =2 DZPE j =1 j =0 j =1 Lifetime (ps) Lifetime (ps)

17. D ZPE ZPE1618 ZPE1818 Stable O3 16+1818 1618+18 PES 161818 Mechanism of DZPE Isotope Effect : 16O18O18O Babikov et al CPL 372, 686 (2003).

18. Rate: 0.92 Rate: 1.50 Mechanism of DZPE Isotope Effect : 16O18O18O Metastable O3* D ZPE ZPE1618 ZPE1818 Stable O3 16+1818 1618+18 PES Babikov et al CPL 372, 686 (2003). 161818

19. Quantum effect; • (classical trajectories • could not reproduce). • General effect. Mechanism of DZPE Isotope Effect : 16O16O18O “Background” Rate: 0.92 Rate: 1.45 Metastable O3* D ZPE ZPE1616 ZPE1618 Stable O3 16+1618 1616+18 PES Babikov et al CPL 372, 686 (2003). 161618

20. ; Channel-specific rate coefficients: ; Mechanism of DZPE Isotope Effect D. Babikov et al, J. Chem. Phys.119, 2577, 2003 A simple model: 161818 181616 3) correct order of magnitude. 2) is always in the right direction. 1) source of a very large isotope effect.

21. Introducing DZPE into Classical Trajectories (Schinke) • Although it is impossible to introduce rigorously the quantum ZPE into classical trajectory simulations, it appears feasible to mimic the DZPEeffect using a simple trick: • DZPE is introduced into the PES ad hoc; • The time classical trajectories for O + O2 • spend in the O3* region is determined; Smooth dependence, no resonances. • Stabilization probability is defined as: Adjustable parameters describe stabilization and are used to fit the experimental data: - stabilizing collision frequency ( ~ P ), - energy transfer efficiency (model), - symmetry effect (postulated). Schinke and Fleurat-Lessard, JCP 122, 94317 (2005).

22. 3D Coupled-Channel Treatment (Bowman) ● Same approach as Clary, but O3 wave functions are full dimensional (3D), as the PES. A basis set of 100 Legendre polynomials is used for the bending motion. ● Calculations were carried out only for ℓ = 0 (zero impact parameter), single value of collision energy, and three orientations of O3* (head, tail, perpendicular). ● The focus was on the DZPEeffect and the role of van der Waals states. • Results: • Confirmed importance of theDZPErange; • Proposed that the vdW states are important. Xie and Bowman, Chem. Phys. Lett. 412, 131 (2005).

23. J > 0 Calculations of the DZPE Effect (Grebenshchikov) • ● Centrifugally Sudden Approximation for rotation and the symmetric top rotor model are used. • ● PES is simplified by removing the vdW part, leaving only the covalent well. • ● Narrow resonances (G ≤ 1 cm-1) are determined for J ≤ 40, K ≤ 10. • ● First order perturbation theory is used to determine the branching ratios for two channels. • ● Stabilization is not treated, simple model is used (strong collision assumption). The bottom-line: Several different authors/methods show importance of the DZPE range. Although not yet modeled with full rigor, this effect appears to be fairly well understood at this point. Grebenshchikov and Shinke, JCP 131, 181103 (2009).

24. New Round of Theoretical and Experimental Research • Focus on: • Explaining very detailed experimental results on O + O2 scattering • studied in the molecular beam conditions; • Finding origin of the symmetry effect (h-factor).

25. Search for Origin of the Symmetry Effect (Schinke) ● Sudden Collision Approximation for energy transfer in the Ar + O3 collisions; IOSA and the Coupled-Channel formalism. ● Similar to Clary/Bowman, but only the bound states (below dissociation threshold) are taken into consideration. Q.: Inaccurate near threshold? ● The wave functions and the PES are full dimensional, but the PES is simplified by removing the vdWpart, leaving only the covalent well. Red – 686 Black – 866 Results: Expected symmetry effect was not observed... Agreement with classical trajectory simulations was surprisingly good (small effect in wrong direction). Ivanov and Schinke, Mol. Phys. 108, 259 (2010).

26. Quantum Symmetry Effect in a Model System (Pack) 16Ne + 18Ne  1618Ne2* (+ M)  1618Ne2 17Ne + 17Ne  1717Ne2* (+ M)  1717Ne2 State-to-state ( v, j ) rate coefficients in 1618: • This simple problem allows to carry out • very clean VRIOSA calculations: • masses are slightly modified in order • to have exactly the same reduced masses, • same energies and lifetimes of all states. • there is no DZPE here, any difference seen • is due to symmetry in the Ne2* + M collision. • number of states is small and easy to treat. 11% isotope effect ! • Results: • weak (strict) selection rules for state to state transitions in 1618 (1717). Quantum effect, classical trajectories would not reproduce. • - this opens additional pathways for the energy transfer and increases the recombination rates. Pack and Walker, JCP 121, 806 (2011).

27. Scaling Problems in the Quantum Dynamics Calculations • Scaling with J • Size of the Hamiltonian matrix scales linearly with J : 1296 x 1296 for J=0 • 41248 x 41248 for J=31 • Cost of calculations scales as ~J3, J2. Kendrick, JCP 114, 8796 (2001). • Scaling with Number of Atoms • - Number of vibrational degrees of freedom is 3N–6. • Most QM methods use Direct Product Basis sets (DPB) to express the wave function. As result: • exponential scaling problem. B. Poirier and co., J. Chem. Phys. 124, 144107 (2006).

28. Parallelization Problem ( Poirier ) It is found that when the quantum dynamics calculations are parallelized using standard math libraries the efficiency significantly drops after p ~10 or so (communication). ~ 2x107 B. Poirier has shown that using the sparsity pattern of the Hamiltonian matrix, a speedup close to linear and efficiency close to one can be achieved with large number of processors for large systems. • Block Jacoby diagonalization; • Domain decomposition; • Data distribution; • Load balancing. Chen and Poirier, J. Comp. Phys. 219 (2006) 185.

29. Adiabatic Bending Model (Babikov) Approach based on fast vs. slow degrees of freedom (Born-Oppenheimer like): • The bending angle is “relaxed” to convert 3D PES into a 2D PES,V(r1,r2). • The overall topology of the surface is preserved: • De, w, “reef”, vdW wells, channels. • Two channels allow studying the DZPE effect. • The effect of the bending enters through the bending energy correction and the partition function.  16O16O + 18O 16O + 16O18O Ivanov and Babikov, JCP 134, 174308 (2011).

30. Mixed Quantum-Classical Theory of Energy Transfer (Babikov) The internal vibrational motion is treated with QM using the TDSE (wave packet): resonances, DZPE and permutationsymmetry. 16O18O16O J=19, Ka=4, Kb=12 The collisional motion O3* + M and rotation of O3* are treated using classical trajectories: computational advantage. Energy is exchanged between translation, rotation, vibration. Total energy is conserved. Classical degrees of freedom allow intrinsic massive parallelization. Ivanov and Babikov, JCP 134, 144107 (2011). b(a0)

31. Can Quantum Isotope Effects Contribute to S-MIF ? Several gas-phase reactions involve S-atoms and may exhibit the DZPE and symmetry isotope effects: Gao and Marcus, JCP 127, 244316 (2007). CI Recombination by ET: S + S2 (+ M) → S3 Sm + Sn (+ M) → Sn+m SO + O (+ M) → SO2 SO2 + O (+ M) → SO3 SO + H (+ M) → HSO SO2+ OH (+ M) → HSO3 Several other: S + SH → S2+ H S + OH → SO+ H HS + O → H + SO S2 + O → S + SO S + O2 → SO + O Farquhar et al., J. Geophys. Res. 106, 32829 (2001). Pavlov and Kasting, Astrobiology 2, 27 (2002).

32. Relative Rates (Exp.) (Mauersberger and co., PCCP3, 4718, 2001)

33. Can Dynamics Methods be Applied to S-MIF ? Although not really a one-day job, the classical trajectory simulations can be carried out for variety of chemical reactions relatively easily. Size of the molecule - the number of degrees of freedom, is not really a problem. (Well, given the potential energy surface…) Example: Isotope effect in the O + NO (+ M) → NO2 recombination. Classical + ZPE method of Schinke was applied and predicted the isotope effect (larger than that in ozone): Ivanov and Schinke, JCP 126, 54304 (2007). Note: It is relatively straightforward to set up such calculations for SO2

34. Construction of the Potential Energy Surfaces Before the nuclear dynamics is studied, the electronic structure problem is solved for many nuclear configurations (independently). Dependence of electronic energy on nuclear configurations gives the potential energy surface. Motion of the nuclei on this surface (dynamics) is studied next. Two major methods for building a continuous surface from the descrete ab initio data points: - Spline (highly accurate, but practical only for small molecules); - Analytic fit (the only way to go in the case of larger polyatomics). Permutational invariance is important in the context of the isotope effects, which affects the choice of - Coordinates; - Functional form of the fitting function.

35. O O 117° • Spectroscopically accurate at low energies; • wrong behavior in the barrier region. • Dissociation energy: • DVQZ = 1.027 eV, • DEXP = 1.132 eV. O 2.4a0 “ Reef ” Along the Minimum Energy Path • Correction is smooth in 3D; • Only upper part of PES; • Corrects barrier and Dexp. • Van-der-Waals tail. Ground State PES of Ozone Ab initio electronic structure: MOLPRO: icMRCI+Q/cc-pVQZ, CASSCF(12,9). Spline on a 3D grid: 11 x 28 x 20 = 6160 points (Schinke and co., JCP116, 9749, 2002) (Babikov and co., JCP118, 6298, 2003)

36. r E = – 0.012 eV E = – 0.02 eV E = – 0.027 eV E = 0 E = – 0.03 eV E = – 0.1 eV O 1. E = – 0.2 eV O O O I. O O O O O O O O E = – 0.3 eV i. O O O E = – 0.4 eV 2. O O E = – 0.5 eV II. O E = – 0.6 eV ii. O E = – 0.7 eV 3. O O O III. E = – 0.8 eV O O E = – 0.9 eV E = – 1.0 eV Shadow Euler ( r,, ; , ,  ) Size Shape q, f 3D surface: 3D surface 1D Slice along the MEP: E = DZPE

37. Potential Energy Surface of S3 Isotopic shifts predicted for 32S3/34S3 mixture: There is no global potential energy surface available for S3 at this time. Exploratory work showed many similarities to O3: Analog of the Hartley band in ozone: • - Isoelectronic, structural similarities; • Two isomers, cyclic-S3 is at much • lower energies, 4.39 kcal/mol. • Calculations at MRCI+Q/CBS are • needed to reproduce spec. const. • - Covalent well is much deeper, 2.7 eV, • vibrational frequencies are smaller • (translates into density of states and • number of coupled channels). Their hypothesis: l~272 nm +M S + S2 →*S3 +hv *S3 → S2+ *S(1D) Francisco and co., JCP123, 54302 (2005). +R *S(1D) → OC*S, *SO2 Francisco and co., JCP125, 84314 (2006).

38. Potential Energy Surfaces of SO2 - accurate empirical PES for the ground X1A1 state. ~ ~ Ab initio PES for the C1B2 state: H. Guo, Chem. Phys. Lett., 329, 503 (2000).

39. Photo-absorption Spectra of SO2 Isotopomers (Gua) Calculated absorption spectra: Vibrational states of the excited PES (adiabatic calculations) bending mode progressions … … … … … • Intensities are quite similar • among the isotopomers. • Frequency shifts are regular. H. Guo, Chem. Phys. Lett. 439, 280 (2007).

40. Conclusions Very neat quantum mechanical effects lead to MIF in O3: - DZPE effect; - symmetry effect; - scattering resonances. Challenges for theory and experiment on ozone: - spectroscopically accurate PES near O + O2 threshold; - collisional stabilization of O3* and the symmetry effect; - dynamics and spectroscopy near threshold. Work done on the S-containing species: - Predictions of statistical theory for SO2; - PESs of SO2; preliminary work on S3; - Photo-absorption spectra of SO2 . In the near future: - Accurate PESs for S3 ( MRCI-CBS level ); - Classical trajectory studies for S3 and SO2( aka Schinke ); - Energy transfer in the S3 + Mcollisions ( mixed Q-C). Acknowledgments: NSF Atmospheric Chemistry Program ($$$)

41. j =180 deg. j =117 deg. j ~ 80 deg. j =0.