Study of the Vibrational Predissociation of the NeBr₂ Complex by Computational Simulation Using the Trajectory Surface Hopping Method

Abstract

The vibrational predissociation of NeBr${}_2$ has been studied using a variety of theoretical and experimental methods, producing a large number of results. It is therefore a useful system for comparing different theoretical methods. Here, we apply the trajectory surface hopping (TSH) method that consists of propagating the dynamics of the system on a potential energy surface (PES) corresponding to quantum molecular vibrational states with possibility of hopping towards other surfaces until the van der Waals bond dissociates. This allows quantum vibrational effects to be added to a classical dynamics approach. We have also incorporated the kinetic mechanism for a better compression of the evolution of the complex. The novelty of this work is that it allows us to incorporate all the surfaces for ($v=16,17,\dots ,29$) into the dynamics of the system. The calculated lifetimes are similar to those previously reported experimentally and theoretically. The rotational distribution, the rotational energy and j${}_{max}$ are in agreement with other works, providing new information for this complex.

1. Introduction

The field of the chemical physics is very broad. Among the main goals of the field is to study the properties and dynamics of molecular systems, including intra-and inter-molecular energy transfer processes leading to dissociation of excited molecules. This goal has been achieved via a large number of theoretical and experimental studies. Time resolved spectroscopic and pump-probe methods, in both the frequency and the time domains, have been especially useful for providing data to test theoretical methods [1,2]. These methods have been applied to the vibrational predissociation of an extensive variety of van der Waals (vdW) complexes composed of three or more atoms with a range of bond energies, atomic masses, and vibrational frequencies. Due to the weakness of the vdW interactions in these systems, the constituents retain their chemical integrity upon complex formation so the energy transfer mechanism can be easily identified and studied at the state-to-state level [3,4,5,6,7,8,9,10,11,12,13,14,15,16]. The NeBr${}_2$ van der Waals molecule has been particularly useful for these studies since the dynamics is on the boarder line for which classical and quantum methods are comparably useful. This allows us to investigate in particular the applicability of the trajectory surface hopping (TSH) method [8,14,15,16]. The TSH method has been widely used [17,18,19,20,21,22,23], demonstrating its validity and efficiency for a variety of molecular dynamics problems. For this study of the vibrational predissociation of NeBr${}_2$, the diabatic potential energy surfaces are formed by the interaction of the Ne atoms with the $v=16,17,\dots ,29$ Br${}_2$ vibrational levels. The couplings between surfaces are provided by the van der Waals potential. The TSH results obtained for NeBr${}_2$ vibrational predissociation are compared to previous work using other methods [8,9,10]. We also implement the “kinetic mechanism” [15] to interpret the results of the TSH simulation. This method considers two mechanism for transferring energy from the vibration of the diatom to the van der Waals modes. The first mechanism corresponds to a direct vibrational predissociation (VP) transfer of the dissociative coordinate. The second mechanism involves preliminary energy transfer to the non-dissociative van der Waals modes, followed by intramolecular vibrational redistribution (IVR). After IVR has taken place, the cluster is cooled by expelling the rare gas atoms, a process called IVR-evaporative cooling (EC). The values of the kinetic rate constants which characterize these elementary steps are determined by fitting the results of the TSH simulation to the analytical expressions for the time evolution of the NeBr${}_2$ concentration. The simplicity of our treatment and its relatively low computational cost will allow it to be extended to systems with more degrees of freedom (e.g., more rare gas atom), thus offering an attractive and different alternative to purely classical treatments.

The paper is organized as follows: in Section 2, we discuss the procedure of the TSH method as well as the computational details. We also describe our implementation of the kinetic mechanism. In Section 3, results are presented and discussed through figures and tables. Section 4 summarizes the conclusions that can be drawn from this work. In the Appendix A, we describe our procedure to compute the transition probabilities in considerable detail.

2. Theory and Methods

We use Jacobi coordinates $(r,R,\theta )$ to describe the NeBr${}_2$ triatomic complex, with r being the bond length of Br${}_2$, R the intermolecular distance from the Ne atom to the center of mass of the dihalogen and the angle between $\overrightarrow{r}$ and $\overrightarrow{R}$ vectors. These definitions are shown in Figure 1. Calculations are performed considering total angular momentum null ( $\overrightarrow{J}=\overrightarrow{j}+\overrightarrow{l}=\overrightarrow{0}$), with $\overrightarrow{j}$ being the angular momentum of Br${}_2$ and $\overrightarrow{l}$ the orbital angular momentum, a well justified constraint while studying photodissociation events.

The classical Hamilton function corresponding to this three degree of freedom model of the complex is written as

