Nonadiabatic fragmentation of H_{2}O^{+} and isotopomers. Wave packet propagation using ab initio wavefunctions†
Received
12th June 2018
, Accepted 16th September 2018
First published on 12th November 2018
1 Introduction
The water cation is naturally found in several environments such as the interstellar medium,^{1,2} and atmospheres of planets^{3} and comets.^{4,5} It can be formed in collisions of electrons, ions or photons with H_{2}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_{2}O with photon energies below 21.2 eV. At these energies, the photoelectron spectrum shows three bands that correspond to the formation of H_{2}O^{+} in the electronic states ^{2}B_{1} (hereafter denoted ), Ã ^{2}A_{1} (Ã) and ^{2}B_{2} (). These electronic states are essentially described by oneelectron removal from the three outermost molecular orbitals (1b_{1}, 3a_{1} and 1b_{2}, respectively) of H_{2}O. For photon energies higher than the dissociation threshold (18.1 eV), the vibrational levels of the state are predissociative and, accordingly, the third band of the H_{2}O photoelectron spectrum exhibits broad peaks. The ionization of H_{2}O above the dissociation threshold leads to the cation fragmentation into H^{+} + OH, H + OH^{+} and O^{+} + H + H. Following the ionization of H_{2}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_{2}O^{+} system include the ab initio studies^{18,19} of the reactions O + H_{2}^{+} → OH^{+} + H and O + H_{2}^{+} → OH + H^{+}, and the corresponding reactions with the isotopomers D_{2}^{+} and HD^{+}. These calculations used the quasiclassicaltrajectory and realwavepacket methods, and stateoftheart potential energy surfaces of the and Ã electronic states. The influence of the Renner–Teller coupling on the O + H_{2}^{+} proton transfer reactions was analyzed in ref. 20.
In a previous paper,^{21} we carried out a theoretical study of the fragmentation of the water cation from the 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.^{12} Our calculations consisted in nuclear wavepacket propagation on three electronic surfaces (, Ã and ), using a modified version of the GridTDSE code^{22} 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_{2}O. The wavepacket propagation showed that the dissociation mechanism involves nonadiabatic transitions from the state to Ã and . The transition –Ã 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 , moves towards the CI and is transferred to Ã with an important excitation of the bending mode. At linear geometries, the states Ã and become the two components of the ^{2}Π_{u} state, and the Ã– 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 , leading to dissociation into OH^{+} + H.
Although the good agreement between our results and the experimental branching ratios indicates that the threestate (–Ã–) mechanism explains the predissociation of the water cation, other possible nonadiabatic transitions must be ruled out. In our previous work,^{21} we checked that the intersystem crossing transitions from and Ã to the dissociative ã ^{4}B_{1} state are very slow compared to –Ã transitions. Here, we also consider the possible role of the ^{2}B_{1} () state in the wave packet evolution, because it exhibits a CI seam with 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_{2}O^{+} to carry out a study of isotopic effects in the fragmentation of D_{2}O^{+}() and HDO^{+}(). Our theoretical fragmentation branching ratios of D_{2}O^{+} compare satisfactorily with the experimental results of Eland,^{11} and they support the small differences found between H_{2}O^{+} and D_{2}O^{+} species in photoelectron–photoioncoincidence experiments.^{13,14} With respect to HDO^{+}, time of flight coincidence experiments^{23,24} in collisions of H^{+} and F^{7+} 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 postcollisional breakdown of HDO^{+} and HDO^{2+}. This interpretation was only supported by wavepacket simulations for HDO^{2+} in a single PES,^{24} 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 4PESs wave packet propagation for HDO^{+}.
In the (e–2e) + (e–e + ion) experiments of Tan et al.,^{12} the ionization of H_{2}O is fast enough, compared to the characteristic rotational and vibrational periods of the molecule, that one can assume that the initial H_{2}O^{+} wave packet is obtained by a FC transition from the ground state of H_{2}O. In this work, we also address the photodissociation experiments of Harbo et al.,^{25} where the fragmentation starts from specific excited vibrational levels in H_{2}O^{+}(). These authors carried out a crossedbeam experiment, where a H_{2}O^{+} beam was crossed by a 532 nm laser. The subsequent dissociation was interpreted as due to the vibrational excitation of metastable lowlying vibrational levels of . The experiment measured the ratio between the cross sections for production of the fragments OH + H^{+} and OH^{+} + H, obtaining a ratio value Γ:

 (1) 
