Nonadiabatic fragmentation of H₂O⁺ and isotopomers. Wave packet propagation using ab initio wavefunctions

Abstract

The fragmentation of the water cation from its B̃ ²B₂ electronic state, allowing the participation of the X̃ ²B₁, Ã ²A₁ and C̃ ²B₁ states in the process, is simulated using the extended capabilities of the collocation GridTDSE code to account for the nonadiabatic propagation of wave packets in several potential energy surfaces connected by nonadiabatic couplings. Molecular data are calculated ab initio. Two initial wave packets are considered to reproduce two different experiments. The isotopic effect in the fragmentation of D₂O⁺ and HDO⁺ is also studied and the results show very good agreement with the experimental cleavage preference in the fragmentation of HDO⁺.

---

1 Introduction

The water cation is naturally found in several environments such as the interstellar medium,^1,2 and atmospheres of planets³ and comets.^4,5 It can be formed in collisions of electrons, ions or photons with H₂O and, given the importance of this molecule, many works have studied its ionization, and the ensuing dynamics of the remaining cation. Several experimental^6–8 and theoretical^9,10 works have addressed the photoionization of H₂O with photon energies below 21.2 eV. At these energies, the photoelectron spectrum shows three bands that correspond to the formation of H₂O⁺ in the electronic states X̃ ²B₁ (hereafter denoted X̃), Ã ²A₁ (Ã) and B̃ ²B₂ (B̃). These electronic states are essentially described by one-electron removal from the three outermost molecular orbitals (1b₁, 3a₁ and 1b₂, respectively) of H₂O. For photon energies higher than the dissociation threshold (18.1 eV), the vibrational levels of the B̃ state are predissociative and, accordingly, the third band of the H₂O photoelectron spectrum exhibits broad peaks. The ionization of H₂O above the dissociation threshold leads to the cation fragmentation into H⁺ + OH, H + OH⁺ and O⁺ + H + H. Following the ionization of H₂O, coincidence experiments^11–14 have measured the branching ratios in the formation of H⁺ + OH and H + OH⁺. Experiments of ref. 15–17 have measured the kinetic energy release distributions of the fragments formed after electron impact on water molecules, which shed light on the possible routes leading to the different fragmentation schemes.

Calculations on the H₂O⁺ system include the ab initio studies^18,19 of the reactions O + H₂⁺ → OH⁺ + H and O + H₂⁺ → OH + H⁺, and the corresponding reactions with the isotopomers D₂⁺ and HD⁺. These calculations used the quasi-classical-trajectory and real-wave-packet methods, and state-of-the-art potential energy surfaces of the X̃ and à electronic states. The influence of the Renner–Teller coupling on the O + H₂⁺ proton transfer reactions was analyzed in ref. 20.

In a previous paper,²¹ we carried out a theoretical study of the fragmentation of the water cation from the B̃ state, obtaining values for the dissociation branching ratios into OH + H⁺ and H + OH⁺ that were in excellent agreement with the experimental data of Tan et al.¹² Our calculations consisted in nuclear wave-packet propagation on three electronic surfaces (B̃, à and X̃), using a modified version of the GridTDSE code²² to reproduce the nonadiabatic quantum nuclear dynamics. The initial condition of our dynamical calculation was the Franck–Condon (FC) wave packet, obtained by vertical ionization of the ground vibrational state of H₂O. The wave-packet propagation showed that the dissociation mechanism involves nonadiabatic transitions from the state B̃ to à and X̃. The transition B̃–à takes place in a few femtoseconds in the vicinity of the conical intersection (CI) between the corresponding potential energy surfaces (PESs). The CI seam appears at angles in the interval 60° ≲ α ≲ 100°, while the minimum of the à PES is located at the linear geometry (α = 180°). As a consequence, the wave packet starts in B̃, moves towards the CI and is transferred to à with an important excitation of the bending mode. At linear geometries, the states à and X̃ become the two components of the ²Π_u state, and the ÖX̃ Renner–Teller (RT) interaction allows the population of the ground electronic state. The fragmentation can take place either by energy transfer from the bending mode into the antisymmetric stretching motion within Ã, leading to dissociation into OH + H⁺, or after the RT transition between à and X̃, leading to dissociation into OH⁺ + H.

Although the good agreement between our results and the experimental branching ratios indicates that the three-state (B̃–ÖX̃) mechanism explains the predissociation of the water cation, other possible nonadiabatic transitions must be ruled out. In our previous work,²¹ we checked that the intersystem crossing transitions from B̃ and à to the dissociative ã ⁴B₁ state are very slow compared to B̃–à transitions. Here, we also consider the possible role of the C̃ ²B₁ (C̃) state in the wave packet evolution, because it exhibits a CI seam with X̃ at α = 180° and relatively large O–H distances (≈5 Bohr). Given the (large) internuclear distances at which this new CI is located, we have enlarged the grid of nuclear coordinates with respect to that of our previous calculation.

Additionally, we take advantage of the PESs and nonadiabatic data obtained for H₂O⁺ to carry out a study of isotopic effects in the fragmentation of D₂O⁺(B̃) and HDO⁺(B̃). Our theoretical fragmentation branching ratios of D₂O⁺ compare satisfactorily with the experimental results of Eland,¹¹ and they support the small differences found between H₂O⁺ and D₂O⁺ species in photoelectron–photoion-coincidence experiments.^13,14 With respect to HDO⁺, time of flight coincidence experiments^23,24 in collisions of H⁺ and F⁷⁺ with HDO have shown a strong preference for breaking the O–H bond after both single and double ionization of the molecular species. In those experiments, it was also found that the isotopic effect is independent of the projectile, suggesting a kinematic effect at play in the post-collisional breakdown of HDO⁺ and HDO²⁺. This interpretation was only supported by wave-packet simulations for HDO²⁺ in a single PES,²⁴ but one expects the fragmentation of HDO⁺ to be dominated by transitions between several electronic states. Here, we study the isotopic dependence in the bond cleavage preference, by carrying out 4-PESs wave packet propagation for HDO⁺.

