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

: 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, . . . , 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.


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, . . . , 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.

Theory and Methods
We use Jacobi coordinates (r, R, θ) 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 r and R vectors. These definitions are shown in Figure 1. Calculations are performed considering total angular momentum null ( J = j + l = 0), with j being the angular momentum of Br 2 and 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 where V Ne,Br 2 (r, R, θ) is the van der Waals interaction between the Ne atom and the Br 2 molecule, and V Br 2 is the Morse's potential between Br-Br atoms. The quantities µ Br 2 = m −1 Br 2 + m −1 Br 2 and µ NeBr 2 = m −1 Ne + (m Br 2 + m Br 2 ) −1 stand for inverse of reduced masses for the Br 2 and NeBr 2 , respectively. The values D 1 , D 2 , α 1 , α 2 ,r andR 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 θ. 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 |ψ(t) .
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.

TSH Method in the Diabatic Representation
For our triatomic system, we can write the Hamiltonian as follow: 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, θ) 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, and Finally, H int (r, R, θ) has the form If we do a variable change |v = |ψ v (r) to simplify notation, the averaged Hamiltonian can be written as: whereĤ BC operator corresponds to the function H Br 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 θ, 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 v |H J=0 |v .
The equations motionscover on the v th surface are defined as follows (taking into account that |v =|v , v|H J=0 |v ): The state vector |ψ(t) describing the vibration of Br 2 is obtained by solving the time-dependent Schrödinger equation |ψ(t) is written as 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: Replacing H q (r)|v = E v |v , multiplying v | on the left, and taking into account ∑ v v |v = δ vv then yields, respectively: and The and for each vibrational level we determine the population as follows: where ρ vv (t) is the density matrix. The transition probabilities from the current state v to all other states v = v during the time interval ∆t are computed using the surface hopping probability (see Appendix A for more detail): 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 where P v and P v are the classical momenta which correspond to the classical coordinate R after and before the transition, respectively. The value of γ vv is obtained by imposing the total energy conservation after the transition. In this way, the angular momentum is also conserved. Energy conservation imposed using: and taking a vv = ( R v |H int (r,R,θ)|v ) 2 2µ NeBr 2 , we obtain: with the energy conservation satisfying: then there is not a real solution for this equation and the hop cannot occur. In this case, it is called a frustrated hop.
the hop can occur, and the rescaling factor (γ vv ) is computed as: when a hop occurs, we reset the wave function employing the "instantaneous decoherence" (ID) approach [26] •

Treatment of Frustrated Hop
When a frustrated hop occurs, we activate "∇V" prescription [28]. Specifically, when a frustrated hop is encountered, the following quantities are computed: where p is the nuclear momentum of the trajectory and "∇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.

1.
p h F h ≥ 0 the momentum keeps its sign; 2.
p h F h < 0 the momentum changes its sign.

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) . . . 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) . . . Ne VP , and a longer-lived contribution coming from the IVR process denoted as Br 2 (v − i) . . . Ne IVR . Therefore, only the sum is taken into account to fit the kinetic rate constants. We include the Br 2 (v − i) . . . 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: where:

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 θ 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 θ 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 × 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 coefficients are the weight of each vibrational level (if v = v ) during the dynamics of the system and c v c * v (if v = 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.

Results and Discussion
The first two columns of Table 1 report the rate constants obtained by fitting the kinetic model to the TSH results for lifetimes and intermediate state dynamics (see Equation (31a)) using  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 E rot 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. Percent % ν 0 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. Experimental results correspond to [30].
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][17][18][19][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 (∆v = −1, −2, −3...).

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  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 ± 1. In that case, we must consider fewer surfaces for dynamics than before.

Appendix A
We suppose that we have a quasiclassical system where the vibrational population for each state in the time t are determines for diagonal elements of ρ vv . For an ensemble of N trajectories that propagate simultaneously, the number of trajectories in the state v will be: In t + ∆t the ocupation state will change to: Assuming that N v (t)>N v (t + ∆t), the minimum number of transitions neccessary for this change of occupation will be: N v (t) − N v (t + ∆t), hops from state |v to another ones, and 0 hops from any state to |v The hopping probability out from state |v is: where (v = v) y ρ v v is the state where is has growing the population, Replacing the last expression for (A3), we get: Taking into account (20) and Re(ic v (t)c * v (t)) = −Im(c v (t)c * v (t)) Changing the notation then gives • if g v→v < 0, then g v→v = 0.
To determine whether a hop from the |v surface is realized, we chose a random number 0 < η < 1. We use a uniform distribution.