in contrast with the value
Γ ≈ 0.31 reported by Tan
et al.^{12} 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 –Ã and – 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 timedependent Schrödinger equation. Our main results are shown in Section 5, which includes the results of a 4state calculation using a FC initial wave packet, the isotopic effect and the calculations with initial conditions that are excited vibrational states of the state. A summary is presented in Section 6.
2 Electronic structure calculations
We have performed ab initio calculations for the electronic states ^{2}B_{1}, Ã ^{2}A_{1}, ^{2}B_{2} and ^{2}B_{1} of H_{2}O^{+} using the MOLPRO package.^{26} 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 augccpvqz and augccpvtz^{27} for oxygen and hydrogen, respectively.
In this paper, we use both internal (r_{1},r_{2},α) and Jacobi (r,R,θ) coordinates to describe the molecular geometry, as depicted in Fig. 1, where r_{1} and r_{2} 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_{1}, 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_{2}O, in the ground state, is located at r_{e} = r_{1} = r_{2} = 1.79 Bohr and α_{e} = 104.6° (θ_{e} = 107.8°). The equilibrium geometry of the ground state ^{2}B_{1} of H_{2}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 Ã ^{2}A_{1} has D_{∞h} symmetry (α_{e} = 180°). The excited ^{2}B_{2} state is again of C_{2v} symmetry, but the r_{e} distance is 2.16 Bohr, larger than in the state, while the bond angle is θ_{e} = 80°, considerably smaller than in the 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_{2}O^{+}. The energy reference is taken as the minimum energy of the H_{2}O^{+} ( ^{2}B_{1})
State 
r
_{e} (Bohr) 
α
_{e} (deg.) 
V
_{min} (eV) 
V
_{∞} (eV) 
Fragments 
gm: global minimum; lm: local minimum; sp: saddle point. 
^{2}B_{1} 
1.90 
110.0 
0.0^{gm} 
5.674 
OH^{+} + H 
Ã ^{2}A_{1} 
1.87 
180.0 
0.952^{gm} 
6.245 
OH + H^{+} 
^{2}B_{2} 
2.07 
76.3 
4.503^{gm} 
7.850 
OH^{+} + H 
^{2}B_{1} 
3.10 
90.0 
10.572^{sp} 
6.245 
OH + H^{+} 
^{2}B_{1} 
2.87 
180.0 
7.850^{sp} 


^{2}B_{1} 
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 abovementioned minima of the and Ã surfaces. The changes in the slopes of the contour lines of the Ã and surfaces at α ≈ 80° are due to the CI between these surfaces. At C_{2v} geometries, the surface shows two relatively shallow minima; the one at α = 180° is a saddlepoint in the (r,R)plane, as seen in the bottomright panel of Fig. 3.

 Fig. 2 Colour maps of the , Ã, and PESs of C_{2v} geometries, as functions of the O–H distance r and the bond angle α. The numbers in the colourbox 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 , Ã, and 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 appears above those of states and Ã, and crosses state at energies with small weight in the initial FC wave packet, illustrated with Gaussian shade in the second panel. At θ = 170°, the states Ã and are almost degenerate up to r ≃ 5 Bohr, where the potential energy curves of states and display an avoided crossing, that turns into a CI for θ = 180°.

 Fig. 4 Potential energy curves of the H_{2}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 and 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 dashedlines correspond to the vibrational energies in the adiabatic 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^{21} to treat the CI between PESs Ã and (henceforth referred as CI_{1}), uses the transformation angle: 
 (2) 
and the unitary transformation 
 (3) 
to 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 matrix 
 (5) 
is transformed into the diabatic potential matrix: 
 (6) 
by 
V^{d} = U(γ_{ij})V^{a}U^{†}(γ_{ij})  (7) 
For twocentre 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

 (8) 
As seen in ref. 21, the Ã– 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,α}

 (9) 