In the (e–2e) + (e–e + ion) experiments of Tan et al.,¹² the ionization of H₂O is fast enough, compared to the characteristic rotational and vibrational periods of the molecule, that one can assume that the initial H₂O⁺ wave packet is obtained by a FC transition from the ground state of H₂O. In this work, we also address the photodissociation experiments of Harbo et al.,²⁵ where the fragmentation starts from specific excited vibrational levels in H₂O⁺(B̃). These authors carried out a crossed-beam experiment, where a H₂O⁺ beam was crossed by a 532 nm laser. The subsequent dissociation was interpreted as due to the vibrational excitation of metastable low-lying vibrational levels of B̃. The experiment measured the ratio between the cross sections for production of the fragments OH + H⁺ and OH⁺ + H, obtaining a ratio value Γ:

in contrast with the value Γ ≈ 0.31 reported by Tan et al.¹² The calculation of Γ with new initial conditions that simulate the experiment of Harbo et al. is a stringent test of both the electronic structure data and the dynamical methods employed.

The paper is organized as follows. In Section 2 we summarize the method employed to calculate the PESs and the dynamical couplings. The treatment of the two conical intersections found between B̃–Ã and X̃–C̃ PESs, which are responsible for the nonadiabatic transitions, is discussed in Section 3. In Section 4, we describe the numerical method applied to solve the time-dependent Schrödinger equation. Our main results are shown in Section 5, which includes the results of a 4-state calculation using a FC initial wave packet, the isotopic effect and the calculations with initial conditions that are excited vibrational states of the B̃ state. A summary is presented in Section 6.

2 Electronic structure calculations

We have performed ab initio calculations for the electronic states X̃ ²B₁, Ã ²A₁, B̃ ²B₂ and C̃ ²B₁ of H₂O⁺ using the MOLPRO package.²⁶ PESs and nonadiabatic couplings have been obtained with multireference configuration interaction calculations in the C_s symmetry point group. Configurations are generated by allocating seven valence electrons in (7a′,1a′′) molecular orbitals, previously obtained with CASSCF calculations in a space of (10a′,1a′′) orbitals. The basis set employed is aug-cc-pvqz and aug-cc-pvtz²⁷ for oxygen and hydrogen, respectively.

In this paper, we use both internal (r₁,r₂,α) and Jacobi (r,R,θ) coordinates to describe the molecular geometry, as depicted in Fig. 1, where r₁ and r₂ are the O–H distances, α is the bond angle, R is the distance from the center of mass of the OH fragment to the remaining hydrogen, r = r₁, and θ is the angle between r and R.

Fig. 1 Internal (left) and Jacobi (right) coordinates for the HOH nuclear framework.

In our calculations, the equilibrium geometry of the neutral H₂O, in the ground state, is located at r_e = r₁ = r₂ = 1.79 Bohr and α_e = 104.6° (θ_e = 107.8°). The equilibrium geometry of the ground state X̃ ²B₁ of H₂O⁺ is also of C_2v symmetry, with no significant changes in the equilibrium coordinates (r_e = 1.90 Bohr and θ_e = 112.5°). In contrast, the equilibrium geometry of à ²A₁ has D_∞h symmetry (α_e = 180°). The excited B̃ ²B₂ state is again of C_2v symmetry, but the r_e distance is 2.16 Bohr, larger than in the X̃ state, while the bond angle is θ_e = 80°, considerably smaller than in the X̃ state. These data, the energies of the corresponding potential energy minima and the dissociation energies, are given in Table 1.

Table 1 C _2v critical points of the first four electronic states of H₂O⁺. The energy reference is taken as the minimum energy of the H₂O⁺ (X̃ ²B₁)

State	r e (Bohr)	α e (deg.)	V min (eV)	V ∞ (eV)	Fragments
gm: global minimum; lm: local minimum; sp: saddle point.
X̃ ^2B1	1.90	110.0	0.0^gm	5.674	OH^+ + H
à ^2A1	1.87	180.0	0.952^gm	6.245	OH + H^+
B̃ ^2B2	2.07	76.3	4.503^gm	7.850	OH^+ + H
C̃ ^2B1	3.10	90.0	10.572^sp	6.245	OH + H^+
C̃ ^2B1	2.87	180.0	7.850^sp
C̃ ^2B1	3.42	30.0	7.388^lm

---

We show in Fig. 2 and 3 the contour plots of the PESs of the four adiabatic electronic states involved in the fragmentation dynamics. One can note the above-mentioned minima of the X̃ and à surfaces. The changes in the slopes of the contour lines of the à and B̃ surfaces at α ≈ 80° are due to the CI between these surfaces. At C_2v geometries, the C̃ surface shows two relatively shallow minima; the one at α = 180° is a saddle-point in the (r,R)-plane, as seen in the bottom-right panel of Fig. 3.

Fig. 2 Colour maps of the X̃, Ã, B̃ and C̃ PESs of C_2v geometries, as functions of the O–H distance r and the bond angle α. The numbers in the colour-box are energies in Hartree.

Fig. 3 Same as in Fig. 2, but for linear geometries and the PESs are plotted as functions of the Jacobi coordinates, r and R.

In Fig. 4 we show the electronic energies of states X̃, Ã, B̃ and C̃ as functions of the Jacobi coordinate R, for fixed r = 1.98 Bohr, and four values of θ. In the panels corresponding to θ = 80° and θ = 100°, the energy curve of state C̃ appears above those of states X̃ and Ã, and crosses state B̃ at energies with small weight in the initial FC wave packet, illustrated with Gaussian shade in the second panel. At θ = 170°, the states à and X̃ are almost degenerate up to r ≃ 5 Bohr, where the potential energy curves of states C̃ and X̃ display an avoided crossing, that turns into a CI for θ = 180°.

Fig. 4 Potential energy curves of the H₂O⁺ system as functions of R for r = 1.98 Bohr and the values of θ indicated in the panels. The electronic states mentioned in the paper are labelled in the left panel, and the CI between X̃ and C̃ is indicated for θ = 180°. In the panel for θ = 100°, the energy distribution of the Franck–Condon wave packet is depicted with a Gaussian shaded area and the horizontal dashed-lines correspond to the vibrational energies in the adiabatic B̃ state.

3 Regularization of the conical intersections

