Correlated proton-electron hole dynamics in protonated water clusters upon extreme ultraviolet photoionization

The ultrafast nuclear and electronic dynamics of protonated water clusters H+(H2O)n after extreme ultraviolet photoionization is investigated. In particular, we focus on cluster cations with n = 3, 6, and 21. Upon ionization, two positive charges are present in the cluster related to the excess proton and the missing electron, respectively. A correlation is found between the cluster's geometrical conformation and initial electronic energy with the size of the final fragments produced. For situations in which the electron hole and proton are initially spatially close, the two entities become correlated and separate in a time-scale of 20 to 40 fs driven by strong non-adiabatic effects.


I. INTRODUCTION
Protonated water clusters of general formula H þ (H 2 O) n constitute an important model of the solvated proton in liquid water. 1 The fact that clusters of the desired size can be prepared and controlled in the gas phase has made possible multitude of spectroscopic and fragmentation dynamics studies on clusters featuring the excess proton motif. [2][3][4][5] Coincidence measurement experiments without time information and based on the detection of fragmentation products on small clusters H 3 O þ and H þ (H 2 O) 2 have revealed the probabilities for different fragmentation channels upon valence photoionization as well as kinetic energy releases of the fragments. 6,7 Time-resolved coincidence techniques to follow the molecular dynamics upon photoionization have been recently applied to systems of similar size. 8,9 Although no time-resolved XUV photoionization experiments have been performed for these systems up to date, full quantum dynamics 10 and mixed quantum-classical molecular dynamics calculations 11 revealed extremely fast correlations between the excess proton and the electron hole for the HðH 2 OÞ 2þ 2 cation related to strong non-Born-Oppenheimer effects.
Chemically, proton transfer processes triggered by the sudden change of electronic structure are known to protect DNA base pairs from photo-damage by UV-visible sunlight absorption. 12,13 Ionized water can be produced in the last stages of decay cascades upon collisions of energetic particles, such as high-energy photons, with biological systems. Not much is yet known about the detailed fragmentation channels and the time-resolved dynamics of ionized protonated water clusters, specially how the positive charges related to an excess proton and an electron hole separate and in what time scale and in what conditions finally the electron hole results in neutral radical fragments while a new solvated proton is released into the system. In this work, we apply mixed quantum-classical surface hopping calculations [14][15][16][17][18][19] to study the ultrafast dynamics and fragmentation products of photoionized H þ (H 2 O) 3 , H þ (H 2 O) 6 , and H þ (H 2 O) 21 cations. The non-adiabatic dynamics are based on ab initio electronic structure potentials and non-adiabatic couplings. We identify several elementary channels characterized by the water cation H 2 O þ and the hydroxyl radical OH • , which are the two main final products that can accommodate the electron vacancy, and hydronium related fragments, e.g., H 3 O þ and H þ (H 2 O) 2 , which carry the excess proton. For small clusters, we find a correlation between the initial amount of electronic energy available in the system, which is directly connected to the binding energy of the photoelectron, and the size of the protonated fragments produced, which becomes smaller for more energetic systems. The shape of the cluster and the initial water molecule which is ionized determine as well whether the electron hole can escape the system in the form of a H 2 O þ cation or whether the charge and spin finally separate yielding OH • or other neutral fragments.
We also focus our attention on the protonated H þ (H 2 O) 21 20,21 with a centrally solvated excess proton inside an almost spherical shell and which can be viewed as a finite model of the solvated proton in liquid water. Removal of an electron from the central water molecule holding the excess proton leads to a fast separation of the two charges and to a delocalization of the hole in the external shell of the cluster, whereas the hydronium remains stable in the center.
The article is organized as follows. Section II presents the methodological and computational details. In Sec. III A, the electronic properties and fragmentation dynamics of the protonated water clusters H þ (H 2 O) 3 and H þ (H 2 O) 6 are discussed. Section III B discusses the dynamics upon photoionization of the HðH 2 OÞ 2þ 21 cation and Sec. IV concludes. In the Appendix, we briefly describe the blockwise coupling scheme for the fewest-switching surface hopping method that was employed in this work and its convergence properties are shown on a simple model Hamiltonian.