$$H^{J=0}=H_{{\mathrm{Br}}_2}+\frac{1}{2{\mu }_{\mathrm{Ne}{\mathrm{Br}}_2}}(p_R^2+\frac{p_{\theta }^2}{R^2})+V_{int}(r,R,\theta ),$$

(1)

where

$$H_{{\mathrm{Br}}_2}=\frac{1}{2{\mu }_{{\mathrm{Br}}_2}}(p_r^2+\frac{p_{\theta }^2}{r^2})+V_{{\mathrm{Br}}_2}(r),$$

(2)

$V_{\mathrm{Ne},{\mathrm{Br}}_2}(r,R,\theta )$ is the van der Waals interaction between the Ne atom and the Br${}_2$ molecule,

$$V_{Br2}(r)=D_1(e^{[-2{\alpha }_1(r-\overline{r})]}-2e^{[-{\alpha }_1(r-\overline{r})]}),$$

(3)

$$V_{\mathrm{Ne}Br2}(r,R,\theta )=V_R(R_1)+V_R(R_2),$$

(4)

$$\begin{gathered} V_R(R)=D_2(e^{[-2{\alpha }_2(R-\overline{R})]}-2e^{[-{\alpha }_2(R-\overline{R})]}) \\ {R}_{1,2}^2=R^2+r^2/4\pm rRcos(\theta ), \end{gathered}$$

and $V_{B{r}_2}$ is the Morse’s potential between Br-Br atoms. The quantities ${\mu }_{{\mathrm{Br}}_2}=m_{{\mathrm{Br}}_2}^{-1}+m_{{\mathrm{Br}}_2}^{-1}$ and ${\mu }_{\mathrm{Ne}{\mathrm{Br}}_2}=m_{\mathrm{Ne}}^{-1}+{(m_{{\mathrm{Br}}_2}+m_{{\mathrm{Br}}_2})}^{-1}$ stand for inverse of reduced masses for the Br${}_2$ and NeBr${}_2$, respectively. The values $D_1$, $D_2$, ${\alpha }_1$, ${\alpha }_2$, $\overline{r}$ and $\overline{R}$ are from [24].

For the quantum treatment, we have used as quantum coordinate the Br${}_2$ vibration $(r)$ and the other variables are treated as classical coordinates; for our case these are R and $\theta$. A potential energy surface (PES) is defined by a single state of the quantum degree of freedom (r). It is on this surface that the trajectories for the classical degrees of freedom evolve. Quantum transitions are modeled by hops of the trajectories from one surface to another. These are governed by the evolution of the multicomponent time-dependent vibrational wave function $|\psi (t)\rangle$.

As we have commented before, our surfaces are diabatic (see Figure 2) and the couplings are provided by the van der Waals potential. This is different from the usual trajectory surface hopping TSH treatment where the transitions occur between electronic adiabatic surfaces, in well defined regions of avoided crossings.

2.1. TSH Method in the Diabatic Representation

For our triatomic system, we can write the Hamiltonian as follow:

$$H(r,R,\theta )=H_q(r)+H_{el}(R)+H_{int}(r,R,\theta )$$

(5)

where $H_q(r)$ is the quantum part depending on the vibrational coordinate r of Br${}_2$, $H_{cl}(R)$ is the Hamiltonian describing the classical degrees of freedom and $H_{int}(r,R,\theta )$ is the interaction operator which couples the quantum and classical degrees of freedom.

In our case, $H_q(r)$ describes the vibrational motion of Br${}_2$ in the B state,

$$H_q(r)=-\frac{{\hslash }^2}{2{\mu }_{B{r}_2}}\frac{{\partial }^2}{\partial r^2}+V_{B{r}_2}(r)$$

(6)

and

$$H_{cl}(R)=\frac{p_R^2}{2{\mu }_{NeB{r}_2}}+\frac{p_{\theta }^2}{2{\mu }_{NeB{r}_2}{R}^2}$$

(7)

Finally, $H_{int}(r,R,\theta )$ has the form

$$H_{int}(r,R,\theta )=\frac{p_{\theta }^2}{2{\mu }_{B{r}_2}{r}^2}+V_{NeB{r}_2}(r,R,\theta )$$

(8)

If we do a variable change $|v\rangle =|{\psi }_v(r)\rangle$ to simplify notation, the averaged Hamiltonian can be written as:

$$\begin{array}{ccc}\hfill \langle v^{\prime }|H^{J=0}|v\rangle & =& (E_v+\frac{p_R^2}{2{\mu }_{NeB{r}_2}}){\delta }_{v^{\prime }v}+\langle v^{\prime }|V_{NeB{r}_2}(r,R,\theta )|v\rangle +\frac{p_{\theta }^2}{2{\mu }_{B{r}_2}}\langle v^{\prime }|\frac{1}{r^2}|v\rangle +\hfill \\ & +& \frac{p_{\theta }^2}{2{\mu }_{NeB{r}_2}{R}^2}{\delta }_{v{v}^{\prime }},\hfill \end{array}$$

where ${\delta }_{v{v}^{\prime }}=$ 1 if $v=v^{\prime }$ and 0 if $v\ne v^{\prime }$ is the Kronecker’s delta and $E_v$ is the vibrational energy for Br${}_2$ obtained by solving

$${\widehat{H}}_{BC}{\psi }_v(r)=E_v{\psi }_v(r),$$

(9)

where ${\widehat{H}}_{BC}$ operator corresponds to the function $H_{B{r}_2}(r)$ defined in Equation (6).

We have used Bode’s method [25] for the integration of these equations. This subroutine is very important because we have the Morse potential depending of $r,R$ and $\theta$, and we want to know the averaged effect of r for a vibrational level. It is the source of the PES and the coupling among them is due to non diagonal elements $\langle v^{\prime }|H^{J=0}|v\rangle$.

The equations motionscover on the $v^{th}$ surface are defined as follows (taking into account that $|v^{\prime }\rangle$=$|v\rangle$, $\langle v|H^{J=0}|v\rangle$):

$$\frac{dR}{dt}=\frac{\partial \langle v|H^{J=0}|v\rangle }{\partial p_R}=\frac{p_R}{{\mu }_{\mathrm{Ne}{\mathrm{Br}}_2}}$$

(10)

$${\displaystyle \frac{d\theta }{dt}=\frac{\partial \langle v|H^{J=0}|v\rangle }{\partial p_{\theta }}=p_{\theta }\left(\frac{\langle v|\frac{1}{r^2}|v\rangle }{{\mu }_{B{r}_2}}+\frac{1}{{\mu }_{NeB{r}_2}{R}^2}\right)}$$

(11)

$$\frac{d{p}_R}{dt}=-\frac{\partial \langle v|H^{J=0}|v\rangle }{\partial R}=\frac{p_{\theta }^2}{{\mu }_{NeB{r}_2}{R}^3}-\langle v|\frac{\partial V_{NeB{r}_2}(r,R,\theta )}{\partial R}|v\rangle$$

(12)

$$\frac{d{p}_{\theta }}{dt}=-\frac{\partial \langle v|H^{J=0}|v\rangle }{\partial \theta }=-\langle v|\frac{\partial V_{NeB{r}_2}(r,R,\theta )}{\partial \theta }|v\rangle .$$

(13)

The state vector $|\psi (t)\rangle$ describing the vibration of Br${}_2$ is obtained by solving the time-dependent Schrödinger equation

$$i\hslash \frac{\partial |\psi (t)\rangle }{\partial t}=[H_q(r)+H_{int}(r,R,\theta )]|\psi (t)\rangle .$$

(14)

$|\psi (t)\rangle$ is written as

$$|\psi (t)\rangle =\sum _v{c}_v(t)e^{-i{E}_{v}t/\hslash }|v\rangle$$

(15)

where the sum is over all v states of $H_q(r)$ with energy $E_v$, $c_v(t)$. The complex variable indicates the amplitude of each vibrational level over the total wave function. This is known as the semiclassical expansion of the electronic wave function. Replacing (Equation (15)) in (Equation (14)), we obtain:

$$i\hslash \frac{\partial }{\partial t}\sum _v{c}_v(t)e^{-i{E}_{v}t/\hslash }|v\rangle =\sum _v{c}_v(t)e^{-i{E}_{v}t/\hslash }[H_q(r)+H_{int}(r,R,\theta )]|v\rangle .$$

(16)

Replacing $H_q(r)|v\rangle =E_v|v\rangle$, multiplying $\langle v^{\prime }|$ on the left, and taking into account ${\displaystyle \sum _v}\langle v^{\prime }|v\rangle ={\delta }_{v{v}^{\prime }}$ then yields, respectively:

$${\displaystyle i\hslash e^{-i{E}_{v}t/\hslash }\left({\dot{c}}_v(t)-\frac{i{E}_v}{\hslash }{c}_v(t)\right)=c_v(t)E_v{e}^{-i{E}_{v}t/\hslash }+\sum _{v^{\prime }}{c}_{v^{\prime }}(t)e^{-i{E}_{v^{\prime }}t/\hslash }\langle v^{\prime }|H_{int}(r,R,\theta )|v\rangle ,}$$