The dynamical couplings are introduced in the numerical treatment by implementing a diabatic potential matrix in the grid representation. In fact, the nonadiabatic transitions take place in the neighbourhood of CIs, where the dynamical couplings are singular and, accordingly, a regularization procedure must be performed to remove these singularities. Given that the CIs are located in different regions of the configuration space, this regularization can be carried out for the CIs separately. The regularization procedure, already applied in our previous paper²¹ to treat the CI between PESs à and B̃ (henceforth referred as CI₁), uses the transformation angle:and the unitary transformationto build the diabatic states χ^d = {χ^d_i,χ^d_j} from the adiabatic ones χ^a = {χ^a_i,χ^a_j}:

χ^d = U(γ_ij)χ^a (4)

The adiabatic potential energy surface matrixis transformed into the diabatic potential matrix:by

V^d = U(γ_ij)V^aU^†(γ_ij) (7)

For two-centre systems the adiabatic–diabatic transformation can be carried out to completely remove the radial couplings 〈χ^d_i|∂/∂R|χ^d_j〉. However, for any three (or more) centre system, the diabatic states do not exist. The scope of the regularization process is to find states with small couplings near the CIs, removing the singularities. In our calculation, this is accomplished by finding a transformation angle, γ_ij, that minimizes the nonadiabatic couplings between the diabatic states in the CI region. Explicitly, if Q is any particular nuclear coordinate, the residual coupling in the diabatic basis set is given by

As seen in ref. 21, the ÖB̃ CI takes place along a seam of C_2v geometries, and the regularization is more easily carried out by employing the symmetry coordinates {x,y,α}

built from the internal coordinates r₁ and r₂ (see Fig. 1). In these coordinates, our calculations indicate that a CI seam is located at y = 0 for angles (in radians):The transformation angle γ_ÃB̃ is written in terms of the internal coordinates, and, similar to previous works,^28,29 the energy difference ΔV_ÃB̃ = V^d_B̃B̃ − V^d_ÃÃ is taken proportional to the distance to the CI seam in the y = 0 plane. As α_CI1′ = dα_CI1/dx is small in the range of x where the interaction is sizable, the energy difference can be approximated by

ΔV^d_ÃB̃(x) = [α(x) − α_CI1(x)]cos(α_CI1′) (11)

The interaction term is proportional to the antisymmetric coordinate y, including a dependence on α that gives an extra flexibility to the transformation angle to minimize the residual couplings:

V^d_ÃB̃(α,y) = (a₀ + a₁α + a₂α²)y (12)

where the values of the a_i parameters are given in Table 2.

Table 2 Parametrization used for CI₁ (eqn (12)) and CI₂ (eqn (13), (17) and (18)). The constants are expressed in atomic units

i	a i	b i	c i	d i
0	2.21	−0.0692	2.82	−1.57
1	0.94	1.11	0.367	−0.245
2	−2.88	2.25	−0.693

---

The second CI considered (CI₂, henceforth) takes place between PESs X̃ and C̃ at angles α ≈ π.

The symbols in Fig. 5 indicate the coordinates of the geometries with zero-energy difference between states X̃ and C̃. The position of the seam in the (x,y) plane is parametrized using the following expression:

in which the values of f_i are obtained by a non-linear least-square fitting procedure of the previous equation to the points of minimum energy difference for eight values of α. The f_i(α) functions are fitted, independently, to the following expressions:

f₀(α) = b₀α + c₀ (14)

f₁(α) = b₁ e^{−c₁(α−π)2} (15)

f₂(α) = b₂ e^{−c₂(α−π)2} (16)

Again, b_i and c_i parameters are obtained by a 1D-non-linear fitting procedure, and are given in Table 2.

Fig. 5 Bullets: position in coordinates (x, y), eqn (9), of the zero-energy differences between states X̃ and C̃ at α = π rad. Solid-line: CI₂ seam fitted-model of eqn (13).

The modelled diabatic potential matrix elements are finally written as

ΔV^d_X̃C̃(α,x,y) = x − x_CI2(α,y) (17)

V^d_X̃C̃(α,y) = d₀ e^d₁y2(α − π) (18)

where d₀ and d₁ are also free parameters obtained by non-linear fitting of this model to the ab initio couplings 〈ψ_X̃|∂/∂x|ψ_C̃〉, for several values of y, and whose values are given in Table 2.

The X̃–C̃ seam described by eqn (13) is plotted in Fig. 5, together with the position of the ab initio zero-energy difference between these states for linear geometries. Also, in Fig. 6, we compare, for α = 175°, the ab initio nonadiabatic couplings (symbols) with the values of the derivatives (solid lines) of the transformation angle γ_X̃C̃ [see eqn (2)] obtained with the fitted model (17) and (18), for a few values of y (those corresponding to the dotted lines of Fig. 5). These figures show the satisfactory fitting of the model to the ab initio data, where the peaks of the lorentzian follow the location of the CI₂ seam and their widths depend on the parameters of eqn (18). A more global view of the topology of the CI₂ is shown in Fig. 7, where we plot the derivative of the transformation angle γ_X̃C̃ with respect to the symmetry coordinates {x,y,α}. These derivatives mimic the behaviour of the nonadiabatic couplings between the states X̃ and C̃ near their CI.

Fig. 6 Comparison between ab initio 〈Ψ_X̃|∂/∂x|Ψ_C̃〉 (symbols) and the model of eqn (13) ∂γ_X̃C̃/∂x (solid-lines) as functions of x, for various values of y, labelled in the graph, and also indicated with dotted lines in Fig. 5. Modelled and ab initio data were obtained with α = 175°.

Fig. 7 Colour-map plot of ∂γ_X̃C̃/∂x (left), ∂γ_X̃C̃/∂y (middle) and ∂γ_X̃C̃/∂α (right) as functions of the symmetric (x) and antisymmetric (y) coordinates for α = 175°. Colour scale values in atomic units.

The treatment of ÖX̃ transitions, explained in our previous paper, is based on the method of Haxton et al.^30,31 that yields a couple of diabatic states. Assuming a small total angular momentum, the Coriolis coupling terms are neglected and the body-fixed nuclear Hamiltonian has the form, in Jacobi coordinates:^31,32