II. COMPUTATIONAL DETAILS
The coupled nuclear-electronic non-Born-Oppenheimer dynamics of the protonated water clusters is described within the fewest switches surface hopping (FSSH) scheme. 14,15 A set of phase space initial conditions was sampled from a Boltzmann distribution at 300 K on the ground electronic state potential energy surface of the non-ionized system. In the FSSH scheme, each nuclear trajectory evolves independently and no quantum correlations among them are considered. The array of probability amplitudes of the electronic eigenstates are obtained with the usual expression for the time-dependent electronic Schr€ odinger equation for a single trajectory where d a JI ¼ h/ J ðr; RÞjr R j/ I ðr; RÞi is the non-Born-Oppenheimer coupling vector and the a index labels individual trajectories. The electronic Hamiltonian matrix H JI ¼ d JI E J is diagonal in the adiabatic representation, and E J is the electronic potential energy of J-th electronic state. The FSSH method introduces electronic state switching events for each trajectory with probability g IJ of hopping from the I-th state to the J-th state within the time interval [t, t þ dt], This switching probability guarantees that the number of switching events required to keep the number of trajectories on each electronic state proportional to the corresponding quantum populations q a J ðtÞ ¼ C Ãa J ðtÞC a J ðtÞ is minimal. 14 The forces on the nuclei are calculated from the potential energy surface gradient on the single electronic state that a given trajectory occupies at a given time.
In order to more efficiently treat systems with a large number of electronic states, we compute the d JI coupling vectors only within a block of states around the electronic state populated at a given time for a particular trajectory. Effectively, this assumes that electronic states not in energetic proximity are decoupled, which is justified in view of the equivalence For H þ (H 2 O) 21 , the block of states for which the non-Born-Oppenheimer is coupling matrix has radius of 5 states, maximally 11 states are taken into consideration. In the Appendix, we apply this block scheme to a model Hamiltonian and show that it quickly converges with the size of the block of states considered. The initial conditions for the three cation models were sampled from a Boltzmann distribution in phase space at T ¼ 300 K on the ground electronic state of the non-ionized system. Along each FSSH trajectory, the gradients and non-adiabatic couplings were calculated ab initio on-the-fly using the complete active space self-consistent field (CASSCF) method implemented in the quantum chemistry package MOLCAS. [22][23][24] Specifically, CASSCF(9,7)/6-31G(d) and CASSCF (11,9) 21 , the electronic states of the singly ionized molecule were obtained from a one-hole (1 h) approximation of the Hamiltonian of the ionized system. 25 This takes advantage of the fact that E J ðRÞ h/ J ðr; RÞjH e ðr; RÞj/ I ðr; RÞi (4) which is Koopmans theorem. 26 Here, j/ J ðr; RÞi ¼ c † J j/ HF ðr; RÞi, where c † J annihilates an electron from the J-th spin-orbital. Therefore, all excited state electronic wavefunctions of the photoionized system j/ J ðr; RÞi are obtained from a single Hartree-Fock calculation of the nonionized system j/ HF ðr; RÞi.

III. RESULTS AND DISCUSSION
A. Water clusters H 1 (H 2 O) 3 and H 1 (H 2 O) 6