(17)

$$i\hslash e^{-i{E}_{v}t/\hslash }{\dot{c}}_v(t)+c_v(t)E_v{e}^{-i{E}_{v}t/\hslash }=c_v(t)E_v{e}^{-i{E}_{v}t/\hslash }+\sum _{v^{\prime }}{c}_{v^{\prime }}(t)e^{-i{E}_{v^{\prime }}t/\hslash }\langle v^{\prime }|H_{int}(r,R,\theta )|v\rangle ,$$

(18)

and

$$i\hslash {\dot{c}}_v(t)=\sum _{v^{\prime }}{c}_{v^{\prime }}(t)e^{i(E_v-E_{v^{\prime }})t/\hslash }\langle v^{\prime }|H_{int}(r,R,\theta )|v\rangle .$$

(19)

The $c_v(t)$ coefficients satisfy

$$\dot{c_v}(t)=-\frac{i}{\hslash }\sum _{v^{\prime }}\langle v^{\prime }|H_{int}(r,R,\theta )|v\rangle e^{[i(E_v-E_{v^{\prime }})t/\hslash ]}{c}_{v^{\prime }}(t)$$

(20)

and for each vibrational level we determine the population as follows:

$${\rho }_{v{v}^{\prime }}(t)=c_v(t)c_{v^{\prime }}^{*}(t),$$

(21)

where ${\rho }_{v{v}^{\prime }}(t)$ is the density matrix.

The transition probabilities from the current state v to all other states $v^{\prime }\ne v$ during the time interval $\Delta t$ are computed using the surface hopping probability (see Appendix A for more detail):

$$g_{v\to v^{\prime }}=\frac{2\Delta t}{\hslash }\frac{{\displaystyle Im\left(c_v^{*}(t)c_{v^{\prime }}(t)e^{[i(E_v-E_{v^{\prime }})t/\hslash ]}\right)}}{|c_v{(t)|}^2}\langle v^{\prime }|H_{int}(r,R,\theta )|v\rangle .$$

(22)

The initial conditions for the classical trajectories are selected randomly for a total energy corresponding to the zero-point of the complex NeBr${}_2$, for a particular vibrational state of the Br${}_2$ molecule. The component of the momentum that is parallel to the quantum state coupling vector is only taken into account to adjust, in order to conserve, the total energy [26,27]. In our case, quantum transitions occur between diabatic surfaces, defined for each vibrational level Br${}_2(B)$. These surfaces are coupled by the NeBr${}_2$ potential. Momenta are adjusted during a surface hop using

$${\mathbf{P}}_{v^{\prime }}={\mathbf{P}}_v-{\gamma }_{v{v}^{\prime }}{\nabla }_R\langle v|H_{int}(r,R,\theta )|v^{\prime }\rangle$$

(23)

where ${\mathbf{P}}_{v^{\prime }}$ and ${\mathbf{P}}_v$ are the classical momenta which correspond to the classical coordinate R after and before the transition, respectively. The value of ${\gamma }_{v{v}^{\prime }}$ is obtained by imposing the total energy conservation after the transition. In this way, the angular momentum is also conserved. Energy conservation imposed using:

$$\frac{P_{v^{\prime }}^2}{2{\mu }_{NeB{r}_2}}+\langle v^{\prime }|H_{int}(r,R,\theta )|v^{\prime }\rangle +E_{v^{\prime }}=\frac{P_v^2}{2{\mu }_{NeB{r}_2}}+\langle v|H_{int}(r,R,\theta )|v\rangle +E_v$$

(24)

Taking into account $|{\overrightarrow{P}}_{v^{\prime }}(R)|=P_{v^{\prime }}$, $|{\overrightarrow{P}}_v(R)|=P_v$, replacing

$$E_0=\langle v|H_{int}(r,R,\theta )|v\rangle +E_v$$

(25)

$$E_f=\langle v^{\prime }|H_{int}(r,R,\theta )|v^{\prime }\rangle +E_{v^{\prime }}$$

(26)

and taking $a_{v{v}^{\prime }}=\frac{({▽}_R\langle v^{\prime }|H_{int}{(r,R,\theta )|v\rangle )}^2}{2{\mu }_{NeB{r}_2}}$ and $b_{v{v}^{\prime }}=\frac{P_v{▽}_R\langle v^{\prime }|H_{int}(r,R,\theta )|v\rangle }{{\mu }_{\mathrm{Ne}{\mathrm{Br}}_2}}$, we obtain:

$$\frac{P_{v^{\prime }}^2}{2{\mu }_{NeB{r}_2}}-\frac{P_v^2}{2{\mu }_{NeB{r}_2}}=a_{v{v}^{\prime }}{\gamma }_{v{v}^{\prime }}^2-b_{v{v}^{\prime }}{\gamma }_{v{v}^{\prime }}$$

(27)

with the energy conservation satisfying:

$$a_{v{v}^{\prime }}{\gamma }_{v{v}^{\prime }}^2-b_{v{v}^{\prime }}{\gamma }_{v{v}^{\prime }}-(E_0-E_f)=0$$

(28)

• If $b_{v{v}^{\prime }}^2+4a_{v{v}^{\prime }}(E_0-E_f)<0$, then there is not a real solution for this equation and the hop cannot occur. In this case, it is called a frustrated hop.
• If $b_{v{v}^{\prime }}^2+4a_{v{v}^{\prime }}(E_0-E_f)\ge 0$, the hop can occur, and the rescaling factor (${\gamma }_{v{v}^{\prime }}$) is computed as:

$${\gamma }_{v{v}^{\prime }}=\frac{b_{v{v}^{\prime }}+\sqrt{b_{v{v}^{\prime }}^2+4a_{v{v}^{\prime }}(E_0-E_f)}}{2a_{v{v}^{\prime }}}if{b}_{v{v}^{\prime }}<0$$

(29)

$${\gamma }_{v{v}^{\prime }}=\frac{b_{v{v}^{\prime }}-\sqrt{b_{v{v}^{\prime }}^2+4a_{v{v}^{\prime }}(E_0-E_f)}}{2a_{v{v}^{\prime }}}if{b}_{v{v}^{\prime }}\ge 0$$

(30)

when a hop occurs, we reset the wave function employing the “instantaneous decoherence” (ID) approach [26]

• $c_v=0$∀$v\ne v^{\prime }$
• $c_{v^{\prime }}=1$

2.2. Treatment of Frustrated Hop

When a frustrated hop occurs, we activate “$\nabla V$” prescription [28]. Specifically, when a frustrated hop is encountered, the following quantities are computed:

• $p_h=$$\overrightarrow{p}·\overrightarrow{h}$
• $F_h=-\nabla V_v·\overrightarrow{h}$

where $\overrightarrow{p}$ is the nuclear momentum of the trajectory and “$\nabla V_v$” is the gradient of the target vibrational state v, $p_h$ and $F_h$ are the projection of the nuclear momentum and the force of the target vibrational state along the hopping vector h, respectively. If $p_h$ and $F_h$ have the same sign, the target vibrational state accelerates the trajectory along h. Otherwise, if the two quantities have opposite signs, the target vibrational state causes a delay in the trajectory of the Ne. For that, we use the follow criterion for frustrated hop.

• $p_h{F}_h$$\ge 0$ the momentum keeps its sign;
• $p_h{F}_h$$<0$ the momentum changes its sign.

2.3. Kinetic Mechanism

The kinetic mechanism allows us to understand the path followed by our system as it relaxes and dissociates (see Figure 3). We consider Br${}_2(v-i)\dots Ne$ (intermediate state detected in the TSH simulation) as a sum of two contributions: a short-lived contribution coming from the VP process, which we denote by Br${}_2(v-i)\dots$ Ne${}_{VP}$, and a longer-lived contribution coming from the IVR process denoted as Br${}_2(v-i)\dots$ Ne${}_{IVR}$. Therefore, only the sum is taken into account to fit the kinetic rate constants. We include the Br${}_2(v-i)\dots$ Ne${}_{VP}$ intermediates in the direct VP process, dividing each direct VP step into two processes characterized by the rate constants $k_{VPa}$ and $k_{VPb}$ for the loss of the first vibrational quantum, and $k_{VP2a}$ and $k_{VP2b}$ for the loss of the second one.

We fit the data obtained in the simulation by using the next procedure:

$$[{\mathrm{Br}}_2(v)\dots \mathrm{Ne}]=[{\{{\mathrm{Br}}_2(v)\dots \mathrm{Ne}\}}_{ivr}+{\{{\mathrm{Br}}_2(v)\dots \mathrm{Ne}\}}_{vp}]$$

(31a)

$$\left[{\{{\mathrm{Br}}_2(v)\dots \mathrm{Ne}\}}_{ivr}\right]={\displaystyle \frac{k_{ivr}}{k_v}}{e}^{-k_{v}t}$$

(31b)

$$\left[{\{{\mathrm{Br}}_2(v)\dots \mathrm{Ne}\}}_{vp}\right]={\displaystyle \frac{k_{vpa}}{k_v}}{e}^{-k_{v}t}$$