where j_z is the z component of the angular momentum operator for the rotation of the diatom OH in the body-fixed frame, with R along the ẑ-axis. In this approximation, the nuclear wave functions are assumed to be eigenfunctions of the body-fixed operator j_z. The operator j_z² is substituted³³ into eqn (19) by its eigenvalue, k² (k = K − M_L), which is only defined in the C_∞v or D_∞h symmetries, where K and M_L are the quantum numbers for the z projections of the total and electronic angular momenta, respectively. In the present case, the electronic states à ²A₁ and X̃ ²B₁ correlate, in the limit θ → 180°, with the two components of the ²Π_u state, both with Λ = |M_L| = 1. In this limit, the χ_±1 eigenfunctions of L_z with M_L = ±1 are linear combinations of the electronic wave functions, χ_à and χ_X̃:These two combinations are employed for all values of θ. Considering the transformations due to the regularization of the two CIs and the transformation (20), the unitary matrix, U, for the transformation from the adiabatic state vector χ^a = {χ^a_C̃,χ^a_X̃,χ^a_Ã,χ^a_B̃} to the diabatic one χ^d = {χ^d_C̃,χ₊,χ₋,χ^d_B̃} iswith γ₁ ≡ γ_X̃C̃ and γ₂ ≡ γ_ÃB̃. U relates the adiabatic, V^a, and diabatic, V^d, potentials:

V^d = UV^aU^† (22)

4 Computational dynamics

The dynamical calculation has been carried out by means of nuclear wave packet propagation on the diabatic electronic PESs, employing a computational method (GridTDSE)²² that is adept for massive parallel calculations. This technique consists of a direct numerical collocation integration of the 3D time-dependent Schrödinger equation (TDSE):where H is the Hamiltonian operator of eqn (19) with the approximation of substituting j_z² by k², as discussed in Section 3. The numerical calculation uses a lattice scheme,³⁴ with the nuclear wave function evaluated at the points of a 3D grid. To incorporate nonadiabatic transitions, the wave function, Ψ, is written as a column vector with components {Ψ_i}, where i labels the (diabatic) electronic states included. The nuclear dynamics on each PES are coupled through the interaction matrix elements H^d_ij(R), following a procedure similar to that of Ma et al.³⁵ This leads to the matrix equation

H^dΨ = (T + V^d)Ψ = i[capital Psi, Greek, dot above] (24)

where V^d stores the values of the potentials at the grid points. T, the kinetic energy operator, is a non-diagonal sparse matrix that is obtained by applying the finite differences method (see ref. 34), considering a stencil of 15 points in each dimension.

The calculations were carried out with a grid of 83 points in R ∈ [1,11] Bohr, 80 points in r ∈ [1.3,11] Bohr and 73 points in θ ∈ [ε,π–ε] rad, with ε = 2 × 10⁻⁴. In such a numerical integration scheme, the singularity at θ = 0 can be eluded by evaluating Ψ on the surroundings of the troublesome points.

The wave function vector components {Ψ_i}, each one referring to a different PES, were simultaneously propagated in time using a second order difference scheme (SOD), specially adapted to the present case since it conserves norm and energy, accumulating errors only in the phase.³⁶ The errors are minimized in our calculations by setting the timestep Δt = 0.18 attosecond (7.5 × 10⁻³ a.u.), ensuring the stability of the norm and the correct account of the coupling terms between the evolving wave functions on each PES.

To avoid non-physical reflections of the wave function at the grid walls (L_max = 11 Bohr), we have included a damping function³⁷ in the propagation scheme. This function leads to a smooth decay of the wave packet for q > L_max − δ:

with q ≡ (R,r), δ = 1 Bohr and ξ = 0.05 Bohr⁻².

Once the description of the nuclear wave packet in each electronic state is achieved, the integrated probability density represents the time-dependent population of the electronic state i. Such populations can change either due to nonadiabatic transitions among the states included in the model, or due to the dissociating fluxes. The corresponding fragmentation probabilities ρ_i(t) are obtained from the continuity equation in its integral form: , evaluating the outgoing flux at R, r = 10.0 Bohr for the adiabatic components of the wave packet on each electronic state. Finally, the comparison with the experimental measurements is carried out at long times t > 5 ps, where the fragmentation probabilities reach asymptotic stable values.

5 Results and discussion

5.1 Franck–Condon initial wave packet

The initial nuclear wave function is built taking the lowest vibrational eigenfunction of the electronic ground state of the water molecule. This PES was obtained with the same procedure and atomic basis set as the electronic states of H₂O⁺ described in Section 2. In order to solve the corresponding time-independent Schrödinger equation (TISE), we employed the Lanczos scheme³⁸ implemented in the GridTDSE code, which is based on the description of Guo et al.³⁹ The analysis of this wave function, when projected over the eigenstates of the adiabatic B̃ PES, shows that it can be described by a wave packet that spans a range of vibrational energies of about 0.1 Hartree, as shown by the vertical Gaussian curve in Fig. 4 for θ = 100°. There, the horizontal dashed-lines correspond to the vibrational energies obtained in the adiabatic B̃ PES.

5.2 Effect of increasing the grid-size

As seen in Section 3, CI₂ between the X̃ and C̃ states appears at (R,r) distances between 3 and 5 Bohr. The correct description of the nonadiabatic transitions in the vicinity of this CI requires increasing the grid size from (R_max,r_max) ≈ 7 Bohr (GRID I²¹) to (R_max,r_max) ≈ 11 Bohr (GRID II). To check the effect of size-enlargement of the grid on the time-evolution of the electronic state populations, we have reproduced our previous three-state calculation in the new grid. The comparison of the new results with those of ref. 21 are presented in Fig. 8(a). There are small differences in the probabilities between both calculations, with small changes in the probabilities, the X̃ state being the more affected.

Fig. 8 Population of the electronic states of H₂O⁺ (labelled in the panels) as a function of time. (a) The wave packet propagation includes only three electronic states. Dashed-orange lines: grid I with L_max = 7 Bohr; solid lines: grid II with L_max = 11 Bohr [see eqn (25)]. (b) The wave packet propagation includes four electronic states (dashed lines) in grid II with L_max = 11 Bohr. Solid-orange lines are the solid lines of (a).