built from the internal coordinates
r_{1} and
r_{2} (see
Fig. 1). In these coordinates, our calculations indicate that a CI seam is located at
y = 0 for angles (in radians):

 (10) 
The transformation angle
γ_{Ã} is written in terms of the internal coordinates, and, similar to previous works,
^{28,29} the energy difference Δ
V_{Ã} =
V^{d}_{} −
V^{d}_{ÃÃ} is taken proportional to the distance to the CI seam in the
y = 0 plane. As
α_{CI1}′ = d
α_{CI1}/d
x is small in the range of
x where the interaction is sizable, the energy difference can be approximated by

ΔV^{d}_{Ã}(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}_{Ã}(α,y) = (a_{0} + a_{1}α + a_{2}α^{2})y  (12) 
where the values of the
a_{i} parameters are given in
Table 2.
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_{2}, henceforth) takes place between PESs and at angles α ≈ π.
The symbols in Fig. 5 indicate the coordinates of the geometries with zeroenergy difference between states and . The position of the seam in the (x,y) plane is parametrized using the following expression:

 (13) 
in which the values of
f_{i} are obtained by a nonlinear leastsquare 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_{0}(α) = b_{0}α + c_{0}  (14) 

f_{1}(α) = b_{1}e^{−c1(α−π)2}  (15) 

f_{2}(α) = b_{2}e^{−c2(α−π)2}  (16) 
Again,
b_{i} and
c_{i} parameters are obtained by a 1Dnonlinear fitting procedure, and are given in
Table 2.

 Fig. 5 Bullets: position in coordinates (x, y), eqn (9), of the zeroenergy differences between states and at α = π rad. Solidline: CI_{2} seam fittedmodel of eqn (13).  
The modelled diabatic potential matrix elements are finally written as

ΔV^{d}_{}(α,x,y) = x − x_{CI2}(α,y)  (17) 

V^{d}_{}(α,y) = d_{0}e^{d1y2}(α − π)  (18) 
where
d_{0} and
d_{1} are also free parameters obtained by nonlinear fitting of this model to the
ab initio couplings 〈
ψ_{}∂/∂
x
ψ_{}〉, for several values of
y, and whose values are given in
Table 2.
The – seam described by eqn (13) is plotted in Fig. 5, together with the position of the ab initio zeroenergy 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 γ_{} [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_{2} seam and their widths depend on the parameters of eqn (18). A more global view of the topology of the CI_{2} is shown in Fig. 7, where we plot the derivative of the transformation angle γ_{} with respect to the symmetry coordinates {x,y,α}. These derivatives mimic the behaviour of the nonadiabatic couplings between the states and near their CI.

 Fig. 6 Comparison between ab initio 〈Ψ_{}∂/∂xΨ_{}〉 (symbols) and the model of eqn (13) ∂γ_{}/∂x (solidlines) 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 Colourmap plot of ∂γ_{}/∂x (left), ∂γ_{}/∂y (middle) and ∂γ_{}/∂α (right) as functions of the symmetric (x) and antisymmetric (y) coordinates for α = 175°. Colour scale values in atomic units.  
The treatment of Ã– 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 bodyfixed nuclear Hamiltonian has the form, in Jacobi coordinates:^{31,32}

 (19) 
where
j_{z} is the
z component of the angular momentum operator for the rotation of the diatom OH in the bodyfixed frame, with
R along the
ẑaxis. In this approximation, the nuclear wave functions are assumed to be eigenfunctions of the bodyfixed operator
j_{z}. The operator
j_{z}^{2} is substituted
^{33} into
eqn (19) by its eigenvalue,
k^{2} (
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 Ã
^{2}A
_{1} and
^{2}B
_{1} correlate, in the limit
θ → 180°, with the two components of the
^{2}Π
_{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
χ_{}:

 (20) 
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}_{},
χ^{a}_{},
χ^{a}_{Ã},
χ^{a}_{}} to the diabatic one
χ^{d} = {
χ^{d}_{},
χ_{+},
χ_{−},
χ^{d}_{}} is

 (21) 
with
γ_{1} ≡
γ_{} and
γ_{2} ≡
γ_{Ã}.
U relates the adiabatic,
V^{a}, and diabatic,
V^{d}, potentials:
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)^{22} that is adept for massive parallel calculations. This technique consists of a direct numerical collocation integration of the 3D timedependent Schrödinger equation (TDSE): 
 (23) 
