# Toward Understanding the Roaming Mechanism in H + MgH → Mg + HH Reaction

Frédéric Mauguière, Peter Collins, Stamatis Stamatiadis, Anyang Li, Gregory

Ezra, Stavros Farantos, Zeb Kramer, Barry Carpenter, Stephen Wiggins, Hua

Guo

► To cite this version:

Frédéric Mauguière, Peter Collins, Stamatis Stamatiadis, Anyang Li, Gregory Ezra, et al. Toward Understanding the Roaming Mechanism in H + MgH → Mg + HH Reaction. Journal of Physical Chemistry A, 2016, 120 (27), pp.5145-5154. 10.1021/acs.jpca.6b00682. http://dx.doi.org/10.1021/acs.jpca.6b00682

HAL Id: hal-05102469

https://univ-reims.hal.science/hal-05102469v1

Submitted on 7 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Toward Understanding the Roaming Mechanism in H + MgH → Mg + HH Reaction

Frédéric A. L. Mauguière,*,† Peter Collins,*,† Stamatis Stamatiadis,*,‡ Anyang Li,*,∇,§ Gregory S. Ezra,*,‖ Stavros C. Farantos,*,⊥ Zeb C. Kramer,*,‖ Barry K. Carpenter,*,# Stephen Wiggins,*,† and Hua Guo*,§

*†School of Mathematics, University of Bristol, Bristol BS8 1TW, United Kingdom

*‡Department of Materials Science and Technology, University of Crete, Iraklion 710 03, Greece

*§Department of Chemistry and Chemical Biology, University of New Mexico, Albuquerque, New Mexico 87131, United States

*‖Baker Laboratory, Department of Chemistry and Chemical Biology, Cornell University, Ithaca, New York 14853, United States

⊥Institute of Electronic Structure and Laser, Foundation for Research and Technology—Hellas, and Department of Chemistry, University of Crete, Iraklion 711 10, Greece

#School of Chemistry, Cardiff University, Cardiff CF10 3AT, United Kingdom

**ABSTRACT:** The roaming mechanism in the reaction H + MgH → Mg + HH is investigated by classical and quantum dynamics employing an accurate ab initio three-dimensional ground electronic state potential energy surface. The reaction dynamics are explored by running trajectories initialized on a four-dimensional dividing surface anchored on three-dimensional normally hyperbolic invariant manifold associated with a family of unstable orbiting periodic orbits in the entrance channel of the reaction (H + MgH). By locating periodic orbits localized in the HMgH well or involving H orbiting around the MgH diatom, and following their continuation with the total energy, regions in phase space where reactive or nonreactive trajectories may be trapped are found. In this way roaming reaction pathways are deduced in phase space.

Patterns similar to periodic orbits projected into configuration space are found for the quantum bound and resonance eigenstates. Roaming is attributed to the capture of the trajectories in the neighborhood of certain periodic orbits. The complex forming trajectories in the HMgH well can either return to the radical channel or “roam” to the MgHH minimum from where the molecule may react.

<figure><img src="image_2.png" /></figure>

## 1. INTRODUCTION

Reactions having product energy distributions that cannot be correlated in standard fashion with the most prominent features of the potential energy surface (PES) landscape are naturally of great interest. This is the case for the roaming mechanism¹⁻⁵ demonstrated almost a decade ago in studies of formaldehyde photodissociation, H₂CO → H₂ + CO. Quasiclassical trajectory calculations⁶ revealed a pathway in which formaldehyde first attempts unsuccessful, “frustrated”, dissociation to radical products, H' + HCO, with the hydrogen atom (H') making long-range excursions before turning back to react and form vibrationally excited hydrogen molecules.⁷

Following the pioneering work on formaldehyde, several other molecules have been studied and found to exhibit roaming mechanisms in reactions involving direct dissociation, intermediate isomerization,³ or even multiple PESs.⁸

We have recently analyzed the roaming phenomenon in ion-molecule reactions, a class of systems characterized by long-range charge/induced-dipole interactions.⁹,¹⁰ A useful model for such systems is the Chesnavich potential.¹¹ By investigating the dynamics of this system in its appropriate setting,¹²,¹³ phase space, we elucidated the roaming phenomenon in the presence of well-defined dividing surfaces (DS) acting as transitions states (minimal flux bottlenecks) and

examined related reaction pathways in phase space. We have also studied the isomerization of ketene;¹⁴ this system provides an example of the roaming phenomenon in types of reaction other than those exhibiting radical product channels.¹⁵,¹⁶

The phase space methodology applied to ion-molecule and ketene models has also been employed to elucidate roaming in formaldehyde.¹⁷ Using a 2D model potential extracted from the six-dimensional PES produced by Bowman and co-workers,¹⁸ we have answered the following question: Why do trajectories roam rather than dissociate through the radical channel? We demonstrated that phase space objects analogous to those found in the ion-molecule model define the roaming region and trap roaming trajectories in formaldehyde.

A recently studied molecule exhibiting roaming is magnesium hydride,¹⁹⁻²¹ a promising storage medium for hydrogen molecules and isoelectronic to formaldehyde.²²,²³ Li, Li, and Guo,²⁴ have carried out both quantum and classical trajectory

**Special Issue:** Piergiorgio Casavecchia and Antonio Lagana Feschtschrift

Received: January 21, 2016

Revised: February 26, 2016

Published: February 26, 2016
calculations on an improved ab initio ground electronic state PES¹⁹ to study the reactions

$$H' + MgH \rightarrow Mg + HH' \quad (\text{reactive: direct and roaming}) \quad (1)$$

$$H' + MgH \rightarrow H' + MgH \quad (\text{nonreactive: direct and roaming}) \quad (2)$$