To study this in more detail, we show in Fig. 9 the sum of the populations of states X̃, Ã and B̃ for both box-sizes calculations (lines without symbols). As the size of the grid increases, there is a displacement of the population fall towards longer times, due to the damping function acting later. The inset in that figure shows the difference between both calculations, which would correspond to the wave function probability in the extra region. According to this figure, the wave packet enters this region after 50 fs, where it monotonously grows during 800 fs. After that, the amount of wave packet that leaves the extra region exceeds the amount that enters the region from shorter distances and the wave packet in the extra region fades away.

Fig. 9 Sum of the populations of the states X̃, Ã and B̃ as functions of time for two boxes with L_max = 7 Bohr (– – –) and L_max = 11 Bohr (—) compared to the probability of finding a particle described by the 1D-wave packet of eqn (26), inside a box with L = 7 Bohr (– • –) and L = 11 Bohr (—•—). Inset: Differences between the populations in the boxes with L = 11 Bohr and L = 7 Bohr for the actual calculation (—) and the 1-D model (—•—).

This behaviour can be readily modelled by the free expansion of a 1D wave packet. In this model, we consider the generic wave packet

ψ(x,t) = exp[−ζ(t)(x − λ₁)²]sin(λ₂x) (26)

where the time-dependent exponential parameter ζ(t) is introduced to simulate the widening of the wave packet with energy above the dissociation threshold, and

ζ = (λ₃ + λ₄t^λ₅)⁻¹ (27)

With this wave packet, the probability of finding the system in the 1D box of length L can be obtained analytically and is plotted in Fig. 9 for the same grid size as in the three-state calculation, and with parameters (in atomic units) λ₁ = 2.5, λ₂ = 0.385, λ₃ = 0.10, λ₄ = 2.8 × 10⁻³ and λ₅ = 1.2, obtained to approximately fit the time-dependent probability of the three states of H₂O⁺. These results suggest that the expansion of the H₂O⁺ wave packet is similar to the free wave packet expansion for the region of the space of coordinates R,r > 7 Bohr, as considered in our previous work.

5.3 Four-state calculation

The initial wave function shown in Section 5.1 was also propagated in the extended grid II with a Hamiltonian matrix that included four electronic states {X̃, Ã, B̃, C̃}. Here, the state C̃ is only coupled to X̃ through nonadiabatic couplings described by the transformation angle γ_X̃C̃.

The time evolution of the populations of the four electronic states is displayed in Fig. 8(b), for a specific projection of the total angular momentum, K = 1. For an easy comparison, the populations obtained in the three-state calculation, in the same grid, are also included. These results show that state C̃ gets populated at the cost of X̃ at short times, but in the picosecond timescale, the populations of the three- and four-state calculations converge.

As it was shown in ref. 21, the case K = 1 is a representative of the fragmentation probabilities at low temperatures. Here, the H₂O⁺(B̃) main fragmentation mechanism remains the same as in the previous three-state calculation: a two-step process in which there is a fast population transfer from B̃ to à in the surroundings of CI₁ and, after 30 fs, a population redistribution and decay through the following reactions:

H₂O⁺(Ã) → H⁺ + OH k₁ (28)

H₂O⁺(Ã) ⇌ H₂O⁺(X̃) k_±2 (29)

while the state X̃ dissociates through

H₂O⁺(X̃) → H + OH⁺(³Σ⁻) k₃ (30)

where k_i are the corresponding reaction constants. States X̃ and C̃ are connected by a CI, and we have

H₂O⁺(X̃) ⇌ H₂O⁺(C̃) k_±4 (31)

while C̃ dissociates through

H₂O⁺(C̃) → H⁺ + OH k₅ (32)

The direct fragmentation from B̃ is

H₂O⁺(B̃) → H + OH⁺(¹Σ⁺) k₆ (33)

The probability for fragmentation along each channel is shown in Fig. 10 for the wave packet propagation up to 1600 fs, and these results are used to fit a model for a first-order dissociation mechanism:

p_i(t) = r_i(1 − e^−k_it) (34)

The parameters obtained by fitting the ab initio data from t = 1000 fs to t = 1600 fs are used to draw the thin lines in Fig. 10, and are given in Table 3. There, the fraction of H₂O⁺ that remains stable is r₀ = 1 − r₁ − r₃ − r₅ − r₆.

Fig. 10 Fragment probabilities as functions of time in the four-state calculation. Thick-dashed lines: result of the wave packet propagation. Thin-solid lines: fitted model [eqn (34)] using data within the shaded area. The number in brackets indicates the spin multiplicity of the diatomic species. The fitting parameters are given in Table 3.

Table 3 Branching ratios and rate constants for the first-order dissociation dynamics model (34) applied to reactions (28)–(33) of H₂O⁺(B̃) and D₂O⁺(B̃) fragmentation

i	H2O^+		D2O^+
	r i	k i (fs^−1)	r i	k i (fs^−1)
0	0.091		0.074
1	0.219	2.4 × 10^−3	0.138	1.3 × 10^−3
3	0.658	5.1 × 10^−4	0.786	6.7 × 10^−4
5	0.026	1.2 × 10^−3	—	—
6	0.006	3.7 × 10^−3	0.002	1.3 × 10^−3

---

In Fig. 8(b), one can note the small populations of state C̃. The differences found in the populations of states à and X̃ by the introduction of state C̃ vanish as t increases, which would eventually lead to similar fragmentation branching ratios. However, one must take into account that we have extended the grid of our previous paper²¹ to incorporate transitions in the X̃–C̃ CI. When the grid size is increased, the system remains within the grid limits for a longer time, which, in our calculation, leads to larger populations of states à and X̃ with respect to those obtained in our previous three-state work (see Fig. 8(a)).

Ferreira et al.¹⁶ have suggested that transitions B̃–C̃ near their PESs crossing might explain their Gaussian structure around 1.8 eV in the H⁺ kinetic energy release distribution. However, states B̃ and C̃ are of different symmetry (A′ and A′′, respectively) and, correspondingly, their interaction is zero in the Hamiltonian matrix. Also, we note that their leading electronic configurations are 1a²₁2a²₁1b₂3a²₁1b²₁ (B̃) and 1a²₁2a²₁1b²₂3a₁1b₁4a₁ (C̃), which are not linked by one-electron operators, and thus we expect a small B̃–C̃ coupling. Additionally, the fraction of the initial FC wave packet near the crossings (see Fig. 4) is very small. Therefore, we think that this mechanism is unlikely to be competitive with the formation of OH + H⁺via B̃–Ã transitions.