where H is the Hamiltonian operator of eqn (19) with the approximation of substituting j_{z}^{2} by k^{2}, as discussed in Section 3. The numerical calculation uses a lattice scheme,^{34} 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.^{35} This leads to the matrix equation 
H^{d}Ψ = (T + V^{d})Ψ = i  (24) 
where V^{d} stores the values of the potentials at the grid points. T, the kinetic energy operator, is a nondiagonal 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^{−4}. 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.^{36} The errors are minimized in our calculations by setting the timestep Δt = 0.18 attosecond (7.5 × 10^{−3} 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 nonphysical reflections of the wave function at the grid walls (L_{max} = 11 Bohr), we have included a damping function^{37} in the propagation scheme. This function leads to a smooth decay of the wave packet for q > L_{max} − δ:

 (25) 
with
q ≡ (
R,
r),
δ = 1 Bohr and
ξ = 0.05 Bohr
^{−2}.
Once the description of the nuclear wave packet in each electronic state is achieved, the integrated probability density represents the timedependent 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_{2}O^{+} described in Section 2. In order to solve the corresponding timeindependent Schrödinger equation (TISE), we employed the Lanczos scheme^{38} implemented in the GridTDSE code, which is based on the description of Guo et al.^{39} The analysis of this wave function, when projected over the eigenstates of the adiabatic 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 dashedlines correspond to the vibrational energies obtained in the adiabatic PES.
5.2 Effect of increasing the gridsize
As seen in Section 3, CI_{2} between the and 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^{21}) to (R_{max},r_{max}) ≈ 11 Bohr (GRID II). To check the effect of sizeenlargement of the grid on the timeevolution of the electronic state populations, we have reproduced our previous threestate 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 state being the more affected.

 Fig. 8 Population of the electronic states of H_{2}O^{+} (labelled in the panels) as a function of time. (a) The wave packet propagation includes only three electronic states. Dashedorange 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. Solidorange 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 , Ã and for both boxsizes 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 , Ã and 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 1Dwave 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 1D 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 − λ_{1})^{2}]sin(λ_{2}x)  (26) 
where the timedependent exponential parameter
ζ(
t) is introduced to simulate the widening of the wave packet with energy above the dissociation threshold, and