(31c)

$$[{\mathrm{Br}}_2(v-1)\dots \mathrm{Ne}]=[{\{{\mathrm{Br}}_2(v-1)\dots \mathrm{Ne}\}}_{ivr}+{\{{\mathrm{Br}}_2(v-1)\dots \mathrm{Ne}\}}_{vp}]$$

(32a)

$$\left[{\{{\mathrm{Br}}_2(v-1)\dots \mathrm{Ne}\}}_{ivr}\right]={\displaystyle \frac{k_{ivr}}{k_1-k_v}}(e^{-k_{v}t}-e^{-k_{1}t})$$

(32b)

$$\left[{\{{\mathrm{Br}}_2(v-1)\dots \mathrm{Ne}\}}_{vp}\right]={\displaystyle \frac{k_{vpa}}{k_v-k_{vpb}}}(e^{-k_{vpb}t}-e^{-k_{v}t})$$

(32c)

$$\left[{\mathrm{Br}}_2(v-1)\right]=[{\{{\mathrm{Br}}_2(v-1)\}}_{ivr}+{\{{\mathrm{Br}}_2(v-1)\}}_{vp}]$$

(33a)

$${\displaystyle \left[{\{B{r}_2(v-1)\}}_{ivr}\right]={\displaystyle \frac{k_{ivr}{k}_{ec}}{(k_1-k_v)k_1{k}_v}}\left(k_1(1-e^{-k_{v}t})-k_v(1-e^{-k_{1}t})\right)}$$

(33b)

$${\displaystyle \left[{\{B{r}_2(v-1)\}}_{vp}\right]={\displaystyle \frac{k_{vpa}}{(k_v-k_{vpb})k_v}}\left(k_v(1-e^{-k_{vpb}t})-k_{vpb}(1-e^{-k_{v}t})\right)}$$

(33c)

$$[{\mathrm{Br}}_2(v-2)\dots \mathrm{Ne}]=[{\{{\mathrm{Br}}_2(v-2)\dots \mathrm{Ne}\}}_{ivr}+{\{{\mathrm{Br}}_2(v-2)\dots \mathrm{Ne}\}}_{vp}]$$

(34a)

$$\begin{array}{cc}\hfill \left[{\{B{r}_2(v-2)\dots Ne\}}_{ivr}\right]& ={\displaystyle \frac{k_{ivr}{k}_{ivr2}}{(k_1-k_v)(k_2-k_1)(k_2-k_v)}}\{(k_2-k_1)e^{-k_{1}t}+\\ & +(k_v-k_2)e^{-k_{1}t}+(k_1-k_v)e^{-k_{2}t}\}\end{array}$$

(34b)

$$\begin{array}{cc}\hfill \left[{\{B{r}_2(v-2)\dots Ne\}}_{vp}\right]& ={\displaystyle \frac{k_{ivr}{k}_{vp2a}}{(k_1-k_v)(k_1-k_{vp2b})(k_{vp2b}-k_v)}}\{(k_1-k_{vp2b})e^{-k_{v}t}+\\ & +(k_{vp2b}-k_v)e^{-k_{1}t}+(k_v-k_1)e^{-k_{vp2b}t}\}\end{array}$$

(34c)

$$\left[{\mathrm{Br}}_2(v-2)\right]=[{\{{\mathrm{Br}}_2(v-2)\}}_{ivr}+{\{{\mathrm{Br}}_2(v-2)\}}_{vp}]$$

(35a)

$$\begin{array}{cc}\hfill \left[{\{B{r}_2(v-2)\}}_{ivr}\right]& ={\displaystyle \frac{k_{ivr}{k}_{ivr2}{k}_{ec2}}{(k_1-k_v)(k_2-k_1)(k_2-k_v)(k_1{k}_v{k}_2)}}\{k_2{k}_1(k_2-k_1)(1-e^{-k_{v}t})\\ & +k_v{k}_2(k_v-k_2)(1-e^{-k_{1}t})+k_1{k}_v(k_1-k_v)(1-e^{-k_{2}t})\}\end{array}$$

(35b)

$$\begin{array}{cc}\hfill \left[{\{B{r}_2(v-2)\}}_{vp}\right]& ={\displaystyle \frac{k_{ivr}{k}_{vp2a}}{(k_1-k_v)(k_1-k_{vp2b})(k_{vp2b}-k_v)(k_1{k}_v)}}\{k_1{k}_{vp2b}(k_1-k_{vp2b})(1-e^{-k_{v}t})\\ & +k_{vp2b}{k}_v(k_{vp2b}-k_v)(1-e^{-k_{1}t})+k_v{k}_1(k_v-k_1)(1-e^{-k_{vp2b}t})\}\end{array}$$