5.4 Isotopic effect

Having established in the previous section the small role of state C̃ in the fragmentation ratio, we address the calculation of the branching ratio of a FC wave packet in the B̃ state of D₂O⁺ and HDO⁺ isotopomers including only three electronic states (X̃, à and B̃). For D₂O⁺, the results are compared with those of H₂O⁺ in Fig. 11, also including the experimental data of Eland¹¹ for D₂O⁺ and Tan et al.¹² for H₂O⁺, in excellent agreement with our calculation. The ratio OD⁺/D⁺ (≈5.7), larger than that of OH⁺/H⁺ (≈3.0) (see the values of r₃/r₁ in Table 3), can be understood in terms of kinematic effects: the D₂O⁺ wave packet is slower and less spread than the H₂O⁺ one, and this is particularly important in the region of linear geometries when the RT coupling transfers part of the wave packet from à to X̃. In D₂O⁺, the RT coupling is less effective due to the larger reduced mass of the fragments (see eqn (19)), but the slower motion increases the interaction time, which results in a faster increase of the population of X̃ at the expense of Ã, and the consequent increase in the production of OD⁺ and decrease of that of D⁺.

Fig. 11 Comparison between fragmentation probabilities of H₂O⁺ and D₂O⁺. Circles are the experimental data of Eland¹¹ for D₂O⁺ and triangles those of Tan et al.¹² for H₂O⁺.

Turning to the HDO⁺ system, it can fragment along four channels: OH⁺ + D, OD⁺ + H, OH + D⁺ and OD + H⁺. Here, we start with the FC wave packet as in H₂O⁺, and follow its dynamics in the 4 PESs. The fragmentation probabilities are presented in Fig. 12, where they are compared with the coincidence-time-of-flight measurements,²⁴ and the dynamic parameters of the first order model (34) are given in Table 4. The calculated σ(OD⁺)/σ(OH⁺) and σ(H⁺)/σ(D⁺) asymptotic ratios are 2.0 and 1.6, respectively, in excellent agreement with the experimental²⁴ values of 2.1 ± 0.3 and 1.5 ± 0.2.

Fig. 12 Upper panel: Fragmentation probabilities of HDO⁺(B̃) as functions of time. The numbers in brackets indicate the spin multiplicity of the molecular fragment species. Lower panel: Branching ratios for bond cleavage preference compared with the experimental²⁴ ones.

Table 4 Branching ratios and rate constants for the first-order dissociation dynamics model (34) applied to the fragmentation of HDO⁺(B̃)

Fragments	r	k (fs^−1)
OD^+(^3Σ^−) + H	0.47	6.1 × 10^−4
OH^+(^3Σ^−) + D	0.23	6.2 × 10^−4
OD + H^+	0.11	1.4 × 10^−3
OH + D^+	0.069	1.3 × 10^−3
OD^+(^1Σ^−) + H	0.0022	1.3 × 10^−3
OH^+(^1Σ^−) + D	0.0019	1.2 × 10^−3

---

In the calculation, OH⁺ and OD⁺ are formed through the predissociation of state B̃, which is linked to the single ionization from the 1b₂ molecular orbital. However, as shown by coincidence experiments,¹² the dissociation into H⁺ + OD or D⁺ + OH can also take place by one-electron removal from the innermost valence orbital 2a₁. In fact, the experimental branching ratios of Tan et al.¹² show that the cross section for formation of protons in the photoionization of H₂O is given by

σ(H⁺) = 0.22σ(1b₂) + 0.74σ(2a₁) (35)

where σ(1b₂) and σ(2a₁) are the cross sections for ionization from the molecular orbitals 1b₂ and 2a₁, respectively. The experiments²⁴ were carried out at projectile energies of the order of 1 MeV u⁻¹, where the cross sections for ionization from both shells can be estimated by applying the Bethe–Born expression, which leads to similar contributions of the orbitals 1b₂ and 2a₁ to the formation of H⁺ (see ref. 40). Therefore, the good agreement between our calculation, which only includes the electron removal from 1b₂, and the experiment, supports the explanation that the ratio σ(H⁺)/σ(D⁺) is mainly a kinematic effect,²⁴ which is common to the wave packet propagation in all electronic PESs.

5.5 Time evolution of an excited vibrational state of B̃

In contrast with the experiment of Tan et al.,¹² where the fragmentation of H₂O⁺(B̃) ions is recorded just after their formation, in the experiment of Harbo et al.,²⁵ the H₂O⁺ ions produced by the ionization of H₂O are stored for a relatively long time; most of the ions undergo fragmentation and only those with an energy below the dissociation threshold survive and are, afterwards, irradiated with a 532 nm (2.33 eV) laser, and the ensuing fragmentation branching ratios measured. The B̃ state only holds two vibrational states below the lowest asymptotic energy of state X̃ and the observed fragmentation was explained as due to the photodissociation of these metastable vibrational states of B̃.

To simulate the experimental²⁵ conditions, we have time-propagated an initial wave packet that corresponds to the vibrational eigenfunction whose energy lies 2.33 eV above the vibrational ground state in the adiabatic B̃ surface. This vibrational state, whose probability density is shown in Fig. 14, is the seven symmetric stretch overtone. The low-lying vibrational states of the B̃ electronic state have long lifetimes (≈198 μs),²⁵ which points to ineffective nonadiabatic transitions from those states. Our simulation starts when the interaction of H₂O⁺ ions with the laser field ends. Then, the wavefunction of the excited vibrational state has sizeable values around the ÖB̃ CI and fast nonadiabatic transitions take place, leading to the ion breakdown.

Fig. 13 shows the evolution of this wave packet for K = 1, the quantum number for the z-component of the total angular momentum [see eqn (19)]. The comparison of the populations of Fig. 8 with those of Fig. 13(a) clearly shows similar mechanisms, but now there is a larger fragmentation into H⁺ + OH, which results in a faster decrease of the population of Ã. This precludes the transition to X̃, henceforth a smaller fragmentation into H + OH⁺; i.e., the initial vibrational excitation leads to a faster motion of the wave packet in the à PES and to the ensuing increase of the k₁ rate constant of reaction (28).