Ultrafast electronic energy relaxation dynamics
Initially, we focus our attention on the Coulomb explosion dynamics upon photoionization of the H þ (H 2 O) 3 and two different conformers of the H þ (H 2 O) 6 cluster starting from different electronic states of the clusters characterized by different ionization potentials. The H þ (H 2 O) 3 cation can be seen as an extension of the Zundel cation by adding an extra water molecule. However, its symmetric stabilization of the excess proton in the form a hydronium in the central water monomer is reminiscent of the Eigen cation H þ (H 2 O) n (n ¼ 4), with one of the external water molecules is missing. The H þ (H 2 O) 6 cation has the possibility to form spatially closed structures by hydrogen bonding. The cage structure is considered to be the most stable conformation of the H þ (H 2 O) 6 cation with an optimal network of hydrogen bonds. 27 For the H þ (H 2 O) 6 , the extended Zundel arrangement, in which a Zundel core is surrounded by four other water molecules that act as hydrogen bond acceptors, is also considered. Although this is not the most stable gas phase structure, it represents an important transient structure in the stabilization of the excess proton in solution. 1 Fig. 1 illustrates stable geometries of these protonated water clusters, H þ (H 2 O) n (n ¼ 3, 6), in their ground electronic states before ionization. The orbitals shown correspond to the ones missing an electron (the photoelectron) in the leading electronic configuration of the corresponding electronic states of the ionized cluster. Fig. 2 shows the adiabatic potential energy curves of the HðH 2 OÞ 2þ 3 dication, as the proton moves between the oxygen atom of the central hydronium cation (z ¼ 0.5) and the oxygen atom of the side water monomer (z ¼ À0.5). 28 The topology of the adiabatic potential energy curves, with multiple true or avoided crossings among them, is indicative of the strong non-adiabatic coupling between the proton coordinates along hydrogen bonds and the electronic subsystem, which will characterize the ultrafast electronic relaxation dynamics in those systems. The scaled coordinate z of proton is defined as where R OO ¼ 2.45 Å is the equilibrium distance between the two oxygen atoms, D 0 ¼ 0.95 Å is defined as the distance between the central proton and each of the oxygen atoms when z is found at the edge of its grid, and the origin is set at the middle of two oxygen atoms. These In general, the electronic energy relaxation dynamics taking place during the Coulomb explosion of these small protonated water clusters takes place within $100 fs or even shorter time scales. This is illustrated by the temporal evolution of electronic state populations presented in Fig. 3. A very fast population decay from the initially populated excited electronic state takes place already after $10 fs, and the ground electronic state becomes the most highly populated state within $100 fs, which signals the end of the electronic relaxation process.
After photoionization, the Coulomb repulsion between the excess proton and the hole quickly leads to a separation of these two entities. In small clusters, in which fragmentation inevitably occurs within short time scales, two major fragmentation channels are present: one major channel containing the excess proton and electronically closed shell and another major channel of radical nature containing the electron vacancy

043203-5 Z. Li and O. Vendrell
where n ¼ m þ l. How these major channels give raise to different fragmentation products is found to be strongly dependent on the cluster geometry and on the initial level of electronic excitation present in the photoionized cluster.
As an illustrative example of fragmentation into these major channels, we look at a particular trajectory for the H þ (H 2 O) 6 case initiated from a highly excited electronic state with an energy $8.7 eV higher than the ground electronic state of the ion state. Fig. 4 represents the motion of the electron hole, which starts localized in the region of the central Zundel motif and in the same region occupied by the excess proton. As can be seen, 75 fs after photoionization, the cluster has split through the two central water molecules originally connected by the excess proton. The excess proton is found after fragmentation in a H þ (H 2 O) 3 moiety, while the electron hole is initially found in a ðH 2 OÞ þ 3 radical cation and finally settles down as an OH • fragment plus a H þ (H 2 O) 2 protonated cluster. Accompanying the electronic relaxation and fragmentation, the rapidly growing distance between proton and electron hole due to Coulomb repulsion is shown in Fig. 5.
In the following, we focus on the different fragmentation channels of the excess proton and hole moieties of the photoionized clusters as a function of cluster size and shape and of their initial electronic energy. explosion of the three cationic clusters considered. These fractions are defined as the ratio of trajectories with this fragment to the total number of trajectories, and a fragment is counted when the distance between its oxygen atoms and all other oxygen atoms is larger than 2.9 Å at the end of a trajectory at t ¼ 120 fs. The fragment fractions are given in Table I.
For the HðH 2 OÞ 2þ 3 cation produced in the D 2 state, in which the electron hole is delocalized over the cluster structure, a relatively uniform yield of H 3 O þ and H þ (H 2 O) 2 fragments from the parent ion is found, with a slightly favored H þ (H 2 O) 2 channel. However, when ionized into the D 0 electronic state, in which the electron hole is localized on one of the external water molecules, the diversity of fragmentation channels disappears and only the H þ (H 2 O) 2 fragment is produced. As discussed later, both for the D 0 and D 2 initial states, the electron vacancy localizes on one of the external water molecules. The excess proton is then carried by the other two water molecules in first instance. However, the D 2 state has a higher electronic energy that is sufficient to fragment the H þ (H 2 O) 2