$$H' + MgH \rightarrow H + MgH' \quad (\text{exchange}) \quad (3)$$

Several microscopic reaction channels were identified for the three elementary reactions: namely, direct and roaming abstraction (eq 1), a direct as well as a roaming nonreactive channel (eq 2), and an exchange reaction (eq 3). It was found that reaction 1 is dominated at low collision energies by the direct abstraction channel, whereas the exchange reaction, which involves a highly energetic intermediate complex, [MgH₂]*, plays almost no role at the collision energies studied. As in the formaldehyde dissociation, the two dominant channels (direct and roaming) produce similar highly excited vibrational distributions for the H₂ product. However, it should be noted that the energetic complex is prepared here by collision rather than photoabsorption.

In the case of formaldehyde, attempts had previously been made to identify roaming trajectories¹,⁶ as those trajectories in which hydrogen undergoes large-amplitude motions in the radical dissociation channel and which pass through a large region of configuration space near a so-called “roaming transition state”.²¹,²⁵ On the basis of their classical trajectory studies, Li et al.²⁴ concluded that it is difficult to define such trajectories unambiguously solely on the basis of their passage near the geometry of the putative “roaming transition state” and so defined roaming trajectories as those reactive trajectories in which the H–H' distance extends beyond 8.5 a₀ after the first turning point in the R coordinate. Such a criterion is similar to that used by Christoffel and Bowman²⁶ in identifying various types of trajectories in the H' + HCO reaction. Such criteria for roaming are not unique and, being based solely in configuration space, do not necessarily have a sound basis in dynamics.

Li et al.²⁴ have argued that for a rotationless MgH₂, “roaming” is quantum mechanically associated with a large-amplitude vibrational progression that emerges below the radical reaction threshold and continues into the continuum leading to “roaming resonances”.

In the present article we further investigate the mechanism of roaming dynamics for reaction 1 by seeking the invariant phase space objects that cause roaming. The remainder of the article is organized as follows: section 2 briefly describes the PES and the numerical methods employed. In section 3 both quantum and classical aspects of the dynamics of reactions 1–3 are examined, while in section 4 we discuss the implications of our work in understanding the roaming effect and its generalization in reaction dynamics. Section 5 concludes.

## 2. METHODS