Fig. 13 Time evolution of the nuclear wave packet that is initially equal to the 8th vibrational wave function (given in Fig. 14) of the adiabatic B̃ state. (a) Populations of the electronic states of H₂O⁺ indicated in the figure. (b) Probabilities of the possible fragments.

Fig. 14 Wave function probability density (shade) of the seven symmetric stretch overtone in B̃ adiabatic PES (contour dotted lines), for C_2v symmetries, as functions of the symmetric coordinate x = (r₁ + r₂)/2 and the internal angle α.

Starting with the excited vibrational state mentioned above, our theoretical fragmentation branching ratio is Γ = σ(H⁺)/σ(OH⁺) = 0.96, in satisfactory agreement with the experimental value²⁵ of 1.3. Furthermore, our value increases to 1.04, which falls within the ±0.3 experimental error bars, when considering only the triplet OH⁺ species. In both cases, these ratios are considerably larger than the 0.3 ± 0.01 value obtained in the experiments of Tan et al.¹² and our theoretical²¹ value of 0.33 with the FC initial wave packet.

6 Conclusions

We have carried out the study of the nonadiabatic fragmentation of H₂O⁺(B̃) by means of wave packet propagation in four electronic PESs (X̃, Ã, B̃ and C̃), using a modified version of the GridTDSE computational code²² to simultaneously reproduce the nuclear dynamics on each surface. The values of the potential energy on each electronic state and the nonadiabatic couplings among them were calculated ab initio on exactly the same grid of points where the finite difference scheme of the GridTDSE is implemented, enhancing the accuracy of the calculations. This study extends our previous simulations²¹ by considering one more electronic state (C̃) and a larger configuration space. Here, we have considered not only an initial FC wave packet, obtained by vertical ionization of H₂O⁺ into the B̃ predissociative state of H₂O⁺, but also vibrationally excited H₂O⁺ ions and the fragmentation of water isotopomers.

The new branching ratio for the fragmentation of H₂O⁺(B̃) into OH + H⁺ and OH⁺ + H agrees with previous (X̃, Ã, B̃) three-state calculations,²¹ which supports the three-state mechanism. We have also obtained an isotopic dependence of the branching ratio in agreement with the experimental^11,24 data and we have discussed the unlikely role of the B̃–C̃ transitions in the Gaussian peak of the H⁺ kinetic energy release measured by Ferreira et al.¹⁶

We have carried out 4-PESs propagation of the initial FC wave packet in the B̃ state of HDO⁺. The good agreement between the results of the simulation and the ion–water collision experiments²⁴ suggests that the cleavage preference found is due to kinematic effects in the evolution of the wave packet in the surfaces that lead to dissociation (à and X̃). Also, we have obtained good agreement with the experiments¹¹ in the fragmentation branching ratio of D₂O⁺(B̃), and explained the differences with respect to the fragmentation of H₂O⁺ as due to the larger transitions between X̃ and à states near the linear geometry because of the slower and more compact wave packet in heavy water than in H₂O.

We have also considered the fragmentation of vibrationally excited H₂O⁺(B̃) ions to simulate the experimental conditions of Harbo et al.²⁵ We have found that the initial excitation in the symmetric stretch drives the dynamics of the initial wave packet in a way that favors the dissociation from state Ã, leading to branching ratios quite different from those obtained starting with an initial FC-wave packet. In this case, our results agree with the experiments.²⁵

Finally, the good agreement between our calculations and the experiments paves the way for the application of this wave packet propagation technique to other three-center systems, see e.g.ref. 41.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

This work has been partially supported by Ministerio de Economía and Competitividad (Spain), project ENE2014-52432-R. The Centro de Computación Científica of UAM is acknowledged for the computational hosting facilities.

References