cation into H 3 O þ and H 2 O.
For the extended Zundel HðH 2 OÞ 2þ 6 , we find that the system always initially separates by the central excess proton bridge into two groups of three water molecules, one holding the excess proton and the other holding the hole. This is the case for an initially localized hole, as in the D 0 state, as well as for a delocalized hole as in D 2 . A higher electronic energy favors the production of the smaller H þ (H 2 O) 2 fragment over H þ (H 2 O) 3 . A very similar conclusion is reached for the HðH 2 OÞ 2þ 6 system in cage conformation. As seen in Table I, the fraction of small H 3 O þ fragments increases with the level of initial electronic excitation, whereas the production of H þ (H 2 O) 3 fragments decreases correspondingly.

Open shell fragments with an electron vacancy
The water radical cation is a highly acidic species which will transfer one proton to neutral water according to whenever possible. In this reaction, spin and charge separate. This process is driven by the difference of proton affinity between H 2 O ($7 eV) and OH • ($6 eV). 6 Therefore, an electron hole in water will eventually yield an OH • and hydrated proton fragments whenever structurally possible. In photoionized protonated water clusters, this channel can be avoided if the system manages to eject a H 2 O þ cation that is not able to transfer its proton to any neighbour monomer. This is possible if the hole localizes in a terminal water molecule that acts as a proton acceptor.  water. The OH • yield of 0.23 is generated from the central water molecule, which leads to two H 3 O þ cations. For HðH 2 OÞ 2þ 6 , either in the Zundel or cage configuration, the larger amount of water monomers results in a larger yield of OH • for various initial excited electronic states considered. The yield of OH • , however, decreases with the level of initial electronic excitation. For the Zundel configuration, this correlates with an increase in the yield of the ðH 2 OÞ þ 2 radical cation, which by the final time of 120 fs has not yet separated into H 3 O þ and OH • . Based on energetic arguments, these two fragments should correspond to be the final products for such trajectories. For cage configurations, the decrease in the OH • yield as the level of electronic excitation increases only weakly correlates with an increase in the yields of H 2 O þ and ðH 2 OÞ þ 2 . Indeed, this corresponds to trajectories, where charge and spin have separated yielding neutral radical fragments of the type OH • (H 2 O) n .
In Fig. 6, the fractions of H 2 O þ and OH • fragments are presented as a function of time in order to characterize the dynamics of the electron hole. As illustrated in Fig. 6, the release of a H 2 O þ fragment progresses rapidly for HðH 2 OÞ 2þ 3 after a localized ionization from one of the external water monomers in the D 0 electronic state. OH • ejection is rare for ionized H þ (H 2 O) 3 , because the external water molecules once ionized as a H 2 O þ cannot transfer a proton to a neighboring neutral H 2 O in its vicinity. However, production of OH • is still possible for excited electronic states when the hydroxyl radical emerges from the central water molecule. The OH • fraction oscillates for ionized H þ (H 2 O) 6 in the extended Zundel conformation, corresponding to the oscillation of the proton between H 2 O þ and H 2 O to finally result in From the lower panel of Fig. 6, one can conclude that the initial cage structure is stable for about 60 fs, much longer than the trimer and hexamer with extended Zundel structure, for which the fragmentation already begins at an earlier time of about 10 fs.  6 . An OH • or H 2 O þ fragment is counted when its oxygen atom separates more than 2.9 Å from other oxygen atoms. The H þ (H 2 O) 21 cluster has a spherical shape with a diameter of $7 Å , comparable with that of C 60 , and a hydronium cation enclosed inside the spherical cage. Contrary to what one might have expected, the configuration of an enclosed hydronium H 3 O þ cation inside a cage of water monomers is not the most stable structure of H þ (H 2 O) 21 in the gas phase. Instead, calculations show that the cluster favors a surface-protonated cage with a neutral water inside, 21,27,29 which is confirmed experimentally. 30 Configurations with a surface excess proton were considered by us in the previous works in which we focused on time-resolved x-ray absorption spectroscopy of photoionized water clusters. 25 Here, the less stable H þ (H 2 O) 21 conformer is considered as a closer model of the solvated proton in water.
After photoionization from an electron with a high ionization potential from the vicinity of the central solvated H 3 O þ cation, the cluster remains largely stable, which can be observed in Fig. 7 by looking at the temporal evolution of the proton-hole distance D ph (t) and the radius of the spherical shell formed by the 20 water monomers. The dynamics of the photoionized HðH 2 OÞ 2þ 21 starts in this case in the outer valence excited state D 62 described at the Koopmans' level of electronic structure theory. In contrast to the protonated hexamer, for which the protonhole distance grows monotonically, the proton-hole distance now oscillates periodically after the initial charge separation. Figure 8 shows the molecular orbital holding the electron hole as a function time for a selected trajectory started from the electronic state D 62 . The electron hole quickly delocalizes to the surface of the cluster, where its charge delocalizes, while the excess proton keeps its position on the central hydronium cation. In the simulation, the dication maintains its spherical structure until the end of the trajectory at t ¼ 200 fs. During this time, the very high level of electronic energy excitation of about 15 eV above the ground electronic state of HðH 2 OÞ 2þ 21 is quickly dissipated into vibrational energy. This can be seen in the plot of all electronic state energies in Fig. 9, where the black line represents the populated electronic state at a certain time for a particular surface-hopping trajectory. The energy span of this quasi-band of electronic states is of about 15 eV for HðH 2 OÞ 2þ 21 , and the number of electronic states corresponds to 3N w , where N w is the number of water monomers, each contributing three outer valence molecular orbitals to the cluster. FIG. 7. Relative proton-hole distance D ph (t) À D ph (t ¼ 0) (upper panel) and mean radius of the spherical shell formed by 20 water mononers (lower panel), when photoionization removes an electron with high binding energy from the vicinity of the hydronium at the center of the sphere. The shell radius is estimated as the mean distance between the central hydronium and the water monomers on the outer shell.