ζ = (λ_{3} + λ_{4}t^{λ5})^{−1}  (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 threestate calculation, and with parameters (in atomic units)
λ_{1} = 2.5,
λ_{2} = 0.385,
λ_{3} = 0.10,
λ_{4} = 2.8 × 10
^{−3} and
λ_{5} = 1.2, obtained to approximately fit the timedependent probability of the three states of H
_{2}O
^{+}. These results suggest that the expansion of the H
_{2}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 Fourstate 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 {, Ã, , }. Here, the state is only coupled to through nonadiabatic couplings described by the transformation angle γ_{}.
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 threestate calculation, in the same grid, are also included. These results show that state gets populated at the cost of at short times, but in the picosecond timescale, the populations of the three and fourstate 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_{2}O^{+}() main fragmentation mechanism remains the same as in the previous threestate calculation: a twostep process in which there is a fast population transfer from to Ã in the surroundings of CI_{1} and, after 30 fs, a population redistribution and decay through the following reactions:

H_{2}O^{+}(Ã) → H^{+} + OH k_{1}  (28) 

H_{2}O^{+}(Ã) ⇌ H_{2}O^{+}() k_{±2}  (29) 
while the state
dissociates through

H_{2}O^{+}() → H + OH^{+}(^{3}Σ^{−}) k_{3}  (30) 
where
k_{i} are the corresponding reaction constants. States
and
are connected by a CI, and we have

H_{2}O^{+}() ⇌ H_{2}O^{+}() k_{±4}  (31) 
while
dissociates through

H_{2}O^{+}() → H^{+} + OH k_{5}  (32) 
The direct fragmentation from
is

H_{2}O^{+}() → H + OH^{+}(^{1}Σ^{+}) k_{6}  (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 firstorder dissociation mechanism:

p_{i}(t) = r_{i}(1 − e^{−kit})  (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
_{2}O
^{+} that remains stable is
r_{0} = 1 −
r_{1} −
r_{3} −
r_{5} −
r_{6}.

 Fig. 10 Fragment probabilities as functions of time in the fourstate calculation. Thickdashed lines: result of the wave packet propagation. Thinsolid 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 firstorder dissociation dynamics model (34) applied to reactions (28)–(33) of H_{2}O^{+}() and D_{2}O^{+}() fragmentation
i

H_{2}O^{+} 
D_{2}O^{+} 
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 . The differences found in the populations of states Ã and by the introduction of state 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^{21} to incorporate transitions in the – 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 with respect to those obtained in our previous threestate work (see Fig. 8(a)).
Ferreira et al.^{16} have suggested that transitions – near their PESs crossing might explain their Gaussian structure around 1.8 eV in the H^{+} kinetic energy release distribution. However, states and 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^{2}_{1}2a^{2}_{1}1b_{2}3a^{2}_{1}1b^{2}_{1} () and 1a^{2}_{1}2a^{2}_{1}1b^{2}_{2}3a_{1}1b_{1}4a_{1} (), which are not linked by oneelectron operators, and thus we expect a small – 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 –Ã transitions.
5.4 Isotopic effect
Having established in the previous section the small role of state in the fragmentation ratio, we address the calculation of the branching ratio of a FC wave packet in the state of D_{2}O^{+} and HDO^{+} isotopomers including only three electronic states (, Ã and ). For D_{2}O^{+}, the results are compared with those of H_{2}O^{+} in Fig. 11, also including the experimental data of Eland^{11} for D_{2}O^{+} and Tan et al.^{12} for H_{2}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_{3}/r_{1} in Table 3), can be understood in terms of kinematic effects: the D_{2}O^{+} wave packet is slower and less spread than the H_{2}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 . In D_{2}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 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_{2}O^{+} and D_{2}O^{+}. Circles are the experimental data of Eland^{11} for D_{2}O^{+} and triangles those of Tan et al.^{12} for H_{2}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_{2}O^{+}, and follow its dynamics in the 4 PESs. The fragmentation probabilities are presented in Fig. 12, where they are compared with the coincidencetimeofflight measurements,^{24} 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^{24} values of 2.1 ± 0.3 and 1.5 ± 0.2.

 Fig. 12 Upper panel: Fragmentation probabilities of HDO^{+}() 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^{24} ones.  
Table 4 Branching ratios and rate constants for the firstorder dissociation dynamics model (34) applied to the fragmentation of HDO^{+}()
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 , which is linked to the single ionization from the 1b_{2} molecular orbital. However, as shown by coincidence experiments,^{12} the dissociation into H^{+} + OD or D^{+} + OH can also take place by oneelectron removal from the innermost valence orbital 2a_{1}. In fact, the experimental branching ratios of Tan et al.^{12} show that the cross section for formation of protons in the photoionization of H_{2}O is given by

σ(H^{+}) = 0.22σ(1b_{2}) + 0.74σ(2a_{1})  (35) 
where
σ(1b
_{2}) and
σ(2a
_{1}) are the cross sections for ionization from the molecular orbitals 1b
_{2} and 2a
_{1}, respectively. The experiments
^{24} were carried out at projectile energies of the order of 1 MeV u
^{−1}, 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
_{2} and 2a
_{1} to the formation of H
^{+} (see
ref. 40). Therefore, the good agreement between our calculation, which only includes the electron removal from 1b
_{2}, and the experiment, supports the explanation that the ratio
σ(H
^{+})/
σ(D
^{+}) is mainly a kinematic effect,
^{24} which is common to the wave packet propagation in all electronic PESs.
5.5 Time evolution of an excited vibrational state of
In contrast with the experiment of Tan et al.,^{12} where the fragmentation of H_{2}O^{+}() ions is recorded just after their formation, in the experiment of Harbo et al.,^{25} the H_{2}O^{+} ions produced by the ionization of H_{2}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 state only holds two vibrational states below the lowest asymptotic energy of state and the observed fragmentation was explained as due to the photodissociation of these metastable vibrational states of .
To simulate the experimental^{25} conditions, we have timepropagated an initial wave packet that corresponds to the vibrational eigenfunction whose energy lies 2.33 eV above the vibrational ground state in the adiabatic surface. This vibrational state, whose probability density is shown in Fig. 14, is the seven symmetric stretch overtone. The lowlying vibrational states of the electronic state have long lifetimes (≈198 μs),^{25} which points to ineffective nonadiabatic transitions from those states. Our simulation starts when the interaction of H_{2}O^{+} ions with the laser field ends. Then, the wavefunction of the excited vibrational state has sizeable values around the Ã– 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 zcomponent 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 , 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_{1} 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 state. (a) Populations of the electronic states of H_{2}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 adiabatic PES (contour dotted lines), for C_{2v} symmetries, as functions of the symmetric coordinate x = (r_{1} + r_{2})/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^{25} 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.^{12} and our theoretical^{21} value of 0.33 with the FC initial wave packet.
6 Conclusions
We have carried out the study of the nonadiabatic fragmentation of H_{2}O^{+}() by means of wave packet propagation in four electronic PESs (, Ã, and ), using a modified version of the GridTDSE computational code^{22} 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^{21} by considering one more electronic state () and a larger configuration space. Here, we have considered not only an initial FC wave packet, obtained by vertical ionization of H_{2}O^{+} into the predissociative state of H_{2}O^{+}, but also vibrationally excited H_{2}O^{+} ions and the fragmentation of water isotopomers.
The new branching ratio for the fragmentation of H_{2}O^{+}() into OH + H^{+} and OH^{+} + H agrees with previous (, Ã, ) threestate calculations,^{21} which supports the threestate 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 – transitions in the Gaussian peak of the H^{+} kinetic energy release measured by Ferreira et al.^{16}
We have carried out 4PESs propagation of the initial FC wave packet in the state of HDO^{+}. The good agreement between the results of the simulation and the ion–water collision experiments^{24} 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 ). Also, we have obtained good agreement with the experiments^{11} in the fragmentation branching ratio of D_{2}O^{+}(), and explained the differences with respect to the fragmentation of H_{2}O^{+} as due to the larger transitions between and Ã states near the linear geometry because of the slower and more compact wave packet in heavy water than in H_{2}O.
We have also considered the fragmentation of vibrationally excited H_{2}O^{+}() ions to simulate the experimental conditions of Harbo et al.^{25} 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 FCwave packet. In this case, our results agree with the experiments.^{25}
Finally, the good agreement between our calculations and the experiments paves the way for the application of this wave packet propagation technique to other threecenter 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 ENE201452432R. The Centro de Computación Científica of UAM is acknowledged for the computational hosting facilities.
References
 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. HilyBlant, C. Joblin, R. Kołos, J. Krełowski, J. MartínPintado, 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 CrossRef.
 E. GonzálezAlfonso, 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. HaileyDunsheath, A. Verma, N. Christopher, R. Davies, R. Genzel and L. Tacconi, Astron. Astrophys., 2013, 550, A25 CrossRef.

R. Wayne, Chemistry of Atmospheres: An Introduction to the Chemistry of the Atmospheres of Earth, the Planets, and Their Satellites, Oxford University Press, 2000 Search PubMed.
 M. A. Disanti, U. Fink and A. B. Schultz, Icarus, 1990, 86, 152–171 CrossRef CAS.
 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 CrossRef.
 C. R. Brundle and D. W. Turner, Proc. R. Soc. London, Ser. A, 1968, 307, 27–36 CrossRef CAS.
 J. E. Reutt, L. S. Wang, Y. T. Lee and D. A. Shirley, J. Chem. Phys., 1986, 85, 6928–6939 CrossRef CAS.
 S. Truong, A. Yencha, A. Juarez, S. Cavanagh, P. Bolognesi and G. King, Chem. Phys., 2009, 355, 183–193 CrossRef CAS.
 A. Lorquet and J. Lorquet, Chem. Phys., 1974, 4, 353–367 CrossRef CAS.
 M. Eroms, M. Jungen and H.D. Meyer, J. Phys. Chem. A, 2010, 114, 9893–9901 CrossRef CAS PubMed.
 J. Eland, Chem. Phys., 1975, 11, 41–47 CrossRef CAS.
 K. Tan, C. Brion, P. V. der Leeuw and M. van der Wiel, Chem. Phys., 1978, 29, 299–309 CrossRef CAS.
 I. Powis and D. J. Reynolds, J. Chem. Soc., Faraday Trans., 1991, 87, 921–926 RSC.
 K. Norwood, A. Ali and C. Y. Ng, J. Chem. Phys., 1991, 95, 8029–8037 CrossRef CAS.
 A. L. F. de Barros, J. Lecointre, H. Luna, M. B. Shah and E. C. Montenegro, Phys. Rev. A, 2009, 80, 012716 CrossRef.
 N. Ferreira, L. Sigaud and E. C. Montenegro, Phys. Rev. A, 2017, 96, 012705 CrossRef.
 N. Ferreira, L. Sigaud and E. C. Montenegro, J. Phys. Chem. A, 2017, 121, 3234–3238 CrossRef CAS PubMed.
 M. Paniagua, R. Martínez, P. Gamallo and M. González, Phys. Chem. Chem. Phys., 2014, 16, 23594–23603 RSC.
 R. Martínez, M. Paniagua, J. MaynerisPerxachs, P. Gamallo and M. González, Phys. Chem. Chem. Phys., 2017, 19, 3857–3868 RSC.
 P. Gamallo, P. Defazio, M. González, M. Paniagua and C. Petrongolo, Phys. Chem. Chem. Phys., 2015, 17, 23392–23402 RSC.
 J. Suárez, L. Méndez and I. Rabadán, J. Phys. Chem. Lett., 2015, 6, 72–76 CrossRef PubMed.
 J. Suarez, S. Farantos, S. Stamatiadis and L. Lathouwers, Comput. Phys. Commun., 2009, 180, 2025–2033 CrossRef CAS.
 I. BenItzhak, 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 CrossRef CAS.
 A. M. Sayler, M. Leonard, K. D. Carnes, R. CabreraTrujillo, B. D. Esry and I. BenItzhak, J. Phys. B: At., Mol. Opt. Phys., 2006, 39, 1701 CrossRef CAS.
 L. S. Harbo, S. Dziarzhytski, C. Domesle, G. Brenner, A. Wolf and H. B. Pedersen, Phys. Rev. A, 2014, 89, 052520 CrossRef.
 H.J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schütz, WIREs Comput. Mol. Sci., 2012, 2, 242–253 CrossRef CAS.
 T. H. Dunning Jr., J. Chem. Phys., 1989, 90, 1007–1023 CrossRef.
 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 CrossRef CAS.
 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 CrossRef PubMed.
 D. J. Haxton, T. N. Rescigno and C. W. McCurdy, Phys. Rev. A, 2007, 75, 012711 CrossRef.
 D. J. Haxton, T. N. Rescigno and C. W. McCurdy, Phys. Rev. A, 2007, 76, 049907(E) CrossRef.
 C. Leforestier, J. Chem. Phys., 1991, 94, 6388–6397 CrossRef CAS.
 S. Carter and N. Handy, Mol. Phys., 1984, 52, 1367–1391 CrossRef CAS.
 R. Guantes and S. C. Farantos, J. Chem. Phys., 1999, 111, 10827–10835 CrossRef CAS.
 J. Ma, C. Xie, X. Zhu, D. R. Yarkony, D. Xie and H. Guo, J. Phys. Chem. A, 2014, 118, 11926–11934 CrossRef CAS PubMed.
 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 CrossRef.
 D. Dundas, J. F. McCann, J. S. Parker and K. T. Taylor, J. Phys. B: At., Mol. Opt. Phys., 2000, 33, 3261 CrossRef CAS.
 C. Lanczos, J. Res. Natl. Bur. Stand., Sect. B, 1950, 45, 255–282 CrossRef.
 H. Guo, R. Chen and D. Xie, J. Theor. Comput. Chem., 2002, 01, 173–185 CrossRef CAS.
 C. Illescas, L. F. Errea, L. Méndez, B. Pons, I. Rabadán and A. Riera, Phys. Rev. A, 2011, 83, 052704 CrossRef.
 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 CrossRef PubMed.
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. 

This journal is © the Owner Societies 2018 