(35c)

where:

• $k_v=k_{vpa}+k_{ivr}$
• $k_1=k_{vp2a}+k_{ivr2}+k_{ec}$
• $k_2=k_{ec2}+k^{\prime }$

2.4. Computational Details

In the methodology involved, two important stages contribute:

• First stage

This step consists of propagating the dynamics of the system by evolving classically the nuclear motion on the potential energy surface. For this, employ the Adams Bashfort method [25], initiating with the method of Runge Kutta 4th order. We obtain, from the system of Equations (10)–(13), the momentum and coordinate of the Ne, the $\theta$ angle and the angular momentum of the system in Jacobi coordinates (see Figure 1). These coordinates are defined as follows: r is the distance between the diatomic constituents, R is the distance of the noble gas to the mass of the chemical bound and $\theta$ is the angle between r and R. In the method, some physical considerations of importance are imposed for obtaining a consistent result, such as the conservation of total energy. We consider that our system is dissociated beyond a certain maximum distance. For this, we take R${}_{max}$$=10$Å and a maximum time which the system will remain bounded, t${}_{max}$$=600$ ps. Taking into account the above, we use an integration step equal to $0.1$ ps, which leads to $6\times {10}^6$ integration cycles and ensuring a total energy conservation error of less than ${10}^{-8}$ cm${}^{-1}$.

Simultaneously, in the other simulation thread, the Equation (20) is integrated, obtaining the $c_v$ coefficients. The $c_v{c}_{v^{\prime }}^{*}$ coefficients are the weight of each vibrational level (if $v=v^{\prime }$) during the dynamics of the system and $c_v{c}_{v^{\prime }}^{*}$ (if $v\ne v$) indicate the coherence between states.

In order to apply the population conservation, the sum of the populations has to be equal the unity. We also incorporate a method for the hop decision (it is known as the Fewest Switches algorithm, see Appendix A). Another important aspect is the rescaling of the momentum to preserve the total energy of the system when the hop occurs.

To obtain the average interaction potentials of van der Waals, the Bode method is implemented. It is very important to define the average potentials and crossings between curves (see Figure 2).

• Second stage

This stage consists of performing a fitting to initial, intermediate and final populations during the simulation. The complexity of doing this is that the parameters are shared and therefore obtaining these depends on the statistical behavior for each initial vibrational level. We followed the steps from the reference [29], which is a generalization of the Marquardt method for multiple equations and shared parameters.

3. Results and Discussion

The first two columns of

Table 1. Rate constants (in $n{s}^{-1}$) obtained by fitting the trajectory surface hopping (TSH) simulation results.

${\mathit{v}}_0$	${\mathit{k}}_{\mathbf{ivr}}$	${\mathit{k}}_{\mathbf{vpa}}$	${\mathit{k}}_{\mathit{v}}$	${\mathit{k}}_1$	${\mathit{k}}_{\mathbf{vpb}}$	${\mathit{k}}_{\mathbf{ec}}$	${\mathit{k}}_{\mathbf{ivr}2}$	${\mathit{k}}_2$	${\mathit{k}}_{\mathbf{vp}2\mathit{b}}$	${\mathit{k}}_{\mathbf{ec}2}$	${\mathit{k}}_{\mathbf{vp}2\mathit{a}}$
16	0.024	4.456	4.600	81,019	8,050,000	1,704,730	51,948	2,523,430	604.148	361.258	24,245
17	0.110	5.716	5.887	60,896	882.070	733.634	15,961	846.233	150.962	1948	153,785
18	0.055	5.410	5.620	20,124	2935.5	7713	28,970	818.535	12.646	3335	23,207
19	1.484	12.805	13.991	63,190	223.312	1907	15,356	217.080	499.413	496.263	11,185
20	2.742	14.914	16.970	55,278	185	734.147	7684	171.379	18.734	769.698	1929
21	11.999	35.608	41.82	2475	162.684	14.016	211.989	123.330	169.610	197.089	791.290
22	3.934	14.956	18.821	5723	228.020	26.9789	367.425	39.847	932.066	145.867	3575
23	4.685	16.540	20.996	5263	168.904	29.035	1023	118.771	268.294	102.744	3339
24	3.966	13.345	17.937	138.804	237.482	1.030	157.180	1911	4412	551.141	86.110
25	10.397	19.055	29.327	149.336	142.794	12.052	132.472	1821	1507	261.619	86.482
26	11.937	16.686	29.887	100.351	126.820	0.157	118.667	3502	6501	1001.190	48.893
27	115.150	0.153	105.520	76.882	0.110	25.249	50.949	1826	692,718	34.391	32.183
28	69.295	2.478	71.940	85.525	1250	2.920	1604	12,152	16,751	444.040	5.350
29	133.750	3.642	127.640	81.460	7.287	2.087	1536	12,276	11,827	431.941	5.684