IV. CONCLUSIONS
We have investigated the ultrafast nuclear and electronic dynamics of protonated water clusters H þ (H 2 O) n upon single photoionization by XUV radiation in the photon energy range between 10 and 25 eV, in which outer valence electrons can become ionized, by surfacehopping molecular dynamics simulations.
First, we have focused on the small clusters H þ (H 2 O) 3 , H þ (H 2 O) 6 in extended Zundel configuration, and H þ (H 2 O) 6 in cage-like conformation. After photoionization of these clusters, two positive charges are left in the system, the one initially related to the excess proton and the one related to the electron vacancy. Ionization into excited states of the doubly charged clusters results in ultrafast electronic relaxation driven by non-Born-Oppenheimer effects and taking place in a time scale of 10 to 30 fs. An analysis of the fragmentation products after photoionization reveals various patterns: (1) the higher the electronic excitation of the parent ion, the smaller the protonated water clusters generated by the Coulomb explosion; (2) ionization of a terminal water acting only as a proton acceptor, which is the case for ionization into the D 0 state of HðH 2 OÞ 2þ 3 , results in ejection of a water radical cation H 2 O þ , for which it is not possible to transfer any of its protons to a proton acceptor; and (3) clusters larger than HðH 2 OÞ 2þ have extremely low yields of H 2 O þ because the fragmentation almost always ends up in neutral and radical fragments and two protonated water clusters.
Finally, we have addressed the dynamics of the H þ (H 2 O) 21 cluster upon photoionization. We have not selected the most stable conformation of this cluster, but a local minimum in which the hydrated proton is located in the central water molecule and is practically spherically solvated by the other 20 water molecules. Photoionization of electrons with binding energy of about 25 eV results in a hole initially located on the same monomer as the excess proton. This hole quickly delocalizes to the outer shell of the cluster driven by the Coulombic repulsion between proton and hole. Strong non-Born-Oppenheimer effects are responsible for the ultrafast electronic relaxation of the cluster, which deposits about 10 eV of electronic energy into vibrational motions of the cluster. The mean distance between the proton and the hole increases from 0 to about 3 Å in about 100 fs and the HðH 2 OÞ 2þ 21 cluster executes some breathing dynamics but it remains mostly stable within the 200 fs for which the simulations have been carried on.
The simulations of the HðH 2 OÞ 2þ 21 cluster seem to indicate that photoionization of the water monomer carrying the excess proton leads to a fast delocalization of the electron hole that may result in a delay in the appearance of neutral radical species resulting from the separation between spin and charge. These aspects remain to be better understood on the basis of simulations of larger protonated clusters or of solvated protons in liquid water.