### 2.1. Ab Initio Potential Energy Surface.
An ab initio global PES fitted by 3D spline interpolation²⁷ has been constructed for the ground electronic state (¹A') of MgH₂. Initially, Li et al.¹⁹ carried out calculations using the internally contracted multireference configuration interaction method with Davidson correction (icMRCl+Q) and with the cc-pVnZ (n = 3, 4, and 5) basis sets extrapolated to the complete basis

set limit. This electronic state correlates with both the Mg(¹S₀) + H₂ and MgH(X²Σ⁺) + H' asymptotes, thus facilitating the study of reactions 1, 2, and 3 up to bond length distances of 9 a₀. Recent calculations²⁴ have extended the validity of this PES to longer distances (up to 15 a₀) and have provided a total of 5406 ab initio points.

The configuration space is described by the Jacobi coordinate system most appropriate for the radical channel: R, the distance between the center of mass of MgH and hydrogen atom H', r, the bond length of the diatom MgH, and θ, the angle between the two distances; θ = π is the H'MgH and θ = 0 the MgHH'. The calculations have been performed in atomic units (m_e = e = ħ = 1), and angles are in radians and degrees. Energies are mainly expressed in hartrees, but other units such as kcal/mol and wavenumber are also used.

### 2.2. Periodic Orbit Calculations.
A powerful method to explore the phase space structure of a nonlinear dynamical system for extended ranges of energy (or other system parameters) involves the study of periodic orbits and their continuation as energy or other parameters vary. Periodic orbits have been located with the program POMULT,²⁸ and first and second derivatives of the potential needed in the multiple shooting method have been computed by the program AUTODERIV²⁹ in Jacobi coordinates. The Hamiltonian employed is that of a triatomic molecule in Jacobi coordinates with zero total angular momentum,

$$H(r,R,\theta,p_r,p_R,p_\theta) = \frac{p_r^2}{2\mu_r} + \frac{p_R^2}{2\mu_R} + \frac{p_\theta^2}{2} \left( \frac{1}{\mu_R R^2} + \frac{1}{\mu_r r^2} \right) + V(r,R,\theta) \quad (4)$$

where p_r, p_R, and p_θ are momenta canonically conjugate to coordinates r, R, and θ, respectively, and μ_R and μ_r are associated reduced masses.

Principal families of PO emanate from equilibrium points, both stable and unstable, and their existence has been proved by Weinstein³⁰ and Moser.³¹ At critical values of the energy, bifurcations take place and new families are born. Continuation/bifurcation diagrams are obtained by plotting PO frequency versus energy.

### 2.3. Quantum Dynamical Calculations.
The Schrödinger equation was solved by discretizing the wave function in a mixed representation consisting of a direct product discrete variable representation (DVR)³²,³³ for the two radial degrees of freedom and a finite basis representation (FBR) for the angular degree of freedom. The propagation was initiated with a wave packet located at R = 6.0 a₀.

The Chebyshev autocorrelation function was computed, and the energy dependent wave functions for both bound and resonance states were reconstructed using the following equation³⁴,³⁵

$$\psi(E) = \sum_{k=0}^{K} (2 - \delta_{k0}) \cos[k \arccos(E)] \psi_k \quad (5)$$

based on the order/angle representation of the Chebyshev operator. E is the total energy and δk0 the Kronecker delta function. Here, ψk is the kth order Chebyshev wave packet, defined as follows:

$$\psi_k = 2\hat{H}\psi_{k-1} - \psi_{k-2} \quad (6)$$

with ψ₁ = 2Hψ₀ and H as the normalized Hamiltonian.
The quantum spectrum is obtained directly from a Fourier transform of a truncated Chebyshev correlation function.

**2.4. Sampling the Orbiting Transition State.** The phase space approach to transition state theory (TST)³⁶ of chemical reactions involving polyatomic molecules requires the objects, namely, *normally hyperbolic invariant manifolds* (NHIM),³⁷ for the construction of dividing surfaces that locally minimize the reaction flux.¹²

In our previous studies,⁹,¹⁰,¹⁷ which involve barrierless dissociation (association) channels on the PES, we have shown the existence of a loose transition state in the exit (entrance) channel associated with the centrifugal barrier.³⁸ For two degrees of freedom (DoF) systems, the NHIM is an unstable orbiting periodic orbit that supports a dividing surface, i.e., the outer ("loose") transition state, that separates the dissociation products from the interacting complex. The dividing surface is 2D, and every dissociating (associating) trajectory must cross it. In these studies we have described in detail how to sample this surface at constant energy in order to assign initial conditions to trajectories.

The theory has recently been extended³⁹ to the case of a diatomic molecule weakly coupled to an atom. The NHIM in this case is made up of a one-parameter family of 2D tori, where the parameter of this family determines the distribution of energy between the two DoF subsystem ($R, \theta$), which describes the motion of the atom at large $R$ distances ($R/r \gg 1$) and the elliptic diatomic oscillator ($r$). In other words, it has been shown that each of the 2D tori composing the NHIM is the product of two periodic orbits, that of the unstable orbiting PO at the entrance channel of the atom–diatom system with the PO of the diatomic oscillator. Anchored to this NHIM are dividing surfaces that are transversely crossed by the reactive trajectories. Details of how to sample such a DS are given in ref 39, and here we apply the method in the $H' + MgH$ collisions.

The selected trajectories are integrated in time, and they are terminated when some reaction criteria are satisfied; in the case of $MgH_2$ we examine the interatomic distances to discern among various reactive events. The total integration times are analyzed by producing gap time distributions.¹⁰,⁴⁰ The definitions used to compute gap times are discussed in Appendix A.

## 3. RESULTS

**3.1. Potential Energy Surface and Equilibria.** The potential is examined as a function of the three Jacobi coordinates ($r, R, \theta$). Contour plots of the PES in the ($R, \theta$) plane and for $r = 3.28 a_0$ are shown in Figure 1a. Five stationary points of the potential can be identified in this figure. The absolute minimum, H'MgH, is located at $\theta = \pi$, $R = 3.36 a_0$, and $r = 3.23 a_0$, whereas a metastable minimum at $\theta = 0$, MgHH', can also be seen. The saddle between these two minima at $R = 3.37 a_0$, $r = 3.45 a_0$ and $\theta = 0.962$ rad is of index-1⁴¹ and lies 0.1121 hartree above the absolute minimum (H'MgH).

Along the conventional minimum energy path (MEP) starting at angle $\theta = \pi$, the hydrogen atom $H'$ approaches MgH with $\theta \approx \pi$ and then passes over to the MgHH' minimum at $\theta = 0$ by overcoming the saddle (tS) at $\theta = 0.962$ rad to give products Mg + $H_2$. This saddle is of $C_{2v}$ symmetry and is associated with the conventional transition state for reaction 1. Another abstraction pathway on the PES has $H'$ approaching MgH along $\theta \approx 0$ (see also Figure 1 in ref 24).

<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Potential energy contours and periodic orbits originated from center–center–saddle bifurcations (in brief center–saddle), and projected in the ($R, \theta$) plane are shown. Except those PO that appear above the saddle rS, marked by the red point at about $R = 8 a_0$, such as the orbiting transition state (OTS) (blue line) located at $R = 13.5 a_0$ and the orbiting periodic orbit (red) that passes close to the rS, all the other PO span an energy range from below to above the threshold to radical products, which lies 73.26 kcal/mol (0.1168 hartree) above the H'MgH absolute minimum. (b) Potential energy contours and periodic orbits of Lyapunov families emanated from the two saddles index-1, tS and rS, projected in the ($R, r$) plane. For both saddles the unstable direction is predominantly along $\theta$. Contours are in kcal/mol (right axis).</figcaption></figure>

In Figure 1b we plot contours of the PES in the ($R, r$) plane and for $\theta = 0.962$ rad. In this projection the tS saddle (index-1) appears as a minimum with the unstable direction predominantly in bend. In the same figure we find minima in the entrance and exit channels, which are in fact two symmetry-related index-1 saddles with the unstable direction also predominantly in the bend. The saddle is located at $\theta = 1.013$ rad, $R = 8.24 a_0$, $r = 3.29 a_0$ and with energy 0.1166 hartree above H'MgH. In Figure 1a it is denoted as rS and is shown with a red dot. It is just 0.05 kcal/mol below the dissociation threshold to radical products (0.1168 hartree above the absolute minimum).

Between the two index-1 saddles, tS and rS, there is an index-2 saddle at 0.1194 hartree above the H'MgH minimum (1.7 kcal/mol above the threshold to radical channel) with geometry $R = 4.25 a_0$, $r = 3.35 a_0$, and $\theta = 1.001$ rad.

Guo and co-workers²⁴ as well as Harding et al.²¹ have attributed roaming to the presence of the rS saddle. As we shall see, in our analysis the roaming region is defined by other phase space structures in the entrance channel.

**3.2. Quantum Mechanical Spectrum and Resonances.**
The Schrödinger equation has been solved by propagating an initial Gaussian wave packet centered at $R = 6.0 a_0$ inside the
H'MgH well. Both bound and resonance states near the radical channel threshold have been extracted and some of their wave functions plotted. Figure 2 shows the quantum spectrum, where

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Quantum mechanical spectrum in an energy interval that spans bound and resonance eigenstates of H'MgH. The red arrows assign those eigenstates whose wave functions are plotted in Figure 3.</figcaption></figure>

the intensity is plotted as a function of the energy relative to the asymptotic potential. A rather regular pattern is revealed for energies below the reaction threshold with energy intervals between successive lines decreasing. The spectrum appears less regular close to the dissociation limit, but above this limit we can clearly see three broad peaks separated by approximately 0.3 kcal/mol at 0.21, 0.53, and 0.83 kcal/mol collision energy.

The red arrows mark those bound and resonance states whose wave functions are projected in the (R, θ) plane in Figure 3. In these plots we depict contours of the amplitude of the eigenfunction.

A smooth transition of the bound vibrational states into the continuum is evident. In particular, we note the presence of excitations along both *R* and *θ* coordinates in all wave functions except one, the state with energy -0.02 kcal/mol below the dissociation threshold and above the saddle rS, which has a very

simple nodal pattern. The nodal patterns for the other bound wave functions appear regular, and we can clearly distinguish two parts of the wave function: one localized in the H'MgH minimum, with the other consisting of a branch pointing toward the rS saddle. However, states like -0.09 and +0.05 forming an arc around $θ = π$ are observed.

Similarly, the resonance wave functions above dissociation have two parts; the first remains localized inside the H'MgH well in the same way as for the bound states (but less regular), while the other branch is now extended in configuration space beyond the saddle rS. The state at energy of +0.53 kcal/mol has amplitude spread in the [0, π] domain of the angle θ. Notice that the wave functions are plotted in the angle interval 0 ≤ θ ≤ π, exploiting the symmetry of the system, whereas the PO in Figure 1 is in the range of 2π.

### 3.3. Periodic Orbit Analysis and Phase Space Structures.
A systematic exploration of the periodic orbit continuation/bifurcation diagram requires first the location of the principal families of PO emerging from the equilibria of the dynamical system and then their continuation with the energy.<sup>13</sup> Here, we study PO at energies of interest below and above the dissociation channel to radical products of MgH<sub>2</sub> relevant to roaming mechanism that we want to investigate.

Simple (in the projection of the (R, θ) plane) families of PO for MgH<sub>2</sub> are shown in Figure 1a. The bands represent several PO of each family spanning a range of energies from below to above the dissociation threshold of MgH<sub>2</sub>. The magenta band consists of periodic orbits inside the H'MgH well connecting the two symmetry-related saddles tS. The short red lines at the angle of approximately 1 rad are Lyapunov unstable periodic orbits which emanate from the tS. Two such families exist and are more clearly shown in Figure 1b as red lines representing the symmetric and antisymmetric stretch vibrational modes expected for the C<sub>2v</sub> symmetry saddle tS.

Similarly, the red spot at about R = 8 a<sub>0</sub> corresponds to the saddle rS. In Figure 1b the Lyapunov type unstable periodic orbits appear as orthogonal lines, indicating the decoupling of the *R* and *r* modes at this region of configuration space. The blue line in Figure 1a depicts PO that appear just above the rS saddle and below the dissociation threshold, whereas the red

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Plots of wave functions for bound with negative energies in kcal/mol and resonance with positive energies eigenstates. The wave functions correspond to the eigenstates marked by the red arrows in the quantum spectrum (Figure 2). Distances are in au and angles in degrees.</figcaption></figure>
line depicts rotating periodic orbits, i.e. the angle $\theta$ covers the range of $[0, 2\pi]$.

At $R = 13.5$ $a_0$ the blue line is a rotating periodic orbit that defines the orbiting transition state (OTS)³⁸ and corresponds to a relative equilibrium with energy 0.5 kcal/mol above the dissociation threshold. We must emphasize that a family of such rotating periodic orbits exists and a different PO must be defined according to the partition of the energy between the 2D ($R, \theta$) subsystem and the decoupled MgH oscillator $r$.³⁹ In this work we examine the dynamics of the system by initializing the trajectories with the MgH diatom at its equilibrium bond length.

The thick green band in Figure 1a consists of periodic orbits emanating from a center-saddle bifurcation below the dissociation energy and continuing above. Figure 4 depicts

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Continuation/bifurcation diagram projected in the energy-frequency plane for the green family of periodic orbits shown in Figure 1a.</figcaption></figure>

the continuation/bifurcation diagram of this family of periodic orbits (green) in the energy-frequency projection plane. Notice the abrupt change of frequency, about $20 \text{ cm}^{-1}$, in a short energy interval, indicating high anharmonicity. This is typical for the periodic orbits we have found in the energy region close and above the saddle $rS$. Finally, the violet band consists of PO located in the MgHH' well.

Figure 5 depicts four different type of periodic orbits born inside the H'MgH well and extending above the dissociation threshold. As in Figure 1a, the bands comprise several PO. The asymmetric PO (Figure 5c, Figure 5d) have a lobe extended toward the $rS$. The symmetric PO in Figure 5a are localized in the well, whereas PO in Figure 5b have a branch that extends up to $R = 10$ $a_0$ and along the $\theta = \pi$ reaction pathway. All of these type of periodic orbits are associated with high order resonances among the three degrees of freedom; their shapes are reflected in the patterns of the wave functions shown in Figure 3.

Simple semiclassical arguments correlate the PO frequencies with the energy gap between successive quantum energy levels ($\Delta E = h/T_{PO}$, where $T_{PO}$ is the period), and indeed this relation holds for the green and red periodic orbits shown in Figure 1a. The energy gap of 0.3 kcal/mol seen in the quantum spectrum above dissociation corresponds to approximately 0.3 ps, and this is in the range of periods found for some periodic orbits shown in Figures 1 and 5.

**3.4. Classical Trajectory Simulations.** In the previous subsections we have studied stationary states in classical and quantum mechanics of magnesium hydride by plotting periodic

orbits and eigenfunctions, respectively, and have found a qualitative agreement. While periodic orbits describe the local behavior of nearby trajectories, often in experiments broad regions of phase space are accessible, such as in studied bimolecular collisions. In simulations of such reactions, initial conditions for trajectories, which correspond to initial states of the two colliding species, must be sampled. This requires sampling of the proper dynamically defined dividing surfaces,¹² which are found by locating NHIM.³⁷

To perform a classical trajectory simulation for H' colliding with MgH, we sample initial conditions on the OTS at constant energy, as is described in detail by Mauguiere et al.³⁹ Propagating "incoming" trajectories, we find that three possible reactive outcomes, eqs 1, 2, and 3, may occur. We use colors to distinguish these outcomes: red for reactive trajectories (eq 1), blue for nonreactive trajectories (eq 2), and green for exchange reactions (eq 3). Trajectories that have not reacted within a prespecified maximum integration time (1 ps) are terminated (black squares).

In Figure 6a, we plot the initial conditions of the sampled trajectories on the OTS in the ($\theta_0, p_{\theta_0}$) plane. This plot reveals alternating bands of different trajectory types on the DS as has been found before (see, for example, refs 9, 10, and 40). The arrangement of bands can be very complicated (fractal).⁴⁰ In Figure 6b we plot gap times¹⁰,⁴⁰ versus $\theta_0$ for initial conditions on the OTS at fixed $p_{\theta_0} = 0$. The plot demonstrates the fractal nature at the border of bands associated with different trajectory types and indicates that gap times diverge at the boundary between bands associated with two different trajectory types. Infinite gap times correspond to trajectory initial conditions that are on the stable invariant manifolds of stationary objects such as hyperbolic tori.

In Figures 7–9 trajectories from the three characteristic classes of Figure 6 are plotted in the ($R, \theta$) plane. Figure 7 and Figure 9 cover the range of $\theta_0$ [0, 1] rad and (1.6, 2.14] rad, respectively. Trajectories in the interval (1, 1.6] rad are plotted in Figure 8. In the Figures 7 and 8 all types of trajectories are recorded: direct and roaming reactions (red), nonreactive ones (blue), and a few green that correspond to the hydrogen exchange reaction. In these calculations direct reactive trajectories are found to be those with the hydrogen atom approaching MgH along the $\theta \approx 0$ rad pathway. Roaming trajectories are those for which the H atom approaches the H'MgH well and does not pass over the tS saddle but, instead, returns to the radical channel and rotates around MgH before abstracting the other hydrogen atom to form Mg + H₂. Only a few trajectories were found to pass over the tS to form an energetic complex. In the angle domain (1.6, 2.14] rad (Figure 9) we find that trajectories follow a specific roaming reaction path to pass from the H'MgH well to the MgHH' well resembling that of the green and red periodic orbits in Figure 1a. Initial conditions from the range of (2.14, $\pi$] lead exclusively to nonreactive trajectories.

## 4. DISCUSSION

By systematically sampling³⁹ the OTS in the entrance of the radical product channel, we have recorded reactive (eq 1), nonreactive (eq 2), and exchange (eq 3) trajectories. However, as shown by their projections in the ($R, \theta$) plane in Figures 7–9, qualitatively different types of trajectories may lead to the same reactive or nonreactive event. This implies that different reaction paths are followed. Indeed, in the direct abstraction reaction, hydrogen approaches the other hydrogen along the
<figure><img src="image_7.png" /><figcaption>Figure 5. Periodic orbits originated from center-saddle bifurcations representing high order resonances. Except the upper-left ones, located for energies below the threshold to radical products, all the others span an energy range from below to above the threshold to radical products.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. (a) Distribution of the different type of trajectories on the OTS at collision energy of *E* = 0.5 kcal/mol. (b) Gap times for trajectories sampled from the OTS dividing surface and at *p*<sub>θ₀</sub> = 0. Red signs denote reactive trajectories, blue nonreactive, green the exchange reactions, and black the terminated trajectories after integrating for 1 ps.</figcaption></figure>

range of the small angles, $\theta \in [0, 1]$ rad, as can be seen in Figure 7, whereas roaming reactive trajectories are most likely in the range of angles (1.0, 2.14] rad (Figures 8, 9). In these trajectories, the hydrogen atom preferentially bounces back to the radical channel, rotates partially or completely around MgH, and then reacts. This is the “roaming mechanism” as has been identified in formaldehyde.<sup>3,17</sup>

On the other hand, we have found a few trajectories (green) that can surmount the barrier associated with the tS saddle (see Figure 1a). This barrier lies 2.9 kcal/mol below the dissociation threshold to radical products, and thus, the molecule has enough energy to follow the conventional reaction path, which nevertheless is avoided. This result is in accord with what Guo and co-workers<sup>24</sup> have found in their six-dimensional classical trajectories studies. The low rate of the exchange reaction
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Trajectories sampled from the OTS in the range of $\theta_0$ of [0–1] rad are plotted in the plane $(R, \theta)$. Red areas are for reactive trajectories, blue for nonreactive, and green for exchange reaction.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Trajectories sampled from the OTS in the range of angles (1–1.6] rad. Red areas are for reactive trajectories, blue for nonreactive, and green for exchange reaction.</figcaption></figure>

reveals inefficiency in energy transfer among the *r* vibrational mode and the other DoF and, thus, nonstatistical behavior of MgH₂.

We observe that configuration space plots of the trapped trajectories closely match those of the resonant PO shown in Figures 1 and 5. Representative trajectories are depicted in Figure 10. Panels a, e, f depict terminated trajectories after reaching the maximum integration time. Panels b, c, and d show nonreactive trajectories, and panels g, h, i, j depict reactive trajectories. Finally, panels k and l show an exchange reactive trajectory.

The PO are located in distinct regions of phase space, and they can trap trajectories in their vicinity. We demonstrate this

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Trajectories sampled from the OTS in the range of angles (1.6–2.14] rad. Red areas are for reactive trajectories.</figcaption></figure>

trapping by constructing Poincaré surfaces of section (PSS). In Figure 11a and Figure 11b we plot projections of PSS for bound trajectories sampled close to the saddle rS and to the green periodic orbits shown in Figure 1a. The projected points are found to lie in bands, indicating nonergodic behavior of the trajectories. Thus, the void in the center of Figure 11a shows that energy does not flow to the *r* mode, which explains the low rate of overcoming the tS saddle.

In support of this conclusion in Figure 12 we show the projection in the (*R*, *r*) plane of two periodic orbits of blue and green type depicted in Figure 1a. Although the MgH bond length (*r*) varies at small *R* values, not enough energy is transferred into the *r* mode to overcome the tS barrier.

In previous work⁹,¹⁰,¹⁴,¹⁷ we have described roaming as a dynamical effect in which the molecule explores alternative reaction pathways in phase space through trapping of trajectories in a certain domain of phase space marked by resonant periodic orbits. It turns out that a similar scenario exists for magnesium hydride with the orbiting transition state in the radical channel being the entrance “portal”. Nevertheless, for magnesium hydride we find that several different resonant conditions are involved in the mechanism of energy transfer from the radial (R) to the angular DoF. The important outcome of the present work is that by systematically sampling trajectory initial conditions on the OTS, we are able to unveil the structure of phase space.

The correspondence between periodic orbits and quantum mechanical eigenfunctions has been extensively demonstrated.⁴² The same is found for MgH₂. The computed eigenfunctions²⁰,²⁴ show localization with a gradual extension of the amplitude into the radical channel as energy increases and even above the dissociation threshold, giving rise to “roaming resonances”.²⁴

## 5. CONCLUSIONS

The mechanisms of the H⁺ + MgH ↔ [MgH₂]* ↔ Mg + HH' reaction have been investigated by examining classical trajectories. A 4D phase space dividing surface that separates the radical reactants from the complex [MgH₂]* can be rigorously defined and partially sampled in order to run trajectories that result in all possible reactive or nonreactive
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Representative trajectories sampled on the OTS. For the 3D representations the vertical axis depicts the angle $\theta$. Red areas are for reactive trajectories, blue for nonreactive, green for an exchange reaction, and black for terminated trajectories after integrating up to 1 ps.</figcaption></figure>

events.³⁹ The classification of the trajectories into different reactive events demonstrates a nonstatistical dynamical behavior of the highly energetic MgH₂. The same conclusion is drawn by locating families of periodic orbits, which serve to

label the regions in phase space where the three internal degrees of freedom are in resonance. The quantum mechanical spectrum and eigenfunctions at energies close to the threshold to radical products also demonstrate regular behavior and
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Projected Poincaré surfaces of section: (a) at energy *E* = 0.116 67 hartree and plane of section *θ* = *π* and (b) at energy *E* = 0.116 78 hartree and plane of section *r* = 3.28 *a*₀.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Periodic orbits above the saddle *r*S are projected in the (*R*, *r*) plane. Blue areas are PO at energy 0.116 669 hartree and green at 0.118 421 hartree. Distances are in au.</figcaption></figure>

localization of wave functions in correspondence with the phase
space objects such as periodic orbits.

In agreement with previous studies⁹,¹⁰,¹⁴,¹⁷ the roaming
mechanism is attributed to the trapping of trajectories in
specific regions of phase space marked by PO, where energy
can be transferred from *R* to *θ*, thus enabling trajectories to
“roam”. However, for MgH₂ there is no a single roaming
pathway but rather several different pathways by which
trajectories explore phase space. It has been demonstrated
that the sampled dividing surface (OTS) acts as a “lens” to
unveil phase space structures and possible reaction pathways.

In the present article we have explored the phase space of
MgH₂ by locating a number of PO. However, significant
questions remain to be explored. For example, what is the
mechanism of isomerization H'MgH ↔ MgHH'? What are the
roles of the red rotating PO for isomerizing by roaming and of
the magenta PO of Figure 1a for following the conventional
reaction pathway? These will be the subject of future work.

■ APPENDIX A: GAP TIME DEFINITIONS

In this appendix we give the definitions needed to compute gap
time distributions<sup>10,40</sup> for trajectories initiated on the OTS.

Trajectories are stopped when they meet the following
criteria corresponding to the three reactions.

Reaction 1: *R* > *R*<sub>thres</sub> and *r* > 8 *a*<sub>0</sub>.

Reaction 2: R > R<sub>thres</sub> and r < 8 a<sub>0</sub>.

Reaction 3: R < 8 a₀ and r > R<sub>thres</sub>.

The threshold value of the coordinate *R* is taken to be *R*<sub>thres</sub> =
14.5 *a*<sub>0</sub>.

■ AUTHOR INFORMATION

**Corresponding Authors**

*F.A.L.: e-mail, frederic.mauguiere@bristol.ac.uk.

*P.C.: e-mail, peter.collins@bristol.ac.uk.

*S.S.: e-mail, stamatis@materials.uoc.gr.

*A.L.: e-mail, liay@nwu.edu.cn.

*G.S.E.: e-mail, gse1@cornell.edu.

*S.C.F.: e-mail, farantos@iesl.forth.gr.

*Z.C.K.: e-mail, zck3@cornell.edu.

*B.K.C.: e-mail, CarpenterB1@cardi.ac.uk.

*S.W.: e-mail, stephen.wiggins@mac.com.

*H.G.: e-mail, hguo@unm.edu.

**Present Address**

A.L.: College of Chemistry and Materials Science, Northwest
University, Xi'an, 710069, China.

**Notes**

The authors declare no competing financial interest.
This is a theory paper; no raw data exists, and all evidence
underpinning this research is contained within the article.

■ ACKNOWLEDGMENTS

This work is supported by the National Science Foundation
under Grant CHE-1223754 (to G.S.E.). F.A.L.M., P.C., and
S.W. acknowledge the support of the Office of Naval Research
(Grant N00014-01-1-0769), the Leverhulme Trust. B.K.C.,
F.A.L.M, P.C., and S.W. acknowledge the support of the
Engineering and Physical Sciences Research Council (Grant
EP/K000489/1). H.G. thanks the Department of Energy for
support (Grant DF-FG02-05ER15694). S.C.F. is grateful to the
Institute for Advanced Studies, University of Bristol, for the
BENJAMIN MEAKER visiting professorship during September
15 to December 15, 2014.

■ REFERENCES

(1) Suits, A. G. Roaming Atoms and Radicals: A New Mechanism in Molecular Dissociation. *Acc. Chem. Res.* **2008**, *41*, 873–881.

(2) Herath, N.; Suits, A. G. Roaming Radical Reactions. *J. Phys. Chem. Lett.* **2011**, *2*, 642–647.

(3) Bowman, J. M. Roaming. *Mol. Phys.* **2014**, *112*, 2516–2528.

(4) Hase, W. L. Simulations of Gas-Phase Chemical Reactions. Applications to SN₂ Nucleophilic Substitution. *Science* **1994**, *266*, 998–1002.

(5) Hase, W. L. Some Recent Advances and Remaining Questions Regarding Unimolecular Rate Theory. *Acc. Chem. Res.* **1998**, *31*, 659–665.

(6) Townsend, D.; Lahankar, S. A.; Lee, S. K.; Chambreau, S. D.; Suits, A. G.; Zhang, Z.; Rheinecker, J.; Harding, L. B.; Bowman, J. M.
The Roaming Atom: Straying from the Reaction Path in Formaldehyde Decomposition. *Science* **2004**, *306*, 1158-1161.

(7) Lahankar, S. A.; Chambreau, S. D.; Zhang, X.; Bowman, J. M.; Suits, A. G. Energy Dependence of the Roaming Atom Pathway in Formaldehyde Decomposition. *J. Chem. Phys.* **2007**, *126*, 044314.

(8) Grubb, M. P.; Warter, M. L.; Xiao, H.; Maeda, S.; Morokuma, K.; North, S. W. No Straight Path: Roaming in Both Ground- and Excited-State Photolytic Channels of NO₃ → NO + O₂. *Science* **2012**, *335*, 1075-1078.

(9) Mauguière, F. A. L.; Collins, P.; Ezra, G. S.; Farantos, S. C.; Wiggins, S. Multiple Transition States and Roaming in Ion-Molecule Reactions: A Phase Space Perspective. *Chem. Phys. Lett.* **2014**, *592*, 282-287.

(10) Mauguière, F. A. L.; Collins, P.; Ezra, G. S.; Farantos, S. C.; Wiggins, S. Roaming Dynamics in Ion-Molecule Reactions: Phase Space Reaction Pathways and Geometrical Interpretation. *J. Chem. Phys.* **2014**, *140*, 134112.

(11) Chesnavich, W. J. Multiple Transition States in Unimolecular Reactions. *J. Chem. Phys.* **1986**, *84*, 2615-2619.

(12) Waalkens, H.; Schubert, R.; Wiggins, S. Wigner's Dynamical Transition State Theory in Phase Space: Classical and Quantum. *Nonlinearity* **2008**, *21*, R1-R118.

(13) Farantos, S. C. *Nonlinear Hamiltonian Mechanics Applied to Molecular Dynamics: Theory and Computational Methods for Understanding Molecular Spectroscopy and Chemical Reactions*; Springer, 2014.

(14) Mauguière, F. A. L.; Collins, P.; Ezra, G. S.; Farantos, S. C.; Wiggins, S. Roaming Dynamics in Ketene Isomerization. *Theor. Chem. Acc.* **2014**, *133*, 1507-1513.

(15) Ulusoy, I. S.; Stanton, J. F.; Hernandez, R. Effects of Roaming Trajectories on the Transition State Theory Rates of a Reduced-Dimensional Model of Ketene Isomerization. *J. Phys. Chem. A* **2013**, *117*, 7553-7560.

(16) Ulusoy, I. S.; Stanton, J. F.; Hernandez, R. Effects of Roaming Trajectories on the Transition State Theory Rates of a Reduced-Dimensional Model of Ketene Isomerization. *J. Phys. Chem. A* **2013**, *117*, 10567.

(17) Mauguière, F. A. L.; Collins, P.; Kramer, Z. C.; Carpenter, B. K.; Ezra, G. S.; Farantos, S. C.; Wiggins, S. Phase Space Structures Explain Hydrogen Atom Roaming in Formaldehyde Decomposition. *J. Phys. Chem. Lett.* **2015**, *6*, 4123-4128.

(18) Zhang, X.; Zou, S.; Harding, L. B.; Bowman, J. M. A Global Ab Initio Potential Energy Surface for Formaldehyde. *J. Phys. Chem. A* **2004**, *108*, 8980-8986.

(19) Li, H.; Xie, D.; Guo, H. An Ab Initio Potential Energy Surface and Vibrational States of MgH₂ (¹¹A'). *J. Chem. Phys.* **2004**, *121*, 4156-4163.

(20) Takayanagi, T.; Tanaka, T. Roaming Dynamics in the MgH + H → Mg + H₂ Reaction: Quantum Dynamics Calculations. *Chem. Phys. Lett.* **2011**, *504*, 130-135.

(21) Harding, L. B.; Klippenstein, S. J.; Jasper, A. W. Separability of Tight and Roaming Pathways to Molecular Decomposition. *J. Phys. Chem. A* **2012**, *116*, 6967-6982.

(22) Malinova, T.; Guo, Z. X. Artificial Neural Network Modelling of Hydrogen Storage Properties of Mg-Based Alloys. *Mater. Sci. Eng., A* **2004**, *365*, 219-227.

(23) Song, Y.; Guo, Z. X.; Yang, R. Influence of Titanium on the Hydrogen Storage Characteristics of Magnesium Hydride: A First Principles Investigation. *Mater. Sci. Eng., A* **2004**, *365*, 73-79.

(24) Li, A.; Li, J.; Guo, H. Quantum Manifestation of Roaming in H + MgH → Mg + H₂: The Birth of Roaming Resonances. *J. Phys. Chem. A* **2013**, *117*, 5052-5060.

(25) Klippenstein, S. J.; Georgievskii, Y.; Harding, L. B. Statistical Theory for the Kinetics and Dynamics of Roaming Reactions. *J. Phys. Chem. A* **2011**, *115*, 14370-14381.

(26) Christoffel, K. M.; Bowman, J. M. Three Reaction Pathways in the H + HCO → H₂ + CO Reaction. *J. Phys. Chem. A* **2009**, *113*, 4138-4144.

(27) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. *Numerical Recipes in FORTRAN: The Art of Scientific Computing*; Cambridge University Press: Cambridge, U.K., 1992.

(28) Farantos, S. C. POMULT: A Program for Computing Periodic Orbits in Hamiltonian Systems Based on Multiple Shooting Algorithms. *Comput. Phys. Commun.* **1998**, *108*, 240-258.

(29) Stamatiadis, S.; Farantos, S. C. AUTO DERIV: Tool for Automatic Differentiation of a Fortran Code (New Version). *Comput. Phys. Commun.* **2010**, *181*, 1818-1819.

(30) Weinstein, A. Normal Modes for Nonlinear Hamiltonian Systems. *Inv. Math.* **1973**, *20*, 47-57.

(31) Moser, J. Periodic Orbits Near an Equilibrium and a Theorem by Alan Weinstein. *Commun. Pure Appl. Math.* **1976**, *29*, 727-747.

(32) Light, J. C.; Carrington, T. J. Discrete-Variable Representations and Their Utilization. *Adv. Chem. Phys.* **2000**, *114*, 263-310.

(33) Lin, S. Y.; Guo, H. Quantum State-to-State Cross Sections for Atom-Diatom Reactions: A Chebyshev Real Wave Packet Approach. *Phys. Rev. A: At., Mol., Opt. Phys.* **2006**, *74*, 022703.

(34) Chen, R.; Guo, H. Extraction of Resonances via Wave Packet Propagation in Chebyshev Order Domain: Collinear H + H₂, Scattering. *Chem. Phys. Lett.* **1996**, *261*, 605-611.

(35) Chen, R.; Guo, H. Evolution of Quantum System in Order Domain of Chebyshev Operator. *J. Chem. Phys.* **1996**, *105*, 3569-3578.

(36) Pechukas, P. Transition State Theory. *Annu. Rev. Phys. Chem.* **1981**, *32*, 159-177.

(37) Wiggins, S. *Normally Hyperbolic Invariant Manifolds in Dynamical Systems*; Springer-Verlag, New York, 1994.

(38) Wiesenfeld, L. Geometry of Phase Space Transition States: Many Dimensions, Angular Momentum. *Adv. Chem. Phys.* **2005**, *130A*, 217-265.

(39) Mauguière, F. A. L.; Collins, P.; Kramer, Z. C.; Carpenter, B. K.; Ezra, G. S.; Farantos, S. C.; Wiggins, S. Phase Space Barriers and Dividing Surfaces in the Absence of Critical Points of the Potential Energy: Application to Roaming in Ozone. *J. Chem. Phys.* **2016**, *144*, 054107.

(40) Ezra, G. S.; Waalkens, H.; Wiggins, S. Microcanonical Rates, Gap Times, and Phase Space Dividing Surfaces. *J. Chem. Phys.* **2009**, *130*, 164118.

(41) Consider a potential energy function $V = V(q_1, ..., q_n)$ that is a function of $n$ coordinates {$q_k$}. (Coordinates describing translation and rotation are excluded.) At a nondegenerate critical point of $V$, where $\partial V / \partial q_k = 0$, $k = 1, ..., n$, the Hessian matrix $\partial^2 V / (\partial q_i \partial q_j)$ has $n$ nonzero eigenvalues. The index of the critical point is the number of negative eigenvalues.

(42) Weiss, J.; Hauschildt, J.; Grebenshchikov, S. Yu.; Düren, R.; Schinke, R.; Koput, J.; Stamatiadis, S.; Farantos, S. C. Saddle-Node Bifurcations in the Spectrum of HOCl. *J. Chem. Phys.* **2000**, *112*, 77-93.