report the rate constants obtained by fitting the kinetic model to the TSH results for lifetimes and intermediate state dynamics (see Equation (31a)) using

$$\tau ={\displaystyle \frac{1}{k_{ivr}+k_{vpa}}}$$

(36)

Table 1 shows the rate constants obtained from the fitting of the TSH results to the kinetic mechanism. From this table we get much information about the path followed and the time spent in each initial, intermediate and final state. If ($k_{vpa}>k_{vpb}$) or ($k_{vp2a}>k_{vp2b}$), the system remains more time in the intermediate state before the dissociation occurs. However, if ($k_{vpb}>k_{vpa}$) or ($k_{vp2b}>k_{vp2a}$) occurs, the complex breaks its bond faster than otherwise. Moreover, for the lower vibrational levels, we obtained larger values for $k_{vpb}$. This means that just after hopping, the molecule is dissociated. On the other hand, for the higher vibrational level there is a great probability of following an IVR process ($k_{ivr}>k_{vpa}$), and the others follow a VP process $(k_{vpa}>k_{ivr})$. In addition, there is a competition between $k_{ec}$ and $k_{ivr2}$. When $k_{ec}>k_{ivr2}$, the system is dissociated, losing a vibrational quantum number. In another case, the system is submitted to an IVR process again. In Figure 4, we show the lifetimes obtained.

In Figure 5, we show the process which the system followed until its dissociation.

In Figure 5, we can see that from vibrational level $v=21$ the loss of two quanta through the VP $v_0-2$ process is significant. This causes the system to take longer to dissociate. Particularly for vibrational levels $v=28,29$, the predominant process is IVR-EC $v_0-2$.

Through the simulation, we can know at what vibrational level the system is, and we do not need algorithms to obtain it. This is very good because the TSH method itself gives us that information. In the following figure, we present the exit channel (statistically) for each initial vibrational level.

As we can see in Figure 6 that the TSH results are in agreement with previous theoretical and experimental results. Other very important observables are $\langle E_{rot}\rangle$ and $j_{max}$. The first one is the averaged rotational energy. The second one is the Br-Br angular momentum. Both quantities are calculated when the molecule is dissociated.

As we can appreciate in Figure 7, we report our results in comparison with experimental ones and previous results.

In Figure 8, we represent the rotational distributions of Br${}_2$ after dissociation. We can see that there is a peak in the range of j = 20–24 (Figure 8c). This is an effect that we can see from $v_0=$ 21 to $v_0=$ 26, with the loss of two vibrational quantum numbers. This means that, for these vibrational levels having greater rotational excitation, the system takes a longer time to dissociate. Therefore, it does not have enough energy to break the vdW bond. They are non-reactive van der Waals modes and, as a direct consequence, the lifetimes are longer.

To analyze the conservation relation, we compute the correlation matrix for all vibrational levels under study. As we discussed before for the lowest vibrational levels ($v=$ 16–20), the fragmentation process occurs when the system loses one quantum energy. For those vibrational levels, more correlation is found for the $k_{vpa}$ and $k_{ivr}$ parameters. As we have shown in Figure 9, the parameter $k_{ec}$ shows a strong correlation with those parameters where the process involves two losses of quantum energy. Then, $k_{ec}$ plays an important role because it links the number of quantum energy losses in the fragmentation of the system ($\Delta v=-1,-2,-3...$).

4. Conclusions

We have investigated the vibrational predissociation process for the NeBr${}_2$ system by using TSH method. In our simulation, we studied a range of vibrational levels (from $v=16$ to $v=29$). We found that the larger values for $k_{vpb}$ correspond to the lower vibrational levels. This is in correspondence with the dissociated of the molecule occurs just after hopping surface. As we comment in Section 3, there is a competition between VP and IVR processes. The TSH method is a robust methodology for the study of molecular fragmentation for this kind of system. This affirmation could be justified through Figure 6, Figure 7 and Figure 8. In these figures, we showed our results in comparison with the previous theoretical results and the experimental ones. As we can appreciate, the agreements are very good.

Perspectives

This method is a powerful tool for dealing with this type of system and can be extended to more complex structures (e.g., more vdW interactions). We continue with this work and we are checking that the strongest coupling belongs to $v\pm 1$. In that case, we must consider fewer surfaces for dynamics than before.