ACKNOWLEDGMENTS
Financial support by the Hamburg Centre for Ultrafast Imaging is gratefully acknowledged. Z.L. is thankful to the Volkswagen Foundation for a Paul-Ewald postdoctoral fellowship. The authors are grateful to Dr. Mohamed El-Amine Madjet and to Professor Robin Santra for stimulating discussions.

APPENDIX: BLOCK-WISE COUPLING SCHEME FOR SURFACE HOPPING DYNAMICS
To treat dynamics of the protonated 21-mer, whose valence band consists of 63 electronic states, we applied a block-wise non-Born-Oppenheimer coupling scheme to the fewest-switches surface hopping calculations. In this scheme, only a sub-block of the matrix of non-Born-Oppenheimer couplings _ R a Á d a JI is calculated at each time step, which corresponds to a certain number of states N / above and below the currently populated electronic state index. In total, hence, 2N / þ 1 are considered at each time step, except when the index of the populated state is closer to the highest or lowest electronic state of the window of electronic states than N / positions. The rest of matrix elements of the non-Born-Oppenheimer coupling matrix are set to zero, assuming non-Born-Oppenheimer interactions with the excluded states have relatively little effects on the evolution of the electronic state amplitudes.
In order to test the block-wise non-Born-Oppenheimer coupling scheme, we first applied it to a one-dimensional model system, which is initiated from the highest energy adiabatic state, and relaxes subsequently to the ground adiabatic electronic state. The electronic dynamics of crossing many parallel potential energy surfaces resembles the non-Born-Oppenheimer relaxation process appearing in the examples discussed in the main text. The model Hamiltonian is given in the diabatic representation as kðjii h0j þ j0i hijÞ ; with parameters chosen to be j ¼ 0.5, d ¼ 0.5, k ¼ 0.05, and N ¼ 9. The model hence consists of 10 states. The derivative coupling matrix elements d JI are numerically calculated from the definition ððr x AÞA À1 Þ JI , 31 where A is the diabatic-adiabatic representation transformation matrix that diagonalizes the diabatic potential matrix. The resulting adiabatic potential curves are shown in Fig. 10. The mass of the particle is chosen to be 2000 a.u. We have launched 2000 trajectories with various block sizes. The population evolution in the adiabatic representation from fewest-switches surface-hopping calculations with increasing block sizes indicates that the populations are practically converged for a block size N / ¼ 3, as shown in Fig. 11. The block-wise non-Born-Oppenheimer coupling scheme was tested as well for the dynamics of the protonated water trimer cation H þ (H 2 O) 3 . In Fig. 12, using population dynamics as a benchmark, convergence is shown with fairly small block size.