1. D. A. Neufeld, J. R. Goicoechea, P. Sonnentrucker, J. H. Black, J. Pearson, S. Yu, T. G. Phillips, D. C. Lis, M. De Luca, E. Herbst, P. Rimmer, M. Gerin, T. A. Bell, F. Boulanger, J. Cernicharo, A. Coutens, E. Dartois, M. Kazmierczak, P. Encrenaz, E. Falgarone, T. R. Geballe, T. Giesen, B. Godard, P. F. Goldsmith, C. Gry, H. Gupta, P. Hennebelle, P. Hily-Blant, C. Joblin, R. Kołos, J. Krełowski, J. Martín-Pintado, K. M. Menten, R. Monje, B. Mookerjea, M. Perault, C. Persson, R. Plume, M. Salez, S. Schlemmer, M. Schmidt, J. Stutzki, D. Teyssier, C. Vastel, A. Cros, K. Klein, A. Lorenzani, S. Philipp, L. A. Samoska, R. Shipman, A. G. G. M. Tielens, R. Szczerba and J. Zmuidzinas, Astron. Astrophys., 2010, 521, L10.
2. E. González-Alfonso, J. Fischer, S. Bruderer, H. S. P. Müller, J. Graciá-Carpio, E. Sturm, D. Lutz, A. Poglitsch, H. Feuchtgruber, S. Veilleux, A. Contursi, A. Sternberg, S. Hailey-Dunsheath, A. Verma, N. Christopher, R. Davies, R. Genzel and L. Tacconi, Astron. Astrophys., 2013, 550, A25.
3. R. Wayne, Chemistry of Atmospheres: An Introduction to the Chemistry of the Atmospheres of Earth, the Planets, and Their Satellites, Oxford University Press, 2000.
4. M. A. Disanti, U. Fink and A. B. Schultz, Icarus, 1990, 86, 152–171.
5. S. A. Fuselier, K. Altwegg, H. Balsiger, J. J. Berthelier, A. Bieler, C. Briois, T. W. Broiles, J. L. Burch, U. Calmonte, G. Cessateur, M. Combi, J. De Keyser, B. Fiethe, M. Galand, S. Gasc, T. I. Gombosi, H. Gunell, K. C. Hansen, M. Hässig, A. Jäckel, A. Korth, L. Le Roy, U. Mall, K. E. Mandt, S. M. Petrinec, S. Raghuram, H. Rème, M. Rinaldi, M. Rubin, T. Sémon, K. J. Trattner, C.-Y. Tzou, E. Vigren, J. H. Waite and P. Wurz, Astron. Astrophys., 2015, 583, A2.
6. C. R. Brundle and D. W. Turner, Proc. R. Soc. London, Ser. A, 1968, 307, 27–36.
7. J. E. Reutt, L. S. Wang, Y. T. Lee and D. A. Shirley, J. Chem. Phys., 1986, 85, 6928–6939.
8. S. Truong, A. Yencha, A. Juarez, S. Cavanagh, P. Bolognesi and G. King, Chem. Phys., 2009, 355, 183–193.
9. A. Lorquet and J. Lorquet, Chem. Phys., 1974, 4, 353–367.
10. M. Eroms, M. Jungen and H.-D. Meyer, J. Phys. Chem. A, 2010, 114, 9893–9901.
11. J. Eland, Chem. Phys., 1975, 11, 41–47.
12. K. Tan, C. Brion, P. V. der Leeuw and M. van der Wiel, Chem. Phys., 1978, 29, 299–309.
13. I. Powis and D. J. Reynolds, J. Chem. Soc., Faraday Trans., 1991, 87, 921–926.
14. K. Norwood, A. Ali and C. Y. Ng, J. Chem. Phys., 1991, 95, 8029–8037.
15. A. L. F. de Barros, J. Lecointre, H. Luna, M. B. Shah and E. C. Montenegro, Phys. Rev. A, 2009, 80, 012716.
16. N. Ferreira, L. Sigaud and E. C. Montenegro, Phys. Rev. A, 2017, 96, 012705.
17. N. Ferreira, L. Sigaud and E. C. Montenegro, J. Phys. Chem. A, 2017, 121, 3234–3238.
18. M. Paniagua, R. Martínez, P. Gamallo and M. González, Phys. Chem. Chem. Phys., 2014, 16, 23594–23603.
19. R. Martínez, M. Paniagua, J. Mayneris-Perxachs, P. Gamallo and M. González, Phys. Chem. Chem. Phys., 2017, 19, 3857–3868.
20. P. Gamallo, P. Defazio, M. González, M. Paniagua and C. Petrongolo, Phys. Chem. Chem. Phys., 2015, 17, 23392–23402.
21. J. Suárez, L. Méndez and I. Rabadán, J. Phys. Chem. Lett., 2015, 6, 72–76.
22. J. Suarez, S. Farantos, S. Stamatiadis and L. Lathouwers, Comput. Phys. Commun., 2009, 180, 2025–2033.
23. I. Ben-Itzhak, A. M. Sayler, M. Leonard, J. Maseberg, D. Hathiramani, E. Wells, M. Smith, J. Xia, P. Wang, K. Carnes and B. Esry, Nucl. Instrum. Methods Phys. Res., Sect. B, 2005, 233, 284–292.
24. A. M. Sayler, M. Leonard, K. D. Carnes, R. Cabrera-Trujillo, B. D. Esry and I. Ben-Itzhak, J. Phys. B: At., Mol. Opt. Phys., 2006, 39, 1701.
25. L. S. Harbo, S. Dziarzhytski, C. Domesle, G. Brenner, A. Wolf and H. B. Pedersen, Phys. Rev. A, 2014, 89, 052520.
26. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schütz, WIREs Comput. Mol. Sci., 2012, 2, 242–253.
27. T. H. Dunning Jr., J. Chem. Phys., 1989, 90, 1007–1023.
28. D. Elizaga, L. F. Errea, A. Macías, L. Méndez, A. Riera and A. Rojas, J. Phys. B: At., Mol. Opt. Phys., 1999, 32, L697.
29. P. Barragán, L. F. Errea, A. Macías, L. Méndez, I. Rabadán, A. Riera, J. M. Lucas and A. Aguilar, J. Chem. Phys., 2004, 121, 11629–11638.
30. D. J. Haxton, T. N. Rescigno and C. W. McCurdy, Phys. Rev. A, 2007, 75, 012711.
31. D. J. Haxton, T. N. Rescigno and C. W. McCurdy, Phys. Rev. A, 2007, 76, 049907(E).
32. C. Leforestier, J. Chem. Phys., 1991, 94, 6388–6397.
33. S. Carter and N. Handy, Mol. Phys., 1984, 52, 1367–1391.
34. R. Guantes and S. C. Farantos, J. Chem. Phys., 1999, 111, 10827–10835.
35. J. Ma, C. Xie, X. Zhu, D. R. Yarkony, D. Xie and H. Guo, J. Phys. Chem. A, 2014, 118, 11926–11934.
36. C. Leforestier, R. Bisseling, C. Cerjan, M. Feit, R. Friesner, A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, H.-D. Meyer, N. Lipkin, O. Roncero and R. Kosloff, J. Comput. Phys., 1991, 94, 59–80.
37. D. Dundas, J. F. McCann, J. S. Parker and K. T. Taylor, J. Phys. B: At., Mol. Opt. Phys., 2000, 33, 3261.
38. C. Lanczos, J. Res. Natl. Bur. Stand., Sect. B, 1950, 45, 255–282.
39. H. Guo, R. Chen and D. Xie, J. Theor. Comput. Chem., 2002, 01, 173–185.
40. C. Illescas, L. F. Errea, L. Méndez, B. Pons, I. Rabadán and A. Riera, Phys. Rev. A, 2011, 83, 052704.
41. A. D. Smith, E. M. Warne, D. Bellshaw, D. A. Horke, M. Tudorovskya, E. Springate, A. J. H. Jones, C. Cacho, R. T. Chapman, A. Kirrander and R. S. Minns, Phys. Rev. Lett., 2018, 120, 183003.

---

Footnotes

† Electronic supplementary information (ESI) available: A file containing the four adiabatic potential energy surfaces employed in the calculation is provided. See DOI: 10.1039/c8cp03725f

‡ Present address: Dipartimento di Chimica, Università degli Studi di Milano, 20133 Milano, Italy.

---