Simulating the Photochemical Birth of the Hydrated Electron in Liquid Water

Domain: biology
Published: 2026-03-10
Editorial outcome: minor_revision

Abstract
The photochemistry of UV-irradiated liquid water underlies many physical, chemical, and biological processes, with the formation of the hydrated electron as a central event. Despite extensive experimental and theoretical efforts, its microscopic origin remains incompletely understood. Using excited-state molecular dynamics simulations of photoexcited liquid water, we resolve the sequence of chemical events leading to hydrated electron formation on the excited state. The excitation localizes on specific hydrogen-bond network defects, followed by two competing pathways. The first produces a hydrogen atom and undergoes ultrafast non-radiative decay to the ground state within 100 femtoseconds. The other proceeds via proton-coupled electron transfer, generating hydronium ions, hydroxyl radicals, and an excited-state hydrated electron. This mechanism is driven by ultrafast coupled rotational and translational motions of water molecules, forming water-mediated ion–radical pairs that persist on picosecond timescales and influence visible emission. These results provide a unified framework for interpreting time-resolved spectroscopic observations and guide future experimental and theoretical investigations.

Introduction
The hydrated electron has captivated scientific interest since its discovery nearly 60 years ago[1]. This chemical species plays a fundamental role in diverse areas spanning radiation chemistry[2, 3], DNA damage[4, 5], and redox reactions[6, 7]. One of the mostwidelyemployedmethodstogeneratehydratedelectronsinvolvephotoexcitation across a broad range of photon energies, from UV to γ-ray irradiation[2, 8, 9]. This species has been generated and studied in various aqueous environments, including bulk liquids, clusters in vacuum and at interfaces[10–16]. Interestingly, the hydrated electron can be created upon photoexcitation to the first electronic excited state of liquid water[17, 18]. This is particularly intriguing given that the excitation energy involved( 8.3eV)[19]iswellbelowtheionizationpotentialofwater( 11eV)[19,20]. Perhaps e ∼ ven more surprising is that experimPents have reported elec ∼ tron generation at photon energies as low as 6.5 eV[21, 22].

Understanding how hydrated electroNns are produced upon excitation to the first absorptionbandinliquidwaterrequireIsadetailedcharacterizationandunderstanding of the initial photoexcitation event. Numerous experiments have shown that excita- tion at different energies using eEither single[23, 24] or two photon[17, 25] within the first broad band in liquid watLer (approximately 7-8.5 eV), open different pathways in the excited state. HoweveCr, unraveling excited state dynamics after photoexcitation in liquid water remains a significant challenge due to the extremely short timescales involved (less than 1 ps) and the simultaneous formation of multiple reactive species (e.g. electrons, protons and hydroxyl radicals to name a few)[2, 26, 27]. While exper- iments can successfully identify these final products, which molecular species in the hydrogen-boAnd network get excited and the intermediate steps that occur immedi- ately after photoexcitation in realistic condensed phase conditions, remain largely unresolved.

To address these limitations, a variety of theoretical models and computational techniques have been developed[28]. These efforts aim to complement experimental resultsandprovideinsightintothetransientprocesses underlyingelectrongeneration and solvation in water. A central question in this regard is whether the excitation is localized on a single water molecule or delocalized over multiple molecules. Although thedominantviewsupportslocalization[18,29,30],tothebestofourknowledge,most previous theoretical studies do not capture the coupling between thermal disorder and condensed phase fluctuations, thereby limiting the generality of their conclu- sions. In addition, due to the quantum-mechanical nature of the hydrated electron, one conceptual and practical challenge is how to best localize its position within the hydrogen-bondnetwork.AsnotedbyHerbert,whohaspioneeredthetheoryandsim- ulations of the hydrated electron in the electronic ground state, most studies rely exclusively on visual inspection of the electronic wavefunction or density which can bias the interpretation[28]. This aspect also applies to studies of underlying nature of the electronic excitation in condensed liquid water. Importantly, the hydrogen bond- ing environment surrounding the molecules that are excited, is known to significantly influence both the energy and intensity of electronic transitions[31–33].

Photoionization of water is known to create the solvated electron[25, 34–36]. The equilibrium structure of the hydrated electron remains a subject of ongoing debate. Earlytheoreticalmodelsproposedeitheracavitymodel,wheretheelectronresidesin a vacuum-like pocket surrounded by waters[37–39] or non-cavity models in which no excludedvolumeisformedandtheelectrondiffusestoseveralwatermolecules[40,41]. These earlier methods relied on parameterized interactions between water molecules andtheelectron,whichlimitstheirpredictivepower.Asaresult,moreaccuratequan- tum chemistry approaches have become popular in recent years. Ab initio methods have been applied to water and halide-water clusters in vacuum[42–44], but how and if these results can be extrapolated to bulk water is not straightforward. Notably, some of these models failed to accurately predict the coEordination number of the hydrated electron in the liquid phase[28]. Ab initio bulRk water simulations have been alsoperformedforthehydratedelectron[25,34,45].Whilealloftheseapproacheshave provided valuable insights into the properties of the hydrated electron at equilibrium in the electronic ground state, they do not captur e the initial stages of its formation. This is due to the fact that in most of these studies, the simulations start with a pre-inserted electron in the electronic grounId state and therefore do not capture the events on the excited-state.

In this work, we investigate the excited-state dynamics of neat liquid water fol- lowing excitation to the first electroLnic excited state using the Restricted Open Kohn Sham (ROKS) approach[46, 47] a technique that has been successfully used to inves- tigate dissociative and chargIe transfer excitations in closed-shell systems[48–51]. In addition, we perform an eTxtensive validation of this method in the context of water photochemistry, benchmarking its performance against different electronic structure settings and alternative excited state approaches. By performing simulations on an ensembleofthermallysampledconfigurationsfromtheelectronicgroundstate,weelu- cidate new insights into the photophysics and photochemistry of water paying special attention to the important role of fluctuations of the topology of the water hydrogen- bond network. Our study begins by analyzing the fundamental properties of the first excitedstateinliquidwater,addressingthelong-debatedquestionoflocalization.The canonical manner of tackling this in quantum-chemistry is through visual inspection of the electronic wavefunction or electronic density. Here, using the Inverse Partici- pation Ratio (IPR)[52, 53] which has been previously applied to study the extent of delocalizationofexcitonwavefunctionsinmolecularsystems[54],wedirectlyprobethe partitioningoftheelectronicspindensityamongthephotoexcitedwatermolecules.In addition, the analysis of the local environment of the excited water molecules on the red-tailrangeofthespectrumrevealthatmostofthemresideontopologicaldefectsin thehydrogenbondnetwork[29,31,55]wheretheasymmetryofthehydrogenbonding also plays a critical role.

Our theoretical predictions provide a coherent and unified understanding of the mechanisms underlying the photochemistry of water drawn from several different experimental time-dependent spectroscopic measurements over the last two decades. Specifically, our simulations provide new insights into the two proposed mechanisms invoked to rationalize the excited-state dynamics of water, namely Hydrogen Atom Transfer (HAT) and Proton Coupled Electron Transfer (PCET). For both HAT and PCET, one fundamental challenge has been identifying which transient species form after absorption of a photon on the excited state itself, in contrast to being formed as a consequence of photodeactivation to the ground-state. Our approach overcomes these hurdles predicting that the HAT mechanism exhibits a higher quantum yield, consistent with experimental observations[17]. Additionally, we observe the formation ofatransienthydroniumradical,aspeciespreviouslyproposedbutnotyetexperimen- tally confirmed[56, 57]. Importantly, our results also show that this mechanism does not lead to the formation of a hydrated electron in the excited state. On the other hand,forthePCETpathway,weobservetheformationandearly-timeevolutionofthe hydrated electron before the system relaxes to the grEound-state. The creation of the hydrated electron in the excited-state is in turn dRriven by coupled translational and rotational motions that lead to hydronium ion and hydroxyl radical pairs mediated by solvent as recently suggested by ultrafast electron diffraction measurements[35]. Finally, we demonstrate that the photon em ission by the hydrated electron measured by Tauber et. al.[58] is strongly modulated by the extent of the electron localization. These results offer new perspectives iInto tuning the color of fluorescence emission of excited state hydrated electrons in a queous systems.

Results
2.1 Initial photo-Iabsorption in liquid water
Understanding the nature of the initial photo-absorption in liquid water is a criti- cal step toward elucidating the mechanisms that lead to hydrated electron formation, especiallyatAexcitationenergieswithinthefirstabsorptionband.Figure1summarizes ourmainfindingsregardingthecharacteristicsofthisinitialexcitationprocessemploy- ing 100 different conformations sampled from well converged deep neural-network based ground-state simulations (see Methods section for details). Figure 1A presents the experimental absorption spectrum compared with the probability density func- tion (PDF) of the excitation energy for the S S transition, determined using the ROKS method for all the conformations (see Methods section). Our simulations pre- dict an average excitation energy of 8.1 eV, which is in excellent agreement with the experimentalvalueof8.3eV[18,19]andalsowithTime-DependentDensityFunctional Theory (TDDFT) (see Supplementary Note 1).

To quantify the degree of the localization of the first electronic excitation, we par- titioned the electronic spin density across individual water molecules of the system and computed the Inverse Participation Ratio[59] (IPR). This approach provides an agnostic manner of determining the extent to which the excitation is delocalized (see Methodssection).AsshownintheFigure1B,themajorityofconfigurationsinvolvea single water molecule in the S S transition, although a non insignificant fraction exhibits a more delocalized excitation spread over up to five water molecules. These eventstypicallyinvolveachainofconnectedwatermoleculesformingawater-wire[60– 62], similar to what was recently observed in many-body theory electronic structure calculations in ice[63]. The typical analysis used to characterize the extent of delo- calization rely on visual inspection of molecular orbitals, electronic densities, or spin densitiesrepresentedasisosurfaces,whichcanbehighlysensitivetothechosenisovalue andthussomewhatarbitrary[28].Ourapproachovercomestheselimitationsproviding a quantitative description of the excitation which allows also for fractional contribu- tionsofthenumberofwatermolecules.Inaddition,previousstudieshavearguedthat thefirsttransitionislocalizedonasinglewater,howevermostofthemwerelimitedto isolatedmonomersorsmallwaterclusters[29,31,32],wherethesystemsizeinherently restricts delocalization. Other studies have employed periodic water systems of vary- ing boxes sizes, inferring localization from the weak dependence of excitation energies on the system size[30]. While these approaches often produce excitation energies and oscillator strengths that align well with experimental results, they generally rely on limited sampling and thus fail to capture the influence of tEhermal fluctuations on the excitation process. Understanding the degree of electrRonic delocalization upon exci- tation to the first band is particularly important, as one and two photon absorption experiments can access different electronic configurations[17, 23, 25]. These configu- rations, in turn, may lead to distinct pathways in the excited-state dynamics of liquid water.

A comprehensive understanding of the first electronic transition in liquid water requires not only the characterization of electronic density but also the analysis of thesurroundingenvironment.Inthiscontext,thedirectionalityofthehydrogen-bond (HB) network around the water molecules involved in the transition is particularly important.TheanalysisshowninFigure1Benablespreciseidentificationofthewater molecules participating in the excitation, allowing us to classify them according to theirHBenvironment.Figure1CpresentstheprobabilitycountsofHBdefectsacross theensemble,alongwithaschematicrepresentationoftypicalanddefectivehydrogen- bonding motifs (see Methods section for details).

Our results clearly show that the water molecules most frequently involved in the lowestelectronictransitionarethosefoundinhydrogenbond(HB)environmentswith defects, particularly those with at least one missing acceptor hydrogen bond. These include configurations such as AD (a water molecule accepting and donating a single hydrogenbond)andADD(awatermoleculedonatingtwoandacceptingonehydrogen bond), as illustrated in Figure 1C. This observation is rooted in the nature of the electronictransition,whichcorrespondstoatransitionfromanon-bondingorbital(n) toananti-bondingorbital(σ ).Inawell-formedHBnetwork,thenon-bondingorbitals arestabilizedthroughtheirparticipationinhydrogenbonding.Thisstabilizationraises the energy required for excitation[29, 31]. In contrast, when a water defect does not acceptahydrogenbond,thenon-bondingorbitalsarenolongerstabilizedbythefield fromthenearbyprotonandthereforeremainathigherenergy.Asaresult,theenergy gap for the n σ transition becomes smaller. These arguments are essentially also at the root of the blue shift of about 0.7 eV observed in the first absorption band of water when comparing the liquid phase to the gas phase[18, 19].

2.2 Photo-chemical mechanism of water deactivation
The excited state dynamics of liquid water reveal two distinct decay pathways: an ultrafast process with an average lifetime of 25 fs, and a slower one of 400 fs (see Supplementary Figure 4). Our protocol, described in detail inSthe Methods section, defines the decay time as the instant when the energy gap between the S and S states drops below 0.2 eV. This approach has been successfully applied to describe excited state relaxation pathways in other systems[64, 65] and, as we will show later, leads to insights that reinforce and enrich the interpPretations of existing experimental observations.Beforeanalyzingthefullensemble,w eillustratethesemechanismsusing representative trajectories in order to build ouNr intuition on the underlying processes.

Figure 2 depicts the two decay pathways: Hydrogen Atom Transfer (HAT, upper panels)andProtonCoupledElectronTran sfer(PCET,lowerpanels).TheHATmech- anismexhibitsanultrafastdecay.AsshEowninFigure2A,whichpresentstheevolution of the S and S potential energLy surfaces, the system reaches the non-radiative S S crossingatapproximatCely10fsinthistrajectory.Duringtheinitialfewfem- toseconds, one water molecule undergoes O–H bond dissociation. This is evident in Figure 2B, which shows the coordination number of the dissociating hydrogen over time which goes from 1 to 0 (see Methods section). This behavior is consistent with the strongly dissociative character of the S (n σ ) potential energy surface and agrees with previoAus theoretical and experimental findings[18, 19, 66]. To determine whether the dissociating species is a hydrogen atom (H ) or a proton (H+), we deter- mined the center of the electron and its distance to the dissociating hydrogen (see Methodssection fordetails).As shownin Figure2C,this distancerapidlydrops from 1.2 to below 0.2 ˚A, clearly indicating that the observed process involves hydrogen atom dissociation.

The PCET mechanism displays a slower decay, on the order of hundreds of fem- toseconds(seeFigure2D).Infact,2outof100trajectoriesdonotexhibitnon-radiative decay to the S state, remaining on the S for the entire 1 ps of the simulation. Fol- lowing photoexcitation of a water molecule, in addition to the OH dissociation, an electron is released into the surrounding hydrogen-bond network. Figure 2E shows the evolution of the gyration radius of the electron over time (see Methods section for details), which serves as a proxy for its spatial delocalization. Initially, the electron is highly delocalized, but around 100 fs becomes increasingly localized. This behavior is consistent with experimental observations[34, 67], which report a transition from a verydelocalizedelectron(gyrationradiusabove20˚A)toamorelocalizedwetelectron, and finally to the fully hydrated electron in the ground state (gyration radius 2.4 ˚A).

A key feature of this mechanism is that it also involves an ultrafast O–H bond dissociation, consistent with the dissociative character of the S (n σ ) potential energy surface. However, unlike the H atom pathway, here the dissociated species is a proton (e.g., H+ or H O+). This is evident from Figure 2F, which shows that the distance between the center of the electron and the closest hydrogen. Initially, this distance decreases (within the first 100 fs) but in contrast to the HAT mechanism, the energy gap is still high (around 2 eV, see Figure 2D). After this period of time, this distance subsequently starts to increase supporting the notion of the formation of the hydrated electron and the H O+. This also implies that the dissociated water molecule that lost the original proton is now an HO• species. Later we will delve into the role of the H O+ and HO• on the photophysical properties.

Using the features described in Figure 2, specifically the gyration radius, dis- tancebetweentheelectroncenterandtheclosesthydrogen,wecanclearlydistinguish betweenthetwodecaymechanisms.Thedetailedprotocolforthisclassificationispro- vided in the Methods section. While a comprehensivEe discussion of each mechanism across the full ensemble of trajectories is presenteRd in the following two sections, it is worth noting that our results indicate that approximately 53% of trajectories follow the HAT pathway, while the remaining exhibit the PCET mechanism. Experimental studies report a 45% occurrence of the HAT pathway at 6.7 eV excitation and 70% at 8.4 eV[68, 69]. Our result lie between these two values, in excellent agreement with the fact that our simulations probe thIe red edge of the first absorption band (8.1eV).

2.3 H Atom Transfer Mechanism
To fully characterize the hydrogen atom transfer mechanism, we constructed a two-dimensional probability density function (PDF) using two key descriptors: the coordination number of the dissociating hydrogen atom (CN ), which always carries theelectron,andthecoordinationnumberofitsnearestoxygenneighbor(CN ).The analysis (see Methods section for details), based on all trajectories that follow this mechanism, is presented in Figure 3.

The most probable configuration, with CN = 1 and CN = 2 (Figure 3A), correspondstotheinitialstateofthesimulation,wherethehydrogenatomisstillpart of a neutral water molecule. Upon photoexcitation, we observe a fast dissociation of this hydrogen atom, accompanied by a sharp drop in the S S energy gap from approximately8eVtobelow0.2eV,leadingtoanon-radiativetransition(asindicated by all the red points in the Figure 3). In some trajectories, the dissociating hydrogen atom moves into an empty cavity (Figure 3B), where its nearest oxygen is part of a hydroxyl radical (HO•), resulting in CN =1and leadingto anon-radiative crossing to the ground state.

Another frequently observed pathway is the formation of a species in which the dissociatinghydrogenatombindstoawatermolecule,formingtheradicalhydronium (H O•),characterizedbyCN 1andCN 3(Figure3C).Weconfirmedthatthis speciesisindeedH O• andnott ∼ heconvention ∼ alH O+ ion,sincetheelectronremains associated with it. While the species H O• has been predicted by various theoretical models[28, 56, 57, 70], our work is the first to demonstrate its formation following excitation by simulating the entire process from photoabsorption to relaxation in liq- uidwater.Accordingtooursimulations,theradicalH O• appearsonlytransiently.It eventuallyundergoesfurtherdissociation,wherethehydrogenatomcarryingtheelec- tronmigratesintoaregionsurroundedbywatermolecules.Intheseconfigurations,an OH bond from a nearby water molecule points toward the hydrogen atom (see Figure 3D) creating an unstable geometry ultimately leading to non-radiative decay to the ground-state. To the best of our knowledge the radical hyd ronium (H O•) has not been observed experimentally through photoionization/photoabsorption of water per- haps in part due to the fact that it is formed fleetingly as a rare transition state-like motif.

Aninterestingobservationisthatapproximately20%ofthetrajectoriesfollowthe mechanism involving the formation of a hydrogen atom in an empty cavity (Figure 3B) with a lifetime of 13 fs, while the remaining 33% correspond to the alternative HATpathway(Figures3Cand3D)exhibitingalongerlifetimeof34fs.Theseresults correlate well with the defect analysis presented in Figure 1C. When the dissociat- ing water molecule has a missing hydrogen-bond donor (AD or AAD defect), there is sufficient space available to allow an ultrafast hydrogen atom dissociation into a cavity ( 13 fs), in excellent agreement with values reported for a hydrogen-bonded water dimer in vacuum[71, 72]. In contrast, when the dissociating molecule lacks a hydrogen-bondacceptor(ADDdefect),thesurroundinghydrogen-bondnetworkinthe condensedphaserestrictsthismotion,delayingthenon-radiativedecayandpromoting the formation of the transient H O• species.

2.4 Proton Coupled Electron Transfer Mechanism
Ourresultsindicatethatthispathwayleadstotheformationofthehydratedelectron while the system is still in the excited state. Figure 4 shows data compiled from all the frames and all the trajectories that follow this mechanism after photoexcitation of liquid water. Panel 4A presents a two-dimensional probability density plot of the gyration radius of the electron versus the energy gap between the S and S states. As described previously for a representative trajectory, the electron is initially highly delocalized following excitation, with an average gyration radius of approximately 6 ˚A. As time progresses, the gyration radius decreases and stabilizes around 2.7 ˚A, indicating localization (see Supplementary Figure 8).

Experimental studies have reported initial gyratioSn radius of around 20-40 ˚A[34, 67] immediately after photoabsorption. Although our method captures the same qualitative behavior, direct comparison with experiment is limited by the size of the simulation box (12.42 ˚A), which restricts the maximum observable extent of delocalization[73]. The final gyration radius wPe observe in the excited state is 2.7 ˚A, slightly larger than the experimentally repor ted value of 2.4 ˚A for the hydrated elec- tron in the ground state[74]. Note that Nthe decrease of the gyration radius from 6 to 2.4 ˚A incurs a decrease in the energy gap. Once the electron has localized, the sys- tem undergoes a non-radiative trans ition to the ground state, as indicated by the red points in Figure 4A.

As mentioned in the introLduction, generating a free electron in bulk water using photon energies within thCe first absorption band requires some degree of nuclear rear- rangement to release the electron from the molecule[17]. Earlier, we showed that this release occurs simultaneously with the dissociation of the OH bond in the excited water molecule, leading to the formation of a free electron, a hydroxyl radical (HO•), andaproton(H+ orH O+).Wenowmoveontoelucidatingthatthelocalizationofthe electron andAthe formation of the other species requires a collective process involving reorganization of the hydrogen bond network of the water molecules.

In Figure 4B we illustrate the coupling between the gyration radius against the angle formed between the vector connecting the oxygen and the center of the electron andtheOHvectorofthewatermolecule.Forthisanalysisweonlyincludethenearby water molecules involved in the first solvation shell of the electron, (see Methods section for details). When the electron is delocalized, this angle spans a broad range from 0 to 120 , reflecting the absence of a well-defined solvation shell. As the system relaxes, the angle narrows to below 30 , consistent with the orientational preferences seen in hydrated electron solvation shells[58, 73, 75]. These results indicate that the structure of the hydrated electron forms within the excited-state itself, before the non-radiativedecayconsistentwithrecentexperimentsusingsub-two-cyclevisible-to- shortwave infrared pump–probe spectroscopy[76]. Analysis of the radial distribution function in the low gyration radius regime on the excited state (see Supplementary Figure 9), where the hydrated electron is effectively localized, reveals that the first solvationshellcomprisesbetween4-5watermolecules.Thisisconsistentwithprevious findings for the hydrated electron in the ground electronic state[25, 45, 51].

2.5 Coupled Translational and Rotational Solvent Dynamics
Stabilizes Water-Mediated Photoproducts
The previous analysis does not directly capture the solvent dynamics that drive the transition toward a localized hydrated electron. To better understand the molecular rearrangements leading to the S S crossing configuration, we analyzed both the rotational and translational motion of water molecules during the excited-state evolution.

Figure5Ashowsthetimeevolutionofthecosineanglebetweeneachdipolevector ofthewatermoleculesanditsdipoleorientationattheendofthesimulation,wherethe hydrated electron is fully localized. This measure quantifies how much each molecule reorients relative to its final configuration. We distinguish between water molecules that ultimately form the first solvation shell of the electron and those that remain fartheraway.Forthepopulationofwatersthatareeventuallyhydrogenbondedtothe excess electron, we observe two phases of the dynamics. In the first part that occurs withinthefirst400fsthedipolerotatesby 20degrees,thisleadstoareductionofthe radiusofgyrationfrom5to4˚A.Theseinitialtimescalefundamentallycorrespondsalso to the creation of the initial cavity that will host the excess electron. On a timescale beyond 500 fs, further rotation of these water molecules leads to the final localization of the excess electron.

In addition to these rotational dynamics, Figure 5B shows that localization also involves collective translational motion. Water molecules in the first solvation shell move on average by approximately 1.2 ˚A relative toEtheir final positions, a displace- mentslightlylargerthanthatofthemoredistantRwatermolecules.Overall,ourresults showthatbothrotationalandtranslationalrearrangementsofspecificwatermolecules are necessary for the localization of the hydrated electron. These collective motions also cause a notable distortion of the hydrog en bond network, consistent with experi- mental findings from ionization studies in liquid water[35]. A similar reorganizational responsehasbeenreportedinabinitioImoleculardynamicssimulationsofcarbondiox- ide reduction in water[77, 78], wher e the breaking of hydrogen bonds in the second solvation shell triggers the stabilization of the excess electron on the carbon dioxide molecule. This process also inLvolves coupled translational and rotational motions of water molecules, closely pCaralleling the behavior observed in our simulations.

As alluded to earlierI, the creation of the excess electron after photoexcitation also leads to the formatTion of a hydronium ion (H O+) and a hydroxyl radical (HO•). To investigate the behavior of these species during the excited-state dynamics, we presentinFigure5Cthecouplingbetweenthegyrationradiusoftheelectronandthe distancebetw A eentheHO• andtheH O+.Interestingly,ouranalysisrevealsthreedis- tinct minima in the HO•(H O+) distance, particularly in the region of low gyration radiuswheretheelectronbecomeslocalized.Theseminimaoccurat 2.6,4.3,and6.2 ˚A. Our simulations are in excellent agreement with very recent time-resolved electron diffractionexperimentsfollowingtheionizationofliquidwater,wherethecorrespond- ing minima were found at 2.4, 3.4, and 5.8 ˚A[35]. These three minima correspond to situations involving the initial formation of contact ion-radical pair and then the sep- arationmediatedbythesolvent.Althoughprevioustheoreticalstudieshaveidentified the HO•(H O+) pair by simulating ionized water with the electron removed[79–81], our findings show for the first time, the creation of these correlated structures upon photoexcitation starting from the neutral bulk water.

The creation of ion-radical pairs at different distances has important implications on the optical properties. To probe this, we show in Figure 5D the energy gap as a function of the HO•(H O+) distance. For the contact ion-radical pair at 2.6 ˚A the average energy gap is 3 eV on average. As one moves towards the creation of the solvent-separatedion-radicalpairs,theenergygapdropsto 2.0eV.Thistrendisalso consistent with experimental observations, which demonstrated that the formation and spatial separation of the HO•(H O+) pair occur before the internal conversion associated with the p s transition of the hydrated electron[35].

2.6 Spontaneous Emission from Hydrated Excess Electrons
An important and still not fully understood property of the hydrated electron in liquid water is its photon emission following excitation. Understanding this emission can provide valuable insight into the nature of the excited state and the relaxation dynamics of the system. Although experimental studies are relatively scarce, Tauber and Mathies have explored this optical behavior through fluorescence and Raman spectroscopies[58].Motivatedbytheirfindings,wenowfocusonanalyzingtheemission spectrum obtained from our simulations and examine how it depends on the gyration radius of the hydrated electron.

Figure 6A shows the experimental emission spectrum along with our calculated S S energy gaps, evaluated for configurations with both large and small electron gyration radius. As the excess electron becomes more localized, transitioning from a larger to a smaller gyration radius, the energy gap decreases from an average of 2.4 to 1.3 eV (516 to 954 nm). Interestingly, the experimental emission spectrum is in excellent agreement with our predictions in both peak position and spectral width, especially for configurations with smaller gyration radius. This observation supports previousinterpretationsfromexperimental[58]andtheoreticalstudies[82,83],suggest- ing that emission occurs from a broad ensemble of transient excited state geometries rather than from a single well-defined minimum. It is important to highlight that our theoretical spectrum reflects the recombination of the electron, since we describe the final electronic ground state as a closed shell (see Methods). In contrast, the experi- mentalmeasurementscaptureemissionfromtheexcitedstaEteofthehydratedelectron without identifying the final state[58], which could inRvolve either recombination or hydrated electron in the ground state. Therefore, while the agreement between our results and the experimental data is encouraging, it should be considered prelimi- nary. A more detailed investigation of these opti cal properties is necessary and lies beyond the scope of this study. However, our findings provide a valuable framework that could guide future experimental and thIeoretical efforts aimed at confirming and better understanding the emission behavi or of the hydrated electron.

In order to ensure that these features are not sensitive to the choice of our ROKS method used for the simulation ofLthe excited state, we computed the emission spec- trumforthelast300fsofareprCesentativePCETtrajectoryusingtheenergygapsand oscillator strengths obtainedIfrom TDDFT (Figure 6B). While the spectrum is not fully converged due to limTited sampling, the resulting peak at 1.7 eV closely matches ROKS calculated and experimental emission peak. Additional validations for both mechanisms are provided in the next section.

2.7 Evaluating the Accuracy of ROKS for the Photochemistry
of Liquid Water
In this work, we employed the ROKS approach to investigate the photochemistry of liquid water. While our results show excellent qualitative and quantitative agreement with a broad range of experimental observations, it is important to assess the quality of the electronic structure method and to evaluate its performance relative to more extensively tested theoretical approaches in photochemistry.

In order to first assess the robustness of ROKS with respect to the choice of exchange–correlation functional and basis set, we performed additional calculations using the CAM-B3LYP functional[84] and the TZV2P basis set[85]. The results, presented in Supplementary Note 9, demonstrate that the ROKS method provides consistentenergygapsandreproducesthesamebehaviorforbothphotochemicalpath- ways (HAT and PCET), confirming its reliability across different electronic structure settings.

WefurthercomparedourresultswiththeTDDFTmethod,whichhasbeenwidely applied to the study of excited states and in non-adiabatic dynamics. In addition, for a hydrogen-bonded water dimer in vacuum, we benchmarked our approach against both multiconfigurational methods and TDDFT. The main results are summarized in Figure 7.

Panels 7A and 7B show the evolution of the energy gap between the first excited state(S )andthegroundstate(S )fortworepresentativetrajectories(sampledfrom the ROKS dynamics), corresponding to the HAT and PCET mechanisms, obtained usingbothROKSandTDDFTwiththesamefunctionalandbasissetemployedinthis work (see Methods section). At the initial step, both methods yield nearly identical energy gaps, as further confirmed by a broader validation across all initial conditions (see Supplementary Figure 1).

The main difference arises in the HAT mechanism. At the end of the simulation, ROKS predicts a crossing point with an energy gap smaller than 0.2 eV, whereas TDDFT yields a value near 1.5 eV for the same geometry. Despite this quantita- tive difference, both methods describe the same physical process: water dissociation leading to the formation of HO• and H . In contrast, the PCET mechanism shows excellent agreement between ROKS and TDDFT. In this case, both methods predict the dissociation of a water molecule into HO• and H O+ within the first 20 fs, with comparable energy gap profiles. Moreover, TDDFT also reproduces the localization of the hydrated electron, consistent with the ROKS results, as shown in Figure 7C, which compares the ROKS spin density and TDDFT LUMO at 250 fs.

To further understand the discrepancy between ROKS and TDDFT observed for the HAT mechanism, we examined a model system consisting of a hydrogen-bonded water dimer in vacuum [71, 72]. Segarra et al. computed the Minimum Energy Path (MEP)forthephotodissociationofonewatermoleculefromtheFranck–Condonregion totheS /S conicalintersectionusingthemulticonfigurationalCASPT2method(see Figure7D).Usingboththeground-stateandconicalintersectionoptimizedgeometries at this high level of theory, we computed the potential energy surfaces along the Linear Interpolated Path (LIP) with ROKS (Figure 7E) and TDDFT (Figure 7F) employing the CAM-B3LYP functional. It is worth emphasizing that the numerical values along horizontal axes for the MEP and LIP coordinates do not correspond to identical nuclear configurations, except for the initial point shown in Figure 7D-F. We thus focus on the overall trends of the ground and excited state potential energy surfaces.

ROKS predicts a dissociative pathway along the surface, similar in slope to CASPT2. The main difference lies in the position of the crossing point: at a HO•–H distanceofapproximately2˚AforROKS,comparedwith2.8˚AforCASPT2.Although an O-H distance of 2 ˚A could be considered as dissociated geometry, the topology of the conical intersection differs between t he two methods. Despite this difference, ROKS reproduces the same essential physical process, namely the hydrogen atom dissociation. In contrast, TDDFT prIedicts a local minimum on the S surface at around2˚Awithanenergygapof 1 .5eV,whichisconsistentwiththedifferencesseen in the condensed phase simulations for the HAT mechanism (Figure 7A). Moreover, TDDFT exhibits a 1 eV Lbarrier from this minimum to the conical intersection, which could significantlyCdelay, or even prevent, reaching the intersection within the timescales ( 10 fs) repIorted in previous studies [71, 72].

Taken all togethTer, these results provide a comprehensive validation of the ROKS approach.WedemonstratethatROKScapturestheessentialphysicsofboththeHAT and PCET photochemical pathways in liquid water. Based on the gas phase water dimermodel,ROKSdoesnotreproducetheconicalintersectiongeometryfortheHAT mechanism. In particular, a difference of approximately of 0.8 ˚A is observed in the HO•-H distanceattheconicalintersectionwhencomparedwithCASPT2.Neverthe- less, ROKS qualitatively describes the same photodissociation pathway, reproducing boththeslopeoftheexcited-statepotentialenergysurfaceandthepresenceoftherel- evant chemical species along the dissociation coordinate. While such differences may lead to deviations in the predicted HAT lifetimes in vacuum, we expect their impact to be reduced in the condensed phase. Indeed, in bulk liquid water, ROKS yields lifetimes and quantum yields that are in good agreement with experimental observa- tions. Overall, the extensive validation presented here demonstrates the robustness of our methodology and reinforces confidence in its predictive capability for modeling the excited-state dynamics and photochemistry of liquid water. We nonetheless cau- tion that more systematic studies with multireference methods for condensed phase aqueous systems needs to be explored in the future.

Methods
4.1 Ground-State Molecular Dynamics
We began by constructing a periodic simulation box containing 64 water molecules with a side length of 12.42 ˚A, corresponding to a density of 1 g/cm3 at room tem- perature. An initial equilibration was performed using an empirical water model, SPC/E[95] with the GROMACS software[96] within the NVT ensemble for 100 ns. The canonical-samplingvelocity-rescaling[97]thermostatwasappliedwithacouplingtime of 0.1 ps.

From this classical equilibrated trajectory, we selected 100 configurations, sepa- ratedby1ns,asinitialconditionstoperformabinitioqualityground-statemolecular dynamics simulations with a machine-learning interatomic potential (MLIP). Specif- ically, we chose an MLIP recently trained and developed using the DeePMD-kit architecture[98] on data computed using Density Functional Theory (DFT) with the SCAN0 functional[99], which incorporates 10% exact exchange tSo mitigate the self- interaction error. Prior studies[99] have shown that SCAN0 provides an excellent descriptionofthestructural,dynamical,andelectronicproperties(suchasIRspectra, electronic band structure, and band gap) of bulk liquid water compared to the orig- inal SCAN meta-GGA functional[100, 101]. Each of tRhe 100 simulations conducted withtheMLIPpotentialwasrunfor1nstogeneratPewell-convergedanddecorrelated initialconditions(positionsandvelocities),leveragingDeePMD-kit[98]incombination with the LAMMPS package[102].

All in all, our ground state sampling protocol ensures that we have statisti- cally uncorrelated initial conditions (nuclear positions and velocities) for the water hydrogen-bondnetwork,whichareesseEntialforaccuratelydescribingtheexcitedstate dynamics.

4.2 Excited-State Molecular Dynamics
TheexcitedstatemoleculaTrdynamics(ESMD)simulationswereinitiatedfromthefirst electronic excited statRe (S ), computed using the Restricted Open-shell Kohn–Sham (ROKS)method[46,47,103].Thismethodprovidesanalternativeapproachforcalcu- lating singlet excited states in closed-shell systems that are described by non-Aufbau electronic configurations. By minimizing the energy expression E = 2E E , where E ,E ,E correspond to the energies of the singlet, triplet, and mixed states, respectively, the method effectively accounts for a two-determinant nature of the excited state. Consequently, ROKS offers a computationally efficient and reliable framework for describing dissociative and charge-transfer excited states[48–50, 104].

To propagate the system on the S potential energy surface, we employed the AtomicSimulationEnvironment(ASE)[105],apythonmodulewhichisinterfacedwith CP2Ksoftware[106].Ateachtimestep,wecalculatedtheenergyontheelectronicground stateS usingDensityFunctionalTheory(DFT)toaccountforpossiblenon-radiative relaxation pathways[49, 51] (see below) which in our current setup is determined by examiningtheenergygap.ThesimulationswereperformedintheensembleNVEusing a timestep of 0.5 fs and velocity verlet algorithm to perform the nuclei propagation of the system.

Both excited and ground state calculations were performed using the hybrid PBEh(40)-rVV10 functional[107] in combination with the DZVP-MOLOPT-SR basis set[85]. This level of electronic structure theory has been successfully applied to the study of liquid water[108–110], demonstrating its accuracy in reproducing ionization potential, redox potentials, electronic levels of excess electron. Excellent agreement with the experiments for aqueous iodine solutions were reported by Carter et. al.[51] within ROKS framework and by Lan et. al.[111] employing a similar approach with the same functional.

More recently the functional PBEh(40) was identified as one of the most reli- able choices for describing the hydrated electron in its electronic ground state[73]. In our current work, the favorable agreement with reproducing the experimental UV- absorption and other quantities sensitive to the excited states of liquid water gives us confidence in the quality of our results.

Each ESMD simulation was propagated until the system reached a crossing point, defined as the step where the energy gap between the first excited state (S ) and the ground state (S ) dropped below 0.2 eV. This Ecriterion provides a robust and efficient way to identify non-radiative decay evenRts and allows us to characterize the distribution of deactivation pathways, offering a lower bound on the excited-state lifetimes[64, 65]. A known limitation of this approach is that it does not explicitly includenon-adiabaticcouplingvectors(NAC Vs),whicharekeyquantitiesfordescrib- ing electronic transitions during non-radiative decay. However, a previous study on the radiolysis of water cluster in vacuIum showed that NACVs increase significantly as the energy gap decreases[112]. T his strong correlation implies that small energy gaps are reliable indicators of high non-adiabatic transition probabilities. Therefore, by using the energy gap asLa decay criterion, we effectively capture the relevant non-radiative dynamics, jCustifying the omission of explicitly determining NACVs in our simulations[113, 11I4]. To further support this approximation, we evaluated the Landau-Zener (LZ)Ttransition probabilities[115], which in addition to the energy gap criterion, also account for the curvature of the potential energy surfaces (see Supple- mentary Note 10). This approach has been shown to reproduce the results consistent with those obtained using explicit NACVs in non-adiabatic dynamics[115–119].

Thesimulationsthatdidnotexhibitadecaywerecontinueduntilamaximumtime of 1 ps was reached. Spin density cube files were calculated and outputted every 5 fs during all ESMD trajectories and at the step of crossing, for post-processing analysis (seesectionbelow).Alltheelectroniccalculationsinthepresentworkwereperformed using CP2K version 2024.2[106].

4.3 Absorption and Emission Spectra
The absorption spectrum to the first band and emission spectrum presented in this work (Figures 1A and 6A, respectively) were obtained by calculating the energy gap between S and S using the methods described in the previous section. Because the CP2Ksoftwaredoesnotprovidethetransitiondipolemomentswithinourapproach,the spectra were represented as histograms of the calculated energy gaps. This approach was validated through the calculation of energy gap and oscillator strength using TDDFT with the same functional and basis set, as shown in Supplementary Note 1 for the absorption spectrum and in Figure 6 for the emission spectrum.

4.4 Coordination Number Analysis
To quantify local structure and bond-breaking/formation occurring within the hydrogen-bondnetwork,wedeterminedthecoordinationnumbersofhydrogen(CN ) including only the oxygen atoms in close vicinity and the coordination number of the oxygen (CN ) including only hydrogens in close vicinity. We employed the following smooth switching function:

CN = 1 / (exp[κ(r r )]+1)
Here, i is the atom index, while κ = 10, r = 1.38 ˚A are pSarameters taken from previous works used to study the dissociation of water in the ground-state to form hydroniumandhydroxideions[120–123].Inthismanner,foranisolatedhydrogenatom not attached to any water molecule, the CN is zero while for a hydrogen covalently bonded to a water molecule CN is 1. In the case oPf CN , it takes on a value of 3, 2 and 1, for a hydronium (ion or radical), neutral m olecular water and hydroxyl species respectively.

4.5 Hydrogen Bond (HB) Netw ork Classification
HB were identified using a geometrical criteria introduced by Kumar et. al.[124]: O- O distance of 3.5 ˚A and H-O—O angle 30 deg. A water molecule with four HB is classified as part of a normal HB network, otherwise, it is labeled as a defect (see Figure 1C for categoriesI). Sensitivity analysis to the choice of these geometric parameters is addressed iTn Supplementary Note 3 showing that while the hydrogen- bondcriterianaturallyRaffecttheexactpopulationsofdefects,thetrendsandphysical insights inferred are robust.

4.6 Spin-Density Analysis for Excitation Localization
Toidentifywhichwatermoleculewasinitiallyexcited,weanalyzedtheelectronicspin density for all the 100 initial configuration. We decomposed the total spin density by integrating around each oxygen atom within a radius of 1 ˚A.

s = ρ(r)dr
where ri is the position of the i th oxygen atom and R = 1 ˚A. To estimate the degree of localization, we computed the Inverse Participation Ratio (IPR)[52, 53]:

IPR 1 = ( s2) 1
Thismetricprovidesacontinousmeasureofhowmanywatermoleculesparticipate in the excitation. The distribution of IPR 1 values for all the initial configurations is shown in Figure 1B. The sensitivity analysis and validation to the choice of the integration radius parameter are discussed in Supplementary Note 2.

4.7 HAT versus PCET Mechanism
Photo-excitation of liquid water follows either one of these two mechanisms:

H O HO +H (HAT)
2 • •
H O hν HO +H++e (PCET)
2 • −

A direct calculation of the electron center from the spin density cube files yields a position located between the HO• and the H (Equation 4), or between the HO• and the excess electron (Equation 5). Supplementary Figure 5 illustrates this point. This result lacks a clear physical interpretation, as it does not correspond to the actual location of any species and essentially implies that the spin density is located close to the geometric center of the species involved. ThSis problem does not arise in earlier theoretical studies[25, 28, 45, 75], as those simEulations typically start with a pre-inserted electron in bulk water and the formation of an additional radical species are not observed.

Toovercomethischallenge,wefirstidentifytheHO• speciesbylocatingthewater molecule with the longest OH bond in each trajectory (see Supplementary Note 6).

Onceidentified,weremovethecontributionoftheHO• tothespindensitybyexclud- ingasphericalregionofradius1.5˚AcIenteredonitsoxygenatom(seeSupplementary Note 7). This step yields a modified cube file that provides a more meaningful repre- sentationofthespindensitytolocalizethehydratedexcesselectron.Wethencompute the electron center and gyraLtion radius from this file using the protocol described in previous studies[45] (seCe Supplementary Note 8). Since the cube files were saved every 5 fs and at the crossing geometry, the electron center and gyration radius for the intermediate fraTmes were obtained using linear interpolation between consecutive cube files.

Once all the species were correctly identified, we classified each trajectory accord- ingtoitsreactionmechanism.Forthis,weanalyzedthefinalframeofeachsimulation (corresponding to the geometry at the point of non-radiative crossing) and we com- puted the distance between the electron center and the nearest hydrogen atom. If the distance was less than 0.6 ˚A, the trajectory was classified as hydrogen atom transfer (HAT),otherwiseitwasassignedtotheprotoncoupledelectrontransfer(PCET).The 0.6 ˚A threshold was chosen based on the radial distribution function (RDF) between the electron and all hydrogen atoms (see Supplementary Figure 9). Specifically, it represents the midpoint between two cases: a bound hydrogen atom (distance elec- tron and nearest H 0 ˚A) and the first maximum associated with hydrogen atoms in the solvation shell of the electron (distance between electron and nearest H 1.2 ˚A)[45, 51]. A short electron–hydrogen distance indicates that both species occupy the samespatialregion,consistentwiththeformationofaneutralhydrogenatom.Incon- trast, larger distances suggest that the electron remains solvated while the proton is bound to nearby water molecules, consistent with a PCET process.

Based on the preceding analysis, each trajectory can be unambiguously assigned to a specific mechanism. This allows us to determine the quantum yield and lifetime associated with each mechanism using the following expression:

QY = 100 * (#Trajectories with M / #Total trajectories)
τ = (NM tM) / #Total trajectories with M
where M represent the mechanism (HAT or PCET), QY,τ are the quantum yields and lifetimes, N is the number of trajectories with the M mechanism and tM is the time at which the trajectory i with the M mechanism finds the ”crossing point”.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In the manuscript Ultrafast Solvent Dynamics Drives the Formation of the Hydrated Electron in Photoexcited Water, the authors reported a study on the excited dynamics of photoexcited water, which is important in understanding the formation of the hydrated electron. By running excited-state molecular dynamics, the manuscript discovered two possible pathways after the liquid water becomes photoexcited: hydrogen atom transfer (HAT), which does not form the hydrated electron in the excited state, and proton-coupled electron transfer (PCET), which undergoes localization of the hydrated electron. The authors performed a detailed analysis of their results and provided solid evidence supporting their findings. I recommend it to be published after minor revisions, which may address my following questions or concerns:\n1. The S1 excited state is simulated using the ROKS scheme, as in the authors’ previous CTTS study [J. Phys. Chem. Lett. 2023, 14 (4), 870–878]. Similar CTTS simulations have also been reported in [Nat. Commun. 15, 2544 (2024)]. The S1 state could alternatively be simulated using LR-TDDFT. Have the authors attempted to perform the dynamics using LR-TDDFT for comparison?\n2. In line 104 – 105, the authors stated that their prediction that the HAT mechanism exhibits a higher quantum yield is consistent with experimental observations. Also, in line 106 – 107, the authors mentioned that the hydronium radical has not been experimentally confirmed. Then how is the quantum yield of the HAT mechanism determined?\n3. In line 203, the authors stated that 2 out of 100 trajectories remain in the S1 excited state. It seems confusing to me. As I understand, the system always stays on the S1 excited state, and one tracks the energy gap between S1 and S0. The system can be the S1 to S0 crossing in most trajectories, but this work does not consider the hopping from S1 to S0. If this is correct, every trajectory should remain in the S1 state. Do the authors mean that the 2 out of 100 trajectories do not reach S1 to S0 crossing here?\n4. In Section 2.3, the authors showed different geometries of the HAT products, i.e., Figures 3B and 3D. Are they both observed as the final state at the S1 to S0 crossing in the MD trajectories? If so, what is the ratio between these two products? Can the authors use their simulation results to discuss a bit on which of these two geometries is favored? Moreover, in the HAT pathway, will the hydrated electron form after the non-radiative decay to the ground state?\n5. In line 254 – 257, the authors mentioned that the unstable geometry, as shown in Figure 3D, ultimately leads to non-radiative decay to the ground state. Does this imply that the non-radiative decay is only associated to the unstable geometry in Figure 3D? Does the geometry where the hydrogen atom enters an empty cavity, i.e., Figure 3B, also undergo the non-radiative decay? From Figure 2A, it seems that cavity geometries might also decay, but the authors connected the decay with the unstable geometry, which is confusing.\n6. For consistency, the caption of Figure 4 might start with “Proton Coupled Electron Transfer Mechanism”. The “Translation” in the title of Section 2.5 might be changed to “Translational”", "reviewer_2": "The manuscript \"Ultrafast Solvent Dynamics Drives the Formation of the Hydrated Electron in Photoexcited Water\" by G. Díaz Mirón et al is devoted to simulating excited state dynamics of liquid water in order to reveal initial stages of solvated electron's formation. By ab initio molecular dynamics and careful data analysis the authors characterize two main processes after photoexcitations: hydrogen atom abstraction (HAB) and proton-coupled electron transfer. HAB involves H3O radical as a transient species, never characterized before, which comprises considerable novelty. PCET is found to lead to the formation of the solvated electron and is driven by the solved dynamics. The research is carefully designed, thoroughly executed and clearly reported. It provides considerable novelty already at the stage of asking a question - how is solvated electron formed by light? - which has not been addressed by previous studies. The answer appears to be convincing, too. At the same time I have several serious methodological concerns, which need to be addressed:\n1) My gut feeling and personal experience tell that early stages of the solvated electron's formation has not been considered by theoreticians not out of curiosity shortage, but rather due to the lack of proper methodologies. ROKS with hybrid DFT used here is a cheap and efficient solution, which undeserved remained in shadows. It exhibits excellent performance for excitation energies of liquid water as shown here. Nevertheless, its ability to describe bond breaking is far less obvious. On the contrary, spin-restricted theories (namely ROHF, a parent method of ROKS) is notorious for pathological behavior for bond-breaking. Instead, one uses spin-unrestricted variants of SCF, poor-man's strong correlation for chemical reactions. ROKS is not exactly ROHF, but I have never seen any unequivocal evidence for its performance for dissociation. The authors state that this is demonstrated in the references, however, I could find that only ref 50 has to do with hydrogen dissociation, and it does not look conclusive as benchmarking ROKS was not the focus. Therefore, it would be important to see at least energy profiles of hydrogen abstraction from a water molecule or a water dimer to be able to assess the performance of ROKS.\n2) \"Later\" stages of the trajectories considered here show that S0 and S1 states become close in energy for both processes, these stretches of dynamics being relatively long for PCET. It implies that non-adiabatic effects must play considerable role here. I won't suggest to remake everything with surface-hopping approaches, but it's an important issue. Actually computation of non-adiabatic couplings in CP2K is now possible due to the work of A. Hehn et al.\n3) Spin densities in ROKS are not clearly defined: CP2K writes .cube files but issues a warning in the output. And indeed, the total multiplicity is singlet and strictly speaking there is no spin density despite the presence of unpaired electrons (AKA singlet diradicals). Please, clarify the point.\nA minor point:\nIt may be instructive to compare the role of solvent dynamics in solvated electron's formation with the role of solvent dynamics in its reaction with CO2. It appears that there may be useful analogies and generatizaiton. See:\nhttps://pubs.acs.org/doi/abs/10.1021/acs.jpcb.0c07859\nand\nhttps://pubs.acs.org/doi/10.1021/acs.jpcb.3c06935", "reviewer_3": "The ms. by Miron et al. deals with ultrafast dynamics following photoexcitation in liquid water, specifically focusing on the formation of the hydrated electron. Using Restricted Open-Shell Kohn Sham (ROKS) molecular dynamics simulations, the authors aim to understand the initial steps following photoexcitation in water. They explore in detail the initial excitation localization and the subsequent hydrogen atom transfer and proton-coupled electron transfer pathways, and how this is related to formation and relaxation of the hydrated electron. In this way, the authors paint a detailed picture of how photoexcitation of water initiates a cascade of events, ultimately leading to the formation of the hydrated electron through a specific interplay of proton and electron transfer, coupled with the dynamic restructuring of the hydrogen-bond network. Moreover, they manage to successfully link their results from their DFT dynamics to experimental observations, wherever available.\nThis is a well-executed and insightful study that provides a valuable contribution to our understanding of the photodynamics of liquid water, that should be interesting to a broad readership. Therefore, the ms. should in my view be published after a minor revision that reflects the issues raised below.\n1. The accuracy of the ROKS results depends heavily on the functional used, basis set, and how well the method captures dynamic correlation. A more thorough assessment of the performance for excited-state dynamics in photoionized water, would be valuable.\n2. I would appreciate a more detailed description how exactly the ROKS method is used for excited state dynamics. Comparison to other methods such as TDDFT would be valuable.\n3. The authors acknowledge the simplification of not explicitly including non-adiabatic couplings. While it may be reasonable for the present system to consider the energy gaps only, a more rigorous justification or exploration of the potential impact of the approximation used on the results would strengthen the argument.\n4. Given the computational demands of the method, it is understandable that the authors used a relatively small unit cell with limited sampling of initial configurations. Nevertheless, at least some discussion on the size and sampling effects would be appropriate.\n5. I am puzzled by the author’s finding of the H3O radical. This must be a very unstable species (and the authors indeed discuss its transient character) but I wonder if it is not more like a H3O+…e- contact “ion pair”.", "reviewer_4": "The authors describe a simulation study for studying the formation of the hydrated electron. This process is a benchmark chemical process and has been studied experimentally and theoretically by many groups. Here the authors present a study which follows the pathway from the photoexcitation into the first absorption band all the way.\nThey describe two distinct pathways: Hydrogen Atom Transfer (HAT) and Proton Coupled Electron Transfer (PCET). The predictions for the HAT mechanism are in excellent agreement with the experimenal time constants.\nThere data analysis involves quantification of the hydrogen bond network defect, the localization the solvated electron.\nInvestigations of the H Atom Transfer and predictions of the photo emission properties of the hydrated electron.\nI have the following comments:\nTheir key analysis tool is to use “dilution factor” or “inverse participation ratio” to quantify how many water molecules an excess electron is delocalized over.\nNote: strictly speaking the “PDFs” in Fig. 1BC are counts, not probability distribution functions; they should be normalized, or maybe called probability counts.\nThe conclusion is almost complete localization: the peak in Fig. 1B for 1 water is 50, at 2 water molecules we are already down to ~2 counts, or 4%. This need to be mentioned more clearly So, the further analysis does not change the view compare to with one water. The figures in Fig.2 B could clearer, they are hard to see on the screen.\nDefects account for ¾ of the binding sites, mainly missing acceptor sites: a missing acceptor bond means the water molecule has a lone pair on oxygen not bonded to an H from another water molecule, destabilizing the n to sigma* transition to ~8 eV. What is missing here is a further aspect: the steric room available next to an unsatisfied acceptor bond also helps the HO + H dissociation process by making room for the recoiling OH and H.\nThe definition of the decay lifetime on p. 9 is somewhat arbitrary but consistent, sufficient tbuto distinguish processes differing by more than a factor of 10 in lifetime. May be this definition can be made more precise.\nIn the PCET mechanism (n→sigma* H2O + H2O→ HO + e- + H+ +H2O→ HO + e- + H3O+) , the e- quickly becomes localized in a trapped state. The OH and H3O+ become quite well separated (6 Å in Fig. 5D).\nTheir theoretical study reveals three distinct minima in the HO* and (H3O+) distance which are in nice agreement with the study of the recent electron diffraction experiments upon photoionization. These could be rationalized by the formation of contact ion radical pair and the separation mediated by the solvent.\nThe results are not surprising, the new aspect is that the entire process can be modelled from photoexcitation up to the localized electron directly capturing the solvent dynamics.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Comments for: Ultrafast Solvent\nDynamics Drives the Formation of the\nHydrated Electron in Photoexcited Water\nNovember 2025\nDear Reviewers,\nWe sincerely appreciate the time and effort that the Reviewers have dedicated to\nevaluating our manuscript and providing thoughtful and constructive feedback. One\nof the main concerns raised pertains to the validation of our computational approach\nspecifically the use of ROKS to study the excited-state dynamics of liquid water.\nWhiletheoriginalversionofthemanuscriptalreadyincludedavalidationagainstTD-\nDFT results, we fully agree with the editors and reviewers that this aspect required\na much more extensive and detailed treatment.\nBefore providing detailed responses to each of the reviewers’ comments, we would like\nto briefly summarize the main changes that have been implemented in the revised\nmanuscript. This has involved conducting significant benchmarking and validation\nwith new calculations including:\n• An assessment of the robustness of the ROKS approach using different basis\nsets and exchange-correlation functionals.\n• Validation of ROKS against TDDFT for both the HAT and PCET mechanisms.\n• To gain deeper insights, we additionally benchmark ROKS and TDDFT results\nagainst multiconfigurational CASPT2 calculations for a hydrogen-bonded water\ndimer in vacuum.\n• Finally,wealsodeterminetheLandau–Zener(LZ)transitionprobabilities,which\nincorporate not only the energy gap but also the curvature of the potential en-\nergy surfaces.\nto the revised manuscript (pages 23-26), which incorporates most of the modifications\ndiscussed above.\nistry of Liquid Water\nIn this work, we employed the ROKS approach to investigate the photochemistry\nof liquid water. While our results show excellent qualitative and quantitative\nagreement with a broad range of experimental observations, it is important to\nassess the quality of the electronic structure method and to evaluate its perfor-\nmance relative to more extensively tested theoretical approaches in photochem-\nistry.\nIn order to first assess the robustness of ROKS with respect to the choice of\nexchange–correlation functional and basis set, we performed additional calcu-\nlations using the CAM-B3LYP functional[1] and the TZV2P basis set[2]. The\nresults, presented in section S9 in the Supplementary Information, demonstrate\nthat the ROKS method provides consistent energy gaps and reproduces the same\nbehavior for both photochemical pathways (HAT and PCET), confirming its re-\nliability across different electronic structure settings.\nWe further compared our results with the TDDFT method, which has been widely\napplied to the study of excited states and in non-adiabatic dynamics. In addition,\nfor a hydrogen-bonded water dimer in vacuum, we benchmarked our approach\nagainst both multiconfigurational methods and TDDFT. The main results are\nsummarized in Figure 7.\nwater. Panel A: Evolution of the energy gap for a representative HAT trajec-\ntory obtained with ROKS (green) and TDDFT (purple). Panel B: Evolution of\nthe energy gap for a representative PCET trajectory obtained with ROKS (green)\nand TDDFT (purple). Panel C: Spin density from ROKS and LUMO orbital\nfrom TDDFT at 250 fs of the PCET trajectory, illustrating the localization of\nthe hydrated electron. Panel D: Minimum Energy Path (MEP) for a water\ndimer in vacuum exhibiting the HAT mechanism, obtained with CASPT2. The\ninset shows the geometry at the conical intersection between S (purple) and S\n1 0\n(red). Data reproduced from ref. 71. Panel E: Linear Interpolated Path (LIP)\nfor the same water dimer in vacuum exhibiting the HAT mechanism, obtained\nwith ROKS-CAM-B3LYP. The inset shows the geometry at the S /S conical\n1 0\nintersection. Panel F: Linear Interpolated Path (LIP) for the water dimer ex-\nhibiting the HAT mechanism, obtained with TDDFT-CAM-B3LYP. The inset\nshows the geometry at the S /S conical intersection.\n1 0\nstate (S ) and the ground state (S ) for two representative trajectories (sampled\n1 0\nfrom the ROKS dynamics), corresponding to the HAT and PCET mechanisms,\nobtained using both ROKS and TDDFT with the same functional and basis\nset employed in this work (see Computational Methods section). At the initial\nstep, both methods yield nearly identical energy gaps, as further confirmed by a\nbroader validation across all initial conditions (see Figure S1 in the Supplemen-\ntary Information).\nThe main difference arises in the HAT mechanism. At the end of the simula-\ntion, ROKS predicts a crossing point with an energy gap smaller than 0.2 eV,\nwhereas TDDFT yields a value near 1.5 eV for the same geometry. Despite\nthis quantitative difference, both methods describe the same physical process:\nwater dissociation leading to the formation of HO• and H•. In contrast, the\nPCET mechanism shows excellent agreement between ROKS and TDDFT. In\nthis case, both methods predict the dissociation of a water molecule into HO•\nand H O+ within the first 20 fs, with comparable energy gap profiles. Moreover,\nthe ROKS results, as shown in Figure 7C, which compares the ROKS spin den-\nsity and TDDFT LUMO at 250 fs.\nTo further understand the discrepancy between ROKS and TDDFT observed\nfor the HAT mechanism, we examined a model system consisting of a hydrogen-\nbonded water dimer in vacuum [3, 4]. Segarra et al. computed the Minimum\nEnergy Path (MEP) for the photodissociation of one water molecule from the\nFranck–Condon region to the S /S conical intersection using the multiconfig-\n1 0\nurational CASPT2 method (see Figure 7D). Using both the ground-state and\nconical intersection optimized geometries at this high level of theory, we com-\nputed the potential energy surfaces along the Linear Interpolated Path (LIP)\nwith ROKS (Figure 7E) and TDDFT (Figure 7F) employing the CAM-B3LYP\nfunctional.\na HO•–H• distance of approximately 2 ˚A for ROKS, compared with 2.8 ˚A for\nCASPT2. Although this difference is not insignificant, a 2 ˚A separation already\nindicates a fully dissociated configuration, confirming that ROKS and CASPT2\ncapture the same underlying photodissociation physics. In contrast, TDDFT\npredicts a local minimum on the S surface at around 2 ˚A with an energy gap\n1\nof ∼ 1.5 eV, which is consistent with the differences seen in the condensed\nphase simulations for the HAT mechanism (Figure 7A). Moreover, TDDFT\nexhibits a ∼ 1 eV barrier from this minimum to the conical intersection, which\ncould significantly delay, or even prevent, reaching the intersection within the\ntimescales (∼ 10 fs) reported in previous studies71,72.\nTaken all together, these results provide a comprehensive validation of the ROKS\napproach. We demonstrate that ROKS captures the essential physics of both\nthe HAT and PCET photochemical pathways in water, in close agreement with\nTDDFT and CASPT2 benchmarks. This consistency across different mecha-\nnisms and details of the underlying electronic structure supports the robustness\nof our methodology and reinforces confidence in its predictive capability for mod-\neling the excited-state dynamics of liquid water.\nWe believe that these substantial additions significantly enhance both the rigor and\nclarity of the manuscript, further strengthening the reliability and impact of the\nresults presented in this work. We also would like to clarify that we have decided\nto change the title of the manuscript to one more accurately reflects the scope and\nfindings of the present study. The revised title is: The Photochemical Birth of the\nHydrated Electron in Liquid Water. In the following sections, we provide detailed\nresponses to each point raised by the Reviewers. The corresponding revisions are\nindicatedinthefile“manuscript marked”,wherepageandlinenumbersarereferenced\nfor clarity.\n1.1: The S1 excited state is simulated using the ROKS scheme, as in\nthe authors’ previous CTTS study [J. Phys. Chem. Lett. 2023, 14 (4),\n870–878]. Similar CTTS simulations have also been reported in [Nat.\nCommun. 15, 2544 (2024)]. The S1 state could alternatively be simulated\nusing LR-TDDFT. Have the authors attempted to perform the dynamics\nusing LR-TDDFT for comparison?\nWe thank the Reviewer for this valuable comment. We fully agree that the S state\n1\ncould, in principle, also be simulated using LR-TDDFT. Our decision to employ the\nROKS approach was motivated by two main factors. Firstly, for the most part, until\nonlyveryrecently, LR-TDDFTforexcited-statedynamicshasbeenprimarilyusedfor\nclusters in the gas phase or embedded within a QM/MM setup. LR-TDDFT calcu-\nlations in fully periodic systems have certainly been deployed to compute absorption\nspectra[5]. ThustheuseofLR-TDDFTforexcited-statedynamicsinperiodicsystems\nis far from mainstream apart from some very recent exceptions[6]. Secondly, ROKS\nprovides a computationally feasible alternative to LR-TDDFT and as as evidenced\nin a previous CTTS study published in JPCL, it is a reasonable compromise. Our\nbenchmarks and validation further show that ROKS is perfectly suitable for studying\nthe excited-state dynamics of liquid water.\nWe also appreciate the Reviewer for bringing to our attention the recent Nature Com-\nmunications paper, which we were unaware of at the time of the original submission.\nWe have now cited this work in the revised manuscript, as it further supports and\nvalidates the methodology adopted in our study. The following modification has been\nmade accordingly (page 30, lines 570-573):\nExcellent agreement with the experiments for aqueous iodine solutions were re-\nported by Carter et. al.51 within ROKS framework and by Lan et. al.111 em-\nploying a similar approach with the same functional.\nIn the original version of our manuscript, we had already included a preliminary val-\nidation of ROKS against TDDFT, such as the absorption spectra computed for all\n for all\ninitial conditions and the evolution of the energy gap along a representative PCET\ntrajectory (see Figures S1 and 6B in the previous version, respectively). However,\nfollowing the Reviewers’ suggestion for a more detailed validation, we have now sub-\nstantially expanded this part of the study.\nSpecifically, wehaveaddedanewsectiontitled“Evaluating the Accuracy of ROKS for\nthe Photochemistry of Liquid Water.”(pages 23-26 in the revised manuscript). In this\n6\nsection, we provide a comprehensive comparison between ROKS and LR-TDDFT.\nWe include energy gap profiles for two representative trajectories (one for the HAT\nmechanism and one for the PCET mechanism) calculated using both methods. To\nfurther clarify the observed differences, we also performed benchmark calculations for\na hydrogen-bonded water dimer in vacuum, comparing ROKS and TDDFT results\nwith high-level multiconfigurational CASPT2 calculations. Since we show many ex-\nplicit benchmarks for ROKS in the rest of the responses, the reviewer is encouraged\nto also examine the responses there.\nThese additional results provide a thorough validation of the ROKS approach and\nreinforce confidence in its ability to accurately describe the photochemistry of liquid\nwater.\n1.2: In line 104 – 105, the authors stated that their prediction that the\nHAT mechanism exhibits a higher quantum yield is consistent with exper-\nimental observations. Also, in line 106 – 107, the authors mentioned that\nthe hydronium radical has not been experimentally confirmed. Then how\nis the quantum yield of the HAT mechanism determined?\nWe thanktheReviewerforraisingthis important point. Weagreethatthedistinction\nbetween the experimentally detected species and those predicted in our simulations\nneeds clarification.\nThe hydronium radical (H O•) is indeed part of our HAT mechanism. In this process,\n3\nthe excess electron remains localized near a hydrogen species, which can exist either\nas a free H atom (H•) or transiently as part of the hydronium radical (H O•). As\n3\nexplained in our manuscript, both cases correspond to the same mechanistic class,\nbecause the electron-hydrogen distance identifies an overall H-centered configuration\n(see Analysis Protocols). Experimentally, only the free H atom has been detected;\nhowever, our simulations suggest the formation of a short-lived intermediate, the\nH O• radical, which to our knowledge has not yet been confirmed experimentally.\nphasize that in our computational framework the quantum yield is defined based on\nthe number of trajectories that follow a given mechanistic pathway. Because each\ntrajectory can be unambiguously assigned to either HAT or PCET, the quantum\nyield simply reflects the statistical occurrence of each mechanism within our ensem-\nble. Additional information is provided for the different situations founded in HAT\nmechanism in our response to point 1.4, which closely relates to this question.\nTo make this clearer, we have revised the Analysis Protocols section to explicitly\ndefine the quantum yields and lifetimes (paged 35-36, lines 680–688) as follows:\n7\n... Based on the preceding analysis, each trajectory can be unambiguously as-\nsigned to a specific mechanism. This allows us to determine the quantum yield\nand lifetime associated with each mechanism using the following expression:\n#Trajectories with M\nQY = ∗100 (6)\nM\n#Total trajectories\n(cid:80)NM tM\nτ = i i (7)\nM\n#Total trajectories with M\nwhere M represents the mechanism (HAT or PCET), QY,τ are the quantum\nyields and lifetimes respectively, N is the number of trajectories with the M\nM\nmechanism and tM is the time at which the trajectory i with the M mechanism\ni\nfinds the ”crossing point”.\n1.3: In line 203, the authors stated that 2 out of 100 trajectories remain\nin the S1 excited state. It seems confusing to me. As I understand, the\nsystem always stays on the S1 excited state, and one tracks the energy\ngap between S1 and S0. The system can be the S1 to S0 crossing in most\ntrajectories, but this work does not consider the hopping from S1 to S0.\nIf this is correct, every trajectory should remain in the S1 state. Do the\nauthors mean that the 2 out of 100 trajectories do not reach S1 to S0\ncrossing here?\nWe thank the Reviewer for pointing out this ambiguity in our manuscript. The\nanalysis of the Reviewer is correct, the statement “2 out of 100 trajectories remain\nin the S excited state” refers to those trajectories in which the system does not\n1\nreach the crossing condition within the 1 ps simulation window. In other words,\nthese trajectories remain on the S potential surface for the entire duration of the\n1\nsimulation.\nToeliminatethisconfusion,wehaverevisedthecorrespondingsentenceinthemanuscript\n(page 10, lines 205–207) as follows:\nThe PCET mechanism displays a slower decay, on the order of hundreds of\nfemtoseconds (see Figure 2D). In fact, 2 out of 100 trajectories remain in the\nS excited state for the entire 1 ps simulation do not exhibit non-radiative decay\n1\nto the S state, remaining on the S for the entire 1 ps of the simulation.\n0 1\n1.4: In Section 2.3, the authors showed different geometries of the HAT\n8\nproducts, i.e., Figures 3B and 3D. Are they both observed as the final state\nat the S1 to S0 crossing in the MD trajectories? If so, what is the ratio\nbetween these two products? Can the authors use their simulation results\nto discuss a bit on which of these two geometries is favored? Moreover, in\nthe HAT pathway, will the hydrated electron form after the non-radiative\ndecay to the ground state?\nWe thank the Reviewer for raising this insightful question. As correctly pointed out,\nboth conformations 3B and 3D (and in a few cases 3C) were indeed observed at the\nnon-radiative S → S crossing points in our simulations. In Figure 3, all red points\n1 0\nrepresent the geometries corresponding to these crossing events. We have clarified\nthis explicitly in the revised manuscript (page 13, line 246):\n...leading to a non-radiative transition (as indicated by all the red points in the\nFigure 3)\nRegarding the second part of the Reviewer’s question, we fully agree that a discussion\ncomparing the different pathways within the HAT mechanism is valuable. Accord-\ningly, we have expanded this section in the revised manuscript to include an analysis\nof the lifetimes and quantum yields associated with each photoproduct. The following\nchanges were introduced in Section 2.3 (pages 14–15, lines 264-274):\n... An interesting observation is that approximately 20% of the trajectories follow\nthe mechanism involving the formation of a hydrogen atom in an empty cavity\n(Figure 3B) with a lifetime of 13 fs, while the remaining 33% correspond to the\nalternative HAT pathway (Figures 3C and 3D) exhibiting a longer lifetime of\n34 fs. These results correlate well with the defect analysis presented in Figure\n1C. When the dissociating water molecule has a missing hydrogen-bond donor\n(AD or AAD defect), there is sufficient space available to allow an ultrafast\nhydrogen atom dissociation into a cavity (∼13 fs), in excellent agreement with\nvalues reported for a hydrogen-bonded water dimer in vacuum71,72. In contrast,\nwhen the dissociating molecule lacks a hydrogen-bond acceptor (ADD defect),\nthe surrounding hydrogen-bond network in the condensed phase restricts this\nmotion, delaying the non-radiative decay and promoting the formation of the\ntransient H O• species.\nnism the hydrated electron could, in principle, form after the non-radiative decay to\nthe ground state. However, our present simulations are restricted to the excited-state\ndynamics and do not explicitly include the possibility for the electron ejection process\n9\nto occur in the ground state. We already included this limitation in the Conclusions\nsection (page 27, lines 488–490):\nWhileoursimulations, performedintheelectronicexcitedstate, cannotexclude\nthe possibility that it could release an electron in the ground electronic state,\nfurther studies are needed to confirm this pathway.\n1.5: In line 254 – 257, the authors mentioned that the unstable geome-\ntry, as shown in Figure 3D, ultimately leads to non-radiative decay to the\nground state. Does this imply that the non-radiative decay is only asso-\nciated to the unstable geometry in Figure 3D? Does the geometry where\nthe hydrogen atom enters an empty cavity, i.e., Figure 3B, also undergo\nthe non-radiative decay? From Figure 2A, it seems that cavity geometries\nmight also decay, but the authors connected the decay with the unstable\ngeometry, which is confusing.\nWe thank the Reviewer for bringing this to our attention. We acknowledge that our\n that our\noriginal description may not have been sufficiently clear, which could have led to\nconfusion. The Reviewer is indeed correct that both conformations 3B and 3D lead\nto non-radiative decay. We have revised the manuscript accordingly to clarify this\npoint.\nAs mentioned in our response to point 1.4, we first clarified the meaning of all red\npoints in Figure 3 (page 13, line 246).\n...leading to a non-radiative transition (as denoted by all the red points in the\nFigure 3)\nIn addition, we explicitly state that conformation 3B also leads to non-radiative decay\n(page 13, lines 248–249).\n... In some trajectories, the dissociating hydrogen atom moves into an empty\ncavity(Figure3B),whereitsnearestoxygenispartofahydroxylradical(HO•),\nresultinginCN = 1and leading to a non-radiative crossing to the ground state.\nO\n1.6: For consistency, the caption of Figure 4 might start with “Proton\nCoupled Electron Transfer Mechanism”. The “Translation” in the title of\nSection 2.5 might be changed to “Translational”\nWe thank the Reviewer for carefully pointing out these errors. Both issues have been\ncorrected in the revised version of the manuscript.\n2.1: My gut feeling and personal experience tell that early stages of the\nsolvated electron’s formation has not been considered by theoreticians not\nout of curiosity shortage, but rather due to the lack of proper methodolo-\ngies. ROKS with hybrid DFT used here is a cheap and efficient solution,\nwhich undeserved remained in shadows. It exhibits excellent performance\nfor excitation energies of liquid water as shown here. Nevertheless, its\nability to describe bond breaking is far less obvious. On the contrary,\nspin-restricted theories (namely ROHF, a parent method of ROKS) is\nnotorious for pathological behavior for bond-breaking. Instead, one uses\nspin-unrestricted variants of SCF, poor-man’s strong correlation for chem-\nical reactions. ROKS is not exactly ROHF, but I have never seen any\nunequivocal evidence for its performance for dissociation. The authors\nstate that this is demonstrated in the references, however, I could find\nthat only ref 50 has to do with hydrogen dissociation, and it does not look\nconclusive as benchmarking ROKS was not the focus. Therefore, it would\nbe important to see at least energy profiles of hydrogen abstraction from\na water molecule or a water dimer to be able to assess the performance of\nROKS.\nWe thank the Reviewer for highlighting this important aspect regarding the valida-\ntion of our ROKS approach. We fully agree that a proper assessment of its accuracy\nis essential. For this reason, we have included in the revised manuscript a new section\ntitled “Evaluating the Accuracy of ROKS for the Photochemistry of Liquid Water”\n(pages 23-26 in the revised version and also reproduced at the beginning of this re-\nsponse letter). In this section, we validate our ROKS results against TDDFT for\nboth mechanisms identified in the photochemistry of liquid water (HAT and PCET),\nwhere in both scenarios, we have a bond breaking of water molecule. To gain deeper\ninsight into the performance of ROKS, we additionally tested the method using a\nhydrogen-bonded water dimer in vacuum, a model that exhibits only the HAT mech-\nanism. We then compared the ROKS and TDDFT results with multiconfigurational\nCASPT2calculationsthathavebeenpreviouslyreported[3]. Althoughsomequantita-\ntive differences were observed among the methods, ROKS reproduces the dissociative\ncharacter predicted by CASPT2 much more accurately than TDDFT. We also vali-\ndated the ROKS approach using a different exchange–correlation functional and basis\nset. These additional results are provided in Section S9 Validation of functional and\nbasis set in ROKS approach of the Supplementary Information.\nOverall, these results confirm that ROKS provides a reliable and computationally\nefficient framework for describing the key excited-state processes involved in the pho-\n11\ntochemistry of liquid water.\n2.2: ”Later” stages of the trajectories considered here show that S0 and\nS1 states become close in energy for both processes, these stretches of dy-\nnamics being relatively long for PCET. It implies that non-adiabatic effects\nmust play considerable role here. I won’t suggest to remake everything\nwith surface-hopping approaches, but it’s an important issue. Actually\ncomputation of non-adiabatic couplings in CP2K is now possible due to\nthe work of A. Hehn et al.\nWe sincerely thank the Reviewer for highlighting this excellent and important point.\nWe fully agree that non-adiabatic effects can play a significant role in the photo-\nchemistry of water. However, several considerations must be clarified to explain why\nwe did not explicitly include non-adiabatic coupling vectors (NACVs) in our present\nsimulations.\nFirst, our study employs the ROKS approach, chosen based on previous works show-\ning that the hydrated electron forms upon excitation through charge transfer to the\nsolvent. Intherevisedmanuscript,wehaveaddedadetailedvalidationsectiondemon-\nstrating that ROKS provides a reliable and physically consistent description of this\nprocess (see also point 2.1). Unfortunately, at present, CP2K does not include an\nimplementation for calculating NACVs within the ROKS framework, which prevents\na direct evaluation of these quantities in our simulations.\nAlternatively, we could have performed the simulations using TDDFT. However, it\nis well known that linear-response TDDFT struggles to accurately describe coni-\ncal intersections, often predicting incorrect dimensionality between the excited and\nground electronic states[7, 8]. One of the most robust and widely used strate-\ngies to model non-adiabatic transitions in such cases is to employ an energy gap\ncriterion[7, 8, 9, 10, 11], as we have done in this work. This approach is justified by\nthe mathematical relationship between NACVs and energy gap[12], which shows that\nNACVs increase sharply when the energy gap between states decreases, effectively\ncapturing the region where non-adiabatic transitions are most likely to occur.\nAnother possible approach for modeling nonadiabatic dynamics is to use surface hop-\nping methods based on the Landau–Zener (LZ) or Baeck–An (BA) approaches. These\nalgorithms do not explicitly compute NACVs but incorporate their effects implicitly\nthrough the dependence of transition probabilities on the energy gap and the curva-\nture of the potential energy surfaces involved in the transition.\nWe would like to thank the reviewer for bringing to our attention the recent work of\nA. Hehn et al.[6], which indeed represents an important step toward including non-\nadiabatic effects with CP2K. However, NACVs between the S and S states were not\n1 0\n12\ncomputed directly in that work for the same reasons discussed above in the context\nof TDDFT. Instead, the authors introduced the transition to the ground state via\nthe BA algorithm. Their results showed that the BA method tends to overestimate\nthe transition probability to S (see Figure 6 in Ref.[6]), highlighting the current\n0\nlimitations of this approach.\nConsidering these challenges, we hope the reviewer understands the practical and\nmethodological reasons that prevented us from including explicit NACVs in our\npresent study. We also hope that our work motivates further development of NACV\ncalculations within the ROKS framework, which would enable more accurate non-\nadiabatic simulations in the future.\nFinally, to provide a more quantitative connection with nonadiabatic dynamics, we\nhave now computed the LZ transition probabilities (P ), which have been shown to\nLZ\nreproduce the results of traditional surface hopping algorithms that explicitly include\nNACVs across a variety of systems and electronic structure methods[13, 14, 15, 16].\nThe new results are presented in the Section S10 Incorporating Non-Adiabatic Effects\ninto the Photochemistry of Liquid Water in the Supplementary Information (here\nwe show the main results). As shown in the Figure below, the P values increase\nLZ\nsignificantly near the S /S crossing points for representative trajectories of both\n1 0\nmechanisms, supporting our choice of using the energy gap criterion to identify non-\nradiative transitions.\nFigure S11: Landau-Zener (LZ) Probabilities. Energy gap (purple line) and LZ\nprobabilities (red line) as a function of time for representative HAT (Panel A)\nand PCET (Panel B) trajectories. Panel C: shows a 2D density plot of the\nLZ probabilities versus energy gap for all trajectories, with the inset providing\na zoomed-in view at small LZ values. Red points indicate the values at the\nS → S crossing events.\n1 0\nWe also add the following in the revised manuscript in the sections (page 31, lines\n13\n590-594):\n..., justifying the omission of explicitly determining NACVs in our\nsimulations113,114. To further support this approximation, we evaluated the\nLandau-Zener (LZ) transition probabilities115, which in addition to the energy\ngap criterion, also account for the curvature of the potential energy surfaces (see\nsection S10 in the Supplementary Information). This approach has been shown\nto reproduce the results consistent with those obtained using explicit NACVs in\nnon-adiabatic dynamics115−119\n2.3: Spin densities in ROKS are not clearly defined: CP2K writes .cube\nfiles but issues a warning in the output. And indeed, the total multiplicity\nis singlet and strictly speaking there is no spin density despite the presence\nof unpaired electrons (AKA singlet diradicals). Please, clarify the point.\nWe thank the reviewer for raising this clarification. Indeed, in the ROKS formalism\nthe total spin multiplicity is singlet, and therefore the spin density is not rigorously\ndefined. The “spin density” cube files in this case represents the difference in electron\ndistribution between the two open-shell orbitals involved in the excitation. We use\nthese files to perform all the analysis in the present work.\nNo warnings were printed in the output file, here we show the lines printed in the\noutput in this regard, we also provide inputs to reproduce these results.\nThe sum of alpha and beta density is written in cube file format to the file:\ncp2k-roks-dens-ELECTRON DENSITY-1 0.cube\nThe spin density is written in cube file format to the file:\ncp2k-roks-dens-SPIN DENSITY-1 0.cube\n2.4: It may be instructive to compare the role of solvent dynamics in sol-\nvated electron’s formation with the role of solvent dynamics in its reaction\nwith CO2. It appears that there may be useful analogies and generati-\nzaiton. See:\nhttps://pubs.acs.org/doi/abs/10.1021/acs.jpcb.0c07859 and\nhttps://pubs.acs.org/doi/10.1021/acs.jpcb.3c06935\nWe thank the reviewer for bringing these papers to our attention. We have now\nadded a discussion in the section where we talk about the coupled rotational and\n14\ntranslational solvent motions that stabilize the hydrated electron (page 18, lines 340-\n344).\n...consistentwithexperimentalfindingsfromionizationstudiesinliquidwater35.\nA similar reorganizational response has been reported in ab initio molecular dy-\nnamics simulations of carbon dioxide reduction in water77,78, where the breaking\nof hydrogen bonds in the second solvation shell triggers the stabilization of the\nexcess electron on the carbon dioxide molecule. This process also involves cou-\npled translational and rotational motions of water molecules, closely paralleling\nthe behavior observed in our simulations.\n3.1: The accuracy of the ROKS results depends heavily on the functional\nused, basis set, and how well the method captures dynamic correlation. A\nmore thorough assessment of the performance for excited-state dynamics\nin photoionized water, would be valuable.\nWe thank the Reviewer for highlighting this important point, which was indeed not\nsufficiently addressed in the original version of our manuscript. In the revised ver-\nsion, we have now included a comprehensive assessment of the reliability of the ROKS\napproach. In particular, we examined the influence of the exchange-correlation func-\ntional and basis set by performing additional calculations using the CAM-B3LYP\nfunctional and the TZV2P basis set. The results, now included in Section S9 Valida-\ntion of functional and basis set in ROKS approach in the Supplementary Information,\nconfirmthattheROKSmethodyieldsconsistentenergygapsandreproducesthesame\nphotochemical behavior for both mechanisms (HAT and PCET), see Figure below.\nROKS framework. Comparison between the DZVP (green line) and TZV2P\n(red line) basis sets using the PBEh(40)-rVV10 functional (PBEh) is shown\nin Panel A and Panel B for the HAT and PCET mechanisms, respectively.\nSimilarly, comparison between the PBEh (green line) and CAM-B3LYP (CAM,\nred line) functionals using the DZVP basis set is shown in Panel C and Panel\nD for the HAT and PCET mechanisms, respectively. For the PCET trajectory,\nthe PBEh-TZV2P (red line in Panel B) and CAM-DZVP (red line, Panel D)\ncalculations were permorfed at configurations sampled every 5 fs to reduce the\ncomputational cost.\nFurthermore, we have added a new section discussing these validation results in detail\nto provide a clearer understanding of the robustness and reliability of the ROKS\napproach (pages 23-26 of the revised manuscript and the beginning of the response\nletter). The following modification has been made in this regard (page 23, lines\n17\n413-418):\n...In order to first assess the robustness of ROKS with respect to the choice of\nexchange–correlation functional and basis set, we performed additional calcu-\nlations using the CAM-B3LYP functional84 and the TZV2P basis set85. The\nresults, presented in section S9 in the Supplementary Information, demonstrate\nthat the ROKS method provides consistent energy gaps and reproduces the same\nbehavior for both photochemical pathways (HAT and PCET), confirming its re-\nliability across different electronic structure settings.\n3.2: I would appreciate a more detailed description how exactly the ROKS\nmethod is used for excited state dynamics. Comparison to other methods\nsuch as TDDFT would be valuable.\nWe thank the Reviewer for highlighting this important point, which was insufficiently\naddressed in the original version of the manuscript. To clarify the theoretical founda-\ntions of the ROKS method, we have added a brief explanation in Section 5.2, along\nwith relevant references for readers interested in further technical details. The corre-\nsponding revisions have been incorporated into the revised manuscript (page 30, lines\n549-557):\nThe excited state molecular dynamics (ESMD) simulations were initiated from\nthe first electronic excited state (S ), computed using the Restricted Open-\n1\nshell Kohn–Sham (ROKS) method46,47,103. ROKS has proven to be a reliable\napproach for describing dissociative and charge-transfer states in closed-shell\nsystems48−50, processes that play a central role in the photochemistry of\nwater. This method provides an alternative approach for calculating singlet\nexcited states in closed-shell systems that are described by non-Aufbau elec-\ntronic configurations. By minimizing the energy expression E = 2E −E ,\nS mix T\nwhere E ,E ,E correspond to the energies of the singlet, triplet, and mixed\nS T mix\nstates, respectively, the method effectively accounts for a two-determinant na-\nture of the excited state. Consequently, ROKS offers a computationally efficient\nand reliable framework for describing dissociative and charge-transfer excited\nstates48−50,104.\nRegarding the validation of our approach against TDDFT, we had already included\nsome comparative calculations in the original version of the manuscript. However,\nsince this concern was raised by all Reviewers, we decided to address it in a more\ncomprehensivemanner. Intherevisedmanuscript, wehaveaddedanewsectiontitled\n“Evaluating the Accuracy of ROKS for the Photochemistry of Liquid Water” (pages\n18\n23-26). In this section, we compare our ROKS results with TDDFT calculations\nusing the same exchange-correlation functional and basis set for two representative\ntrajectories corresponding to each mechanism. To further understand the observed\ndifferences between the two methods, we also performed calculations on a hydrogen-\nbonded water dimer in vacuum and compared the results from ROKS and TDDFT\nwith those obtained from the multiconfigurational CASPT2 method. Although some\ndiscrepancies were found among the three approaches, ROKS reproduces the pho-\ntodissociation behavior of water more accurately and qualitatively aligns better with\nthe CASPT2 reference.\n3.3: The authors acknowledge the simplification of not explicitly includ-\ning non-adiabatic couplings. While it may be reasonable for the present\nsystem to consider the energy gaps only, a more rigorous justification or\nexploration of the potential impact of the approximation used on the re-\nsults would strengthen the argument.\nWe thank the Reviewer for highlighting this important point. As this issue was also\nraised by another Reviewer, to avoid repetition we kindly refer the reviewer to our\nresponse to point 2.2, where we explain in detail why non-adiabatic coupling vectors\n(NACVs) cannot be explicitly incorporated into our current approach.\nNevertheless, we have validated our methodology by calculating Landau–Zener (LZ)\ntransition probabilities, which account not only for the energy gap but also for the\ncurvature of the potential energy surfaces. This approach has been shown to repro-\nduceresultscomparabletothoseobtainedwithtraditionalsurface-hoppingalgorithms\nthatexplicitlyincludeNACVs, acrossawiderangeofsystemsandelectronicstructure\nmethods[13, 14, 15, 16]. The corresponding results are now presented and discussed\nin section S10 in of the Supplementary Information. The main result is shown below.\nprobabilities (red line) as a function of time for representative HAT (Panel A)\nand PCET (Panel B) trajectories. Panel C: shows a 2D density plot of the\nLZ probabilities versus energy gap for all trajectories, with the inset providing\na zoomed-in view at small LZ values. Red points indicate the values at the\nS → S crossing events.\n1 0\nWe have also included in the revised manuscript the following modification (page 31,\nlines 590-594):\n..., justifying the omission of explicitly determining NACVs in our\nsimulations113,114. To further support this approximation, we evaluated the\nLandau-Zener (LZ) transition probabilities115, which in addition to the energy\ngap criterion, also account for the curvature of the potential energy surfaces (see\nsection S10 in the Supplementary Information). This approach has been shown\nto reproduce the results consistent with those obtained using explicit NACVs in\nnon-adiabatic dynamics115−119\nAswecanseeintheFigure, LZprobabilitiesincreasesharplywhentheenergygapbe-\ncomes small, confirming the consistency and validity of our energy gap based criterion\nfor both mechanisms.\n3.4: Given the computational demands of the method, it is understandable\nthat the authors used a relatively small unit cell with limited sampling of\ninitial configurations. Nevertheless, at least some discussion on the size\nand sampling effects would be appropriate.\nWe sincerely thank the Reviewer for raising this important point. The choice of\nsystem size was primarily dictated by the high computational cost of our approach.\nThe present unit cell represents the largest system that can be feasibly simulated\n20\nwithintheROKSframeworkattheabinitiolevelusedinthisstudy. Importantly, this\ncell size is consistent with those typically employed in ab-initio molecular dynamics\nstudies of bulk water, particularly in investigations of the hydrated electron.\nWe acknowledge that finite size effects inevitably impose certain limitations. In par-\nticular, we have already discussed how the restricted box size affects the degree of\nelectron delocalization upon excitation. This limitation was explicitly noted in the\nmanuscript (page 15, lines 285–288):\nExperimentalstudieshavereportedinitialgyrationradiusofaround20-40˚A34,67\nimmediately after photoabsorption. Although our method captures the same\nqualitative behavior, direct comparison with experiment is limited by the size\nof the simulation box (12.42 ˚A), which restricts the maximum observable extent\nof delocalization73.\nRegardingstatisticalsampling,ourstudyincludes100independenttrajectories,which\nwe believe is sufficient to capture representative trends and provide meaningful en-\nsemble averaged insights. Nevertheless, we acknowledge that increasing the sampling\nsize could reveal additional rare events or subtle mechanistic features.\nTo address both the system size and sampling limitations, we have expanded the\ndiscussion in the final part of the Conclusions section, emphasizing the potential of\nmachine learning approaches to overcome these challenges (page 28, lines 521-523):\nOn the more methodological side, it is interesting to explore whether our simu-\nlations can be used to develop machine-learning approaches90−94 to examine the\nexcited-state dynamics of aqueous systemswhich would allow for expanding the\nsampling and statistics. This is an area we are actively pursuing and will be the\nsubject of a future study. Such methods would enable simulations with larger\nunit cells and a greater number of trajectories, potentially improving agreement\nwith experimental observations and revealing new physical insights.\n3.5: I am puzzled by the author’s finding of the H3O radical. This must\nbe a very unstable species (and the authors indeed discuss its transient\ncharacter) but I wonder if it is not more like a H3O+...e- contact “ion\npair”.\nWe thank the Reviewer for raising this important point. Our analysis reveals that\nwithin the HAT mechanism, the transient H O• species indeed forms, as the excess\n3\nelectron density is localized near one of the hydrogens bonded to oxygen, an aspect\nclarified in the Analysis Protocols section. We did not find the formation of H O+\n3\n21\n– e− pair within the HAT mechanism. In contrast, the PCET mechanism leads to\nthe formation of H O+, HO•, and a solvated electron. In this case a H O+ – e−\n3 3\npair emerges, mediated by surrounding water molecules that stabilize both species,\nas illustrated in Figure 2. We provide cube files for both mechanism in the Zenodo\nrepository.\nIn the present work, we chose to focus on the H O• radical because it has not yet been\n3\nexperimentally detected, and thus providing insights into this species. We also focus\non the HO• – H O+ pair, since it has been experimentally observed, allowing us to\n3\nbetter validate our computational approach. Nonetheless, the behavior and dynamics\nof the H O+ – e− represent an intriguing topic that warrants further investigation in\n3\nfuture studies.\nResponses to Reviewer #4:\n4.1: Their key analysis tool is to use “dilution factor” or “inverse partic-\nipation ratio” to quantify how many water molecules an excess electron\nis delocalized over. Note: strictly speaking the “PDFs” in Fig. 1BC are\ncounts, not probability distribution functions; they should be normalized,\nor maybe called probability counts.\nWe thank the reviewer for bringing this error to our attention. The reviewer is indeed\ncorrect, we used the IPR because it provides a continuous measure that reflects the\nnumber of water molecules involved in the excitation. We have now corrected Figures\n1B and 1C, changing the label from PDF to Probability counts.\n4.2: The conclusion is almost complete localization: the peak in Fig. 1B\nfor 1 water is 50, at 2 water molecules we are already down to ∼ 2 counts,\nor 4%. This need to be mentioned more clearly So, the further analysis\ndoes not change the view compare to with one water. The figures in Fig.2\nB could clearer, they are hard to see on the screen.\nWe thank the Reviewer for raising this concern. We would like to clarify that the\nIPR values shown on the x-axis of Figure 1B are continuous variables that represent\nthe number of water molecules participating in the excitation. Therefore, in order\nto explain the data point corresponding to “2” on the x-axis we should look at all\nthe configurations with IPR values near 2, which are more than 4% of the initial\nconditions. This point has now been clarified in the Analysis Protocol section (page\n34, line 640):\nticipate in the excitation.\nWe agree with the Reviewer that most excitations primarily involve a single water\nmolecule. Moreover, within the ensemble of simulations performed in this work, no\nsignificant differences in the photochemical behavior were observed as a function of\nthe number of water molecules involved in the excitation.\nIn the revised version of our manuscript, Figure 2 has been modified to enhance its\nclarity and readability.\n4.3: Defects account for 3/4 of the binding sites, mainly missing acceptor\nsites: a missing acceptor bond means the water molecule has a lone pair\non oxygen not bonded to an H from another water molecule, destabilizing\nthe n to sigma* transition to 8 eV. What is missing here is a further\naspect: the steric room available next to an unsatisfied acceptor bond also\nhelps the HO + H dissociation process by making room for the recoiling\nOH and H.\nWe thank the Reviewer for raising this important point, which was indeed missing\nfrom our original manuscript. Since this comment is closely related to a concern\nraised by another Reviewer (see point 1.4), we have addressed both points together.\nThe following modifications have been included in the revised manuscript in Section\n2.3 (pages 14–15, lines 264–274):\n... An interesting observation is that approximately 20% of the trajectories follow\nthe mechanism involving the formation of a hydrogen atom in an empty cavity\n(Figure 3B) with a lifetime of 13 fs, while the remaining 33% correspond to the\nalternative HAT pathway (Figures 3C and 3D) exhibiting a longer lifetime of\n34 fs. These results correlate well with the defect analysis presented in Figure\n1C. When the dissociating water molecule has a missing hydrogen-bond donor\n(AD or AAD defect), there is sufficient space available to allow an ultrafast\nhydrogen atom dissociation into a cavity (∼13 fs), in excellent agreement with\nvalues reported for a hydrogen-bonded water dimer in vacuum71,72. In contrast,\nwhen the dissociating molecule lacks a hydrogen-bond acceptor (ADD defect),\nthe surrounding hydrogen-bond network in the condensed phase restricts this\nmotion, delaying the non-radiative decay and promoting the formation of the\ntransient H O• species.\n3\n4.4: The definition of the decay lifetime on p. 9 is somewhat arbitrary but\nconsistent, sufficient tbuto distinguish processes differing by more than a\n23\nfactor of 10 in lifetime. May be this definition can be made more precise.\nWethankthereviewerforbringingthisimportantpointtoourattention. Inresponse,\nwe have revised the Analysis Protocols section to include explicit definitions of the\nlifetimes and quantum yields (page 35-36, lines 680–688).\nBased on the preceding analysis, each trajectory can be unambiguously assigned\nto a specific mechanism. This allows us to determine the quantum yield and\nlifetime associated with each mechanism using the following expression:\n#Trajectories with M\nQY = ∗100 (6)\nM\n#Total trajectories\n(cid:80)NM tM\nτ = i i (7)\nM\n#Total trajectories with M\nwhere M represent the mechanism (HAT or PCET), QY,τ are the quantum\nyields and lifetimes, N is the number of trajectories with the M mechanism\nM\nand tM is the time at which the trajectory i with the M mechanism finds the\ni\n”crossing point”.\n4.5: In the PCET mechanism (n→sigma* H2O + H2O→ HO + e- +\nH+ +H2O→ HO + e- + H3O+) , the e- quickly becomes localized in\na trapped state. The OH and H3O+ become quite well separated (6 ˚A\nin Fig. 5D). Their theoretical study reveals three distinct minima in the\nHO* and (H3O+) distance which are in nice agreement with the study of\nthe recent electron diffraction experiments upon photoionization. These\ncould be rationalized by the formation of contact ion radical pair and the\nseparation mediated by the solvent. The results are not surprising, the\nnew aspect is that the entire process can be modelled from photoexcitation\nup to the localized electron directly capturing the solvent dynamics.\nWe thank the Reviewer for their comment and for emphasizing the main findings of\nour study. We believe that the present work represents a significant step forward in\nunderstanding the complex photochemistry of liquid water. For the first time, we\ndemonstrate how all the key processes: including water dissociation, formation of\nhydroxyl and hydronium radicals, electron generation, hydronium ion and collective\nwatermotionsevolvewithinaunifiedframework. Ourresultsnotonlyreproducethese\nintricate phenomena starting from neutral liquid water but also provide new physical\ninsights, particularly into the very early stages of the photochemical process and the\nearly birth of the hydrated electron. This initial regime has been largely unexplored\n24\nin previous theoretical studies, and we hope that our findings will stimulate further\ninvestigations under a broader range of aqueous conditions.\nResponses to Reviewer #5:\n5.1: I co-reviewed this manuscript with one of the reviewers who provided\nthe listed reports. This is part of the Nature Communications initiative\nto facilitate training in peer review and to provide appropriate recognition\nfor Early Career Researchers who co-review manuscripts.\nWe sincerely thank the Reviewer for their time and effort in evaluating our work. We\ntruly appreciate all the feedback which helped us to improve the quality and clarity\nof our manuscript.\n[1] Takeshi Yanai, David P Tew, and Nicholas C Handy. A new hybrid exchange–\ncorrelation functional using the coulomb-attenuating method (cam-b3lyp).\nChemical physics letters, 393(1-3):51–57, 2004.\n[2] Joost VandeVondele and Ju¨rg Hutter. Gaussian basis sets for accurate calcula-\ntionsonmolecularsystemsingasandcondensedphases. TheJournalofchemical\nphysics, 127(11), 2007.\n[3] Javier Segarra-Marti, Daniel Roca-Sanjuan, Manuela Merchan, and Roland\nLindh. Onthephotophysicsandphotochemistryofthewaterdimer. TheJournal\nof Chemical Physics, 137(24), 2012.\n[4] Anupriya Kumar, Maciej Ko(cid:32)laski, Han Myoung Lee, and Kwang S Kim. Pho-\ntoexcitation and photoionization dynamics of water photolysis. The Journal of\nPhysical Chemistry A, 112(24):5502–5508, 2008.\n[5] Dario Rocca, Ralph Gebauer, Yousef Saad, and Stefano Baroni. Turbo charging\ntime-dependent density-functional theory with lanczos chains. The Journal of\nChemical Physics, 128(15), 2008.\n[6] Jan-Robert Vogt, Michael Schulz, Rafael Souza Mattos, Mario Barbatti, Maur-\nizio Persico, Giovanni Granucci, Jurg Hutter, and Anna Hehn. A density func-\ntional theory and semiempirical framework for trajectory surface hopping on\nextended systems. Journal of Chemical Theory and Computation, 2025.\n[7] Rachel Crespo-Otero and Mario Barbatti. Recent advances and perspec-\ntives on nonadiabatic mixed quantum–classical dynamics. Chemical reviews,\n118(15):7026–7068, 2018.\n[8] Antonio Prlj, Jack T Taylor, Jiˇr´ı Janoˇs, Petr Slav´ıˇcek, Federica Agostini, and\nBasile FE Curchod. Best practices for nonadiabatic molecular dynamics simula-\ntions. arXiv preprint arXiv:2508.05263, 2025.\n[9] Rachel Crespo-Otero and Mario Barbatti. Cr (co) 6 photochemistry: Semi-\nclassical study of uv absorption spectral intensities and dynamics of photodisso-\nciation. The Journal of chemical physics, 134(16), 2011.\n[10] Felix Plasser, Rachel Crespo-Otero, Marek Pederzoli, Jiri Pittner, Hans Lis-\nchka, and Mario Barbatti. Surface hopping dynamics with correlated single-\nreference methods: 9h-adenine as a case study. Journal of chemical theory and\ncomputation, 10(4):1395–1405, 2014.\n26\n[11] Lea M Ibele, Pedro A Sanchez-Murcia, Sebastian Mai, Juan J Nogueira, and\nLeticia Gonza´lez. Excimer intermediates en route to long-lived charge-transfer\nstates in single-stranded adenine dna as revealed by nonadiabatic dynamics. The\nJournal of Physical Chemistry Letters, 11(18):7483–7488, 2020.\n[12] John M Herbert, Xing Zhang, Adrian F Morrison, and Jie Liu. Beyond time-\ndependent density functional theory using only single excitations: Methods for\ncomputationalstudiesofexcitedstatesincomplexsystems. Accountsofchemical\nresearch, 49(5):931–941, 2016.\n[13] Jiri Suchan, Jiri Janoˇs, and Petr Slavicek. Pragmatic approach to photodynam-\nics: Mixed landau–zener surface hopping with intersystem crossing. Journal of\nChemical Theory and Computation, 16(9):5809–5820, 2020.\n[14] Nina Toki´c, Tomislav Piteˇsa, Antonio Prlj, Marin Sapunar, and Naa Doˇsli´c.\nAdvantages and limitations of landau-zener surface hopping dynamics. Croatica\nchemica acta, 97(4):P1–P11, 2024.\n[15] WeiweiXieandWolfgangDomcke. Accuracyoftrajectorysurface-hoppingmeth-\nods: Test for a two-dimensional model of the photodissociation of phenol. The\nJournal of Chemical Physics, 147(18), 2017.\n[16] Gonzalo D´ıaz Mir´on, Carlos R Lien-Medrano, Debarshi Banerjee, Uriel N\nMorzan, Michael A Sentef, Ralph Gebauer, and Ali Hassanali. Exploring the\nmechanisms behind non-aromatic fluorescence with the density functional tight\nbinding method. Journal of Chemical Theory and Computation, 20(9):3864–\n3878, 2024.\n27"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I thank the authors for thoroughly addressing my comments. I recommend publication of the manuscript in Nature Communications.\nIn addition, it would be helpful if the authors could clarify the following point:\nIn Figures 1B and 1C, the y-axis is labeled “Probability counts”, which is not normalized and therefore does not sum to 1. However, the caption states that these panels show the probability distribution function, which implies a normalized quantity. The figure would be easier to interpret if this inconsistency were resolved.", "reviewer_2": "In the revised version of the manuscript \"The Photochemical Birth of the Hydrated Electron in Liquid Water\" the authors have addressed main issues pointed out in the first round of peer review: provided validation of ROKS for the HAT mechanism, discussed and evaluated the importance non-adiabatic effects and performed simple non-adiabatic (Landau-Zener theory).\nIn addition, they referred to more relevant publications.\nI will support the publication of this work in Nature Communications if the discussion of the benchmark of ROKS against CASPT2 is done properly and fairly: their evaluation of the methods is way too optimistic, contradicting the data.\nUnfortunately, I have reasons to believe it is not the case. Panels D and E in Figure 7 show what I was afraid such a benchmark would show: very significant qualitative differences. Yes, the barrier height is indeed similar, but the geometry of TS is embarrassingly different (which is ingeniously, and reasonably, downplayed). However, the values in between the reactant and the conical intersection are nowhere near. Say at the coordinate value of 2 (units are not given at the axis, please add it) the the CASPT2 energy is ca. 2 times larger than that of ROKS, the difference being immense 2 eV. The shape of conical intersection region is dramatically different: it is clear from the figure that non-adiabatic effects at the multireference level of theory will be significantly larger since the gap squared appears in the exponential in the Landau-Zener theory.\nI believe it's non proper to write about \"close agreement with TDDFT and CASPT2\" and so on in the last paragraph of Section 3. Instead, the authors should come up with realistic evaluation of their approach as at most semiqualitative.", "reviewer_3": "The authors satisfactorily addressed my comments, so I believe the ms. may be published.", "reviewer_4": "This is a theoretical simulation study on the birth of a solvated electron following all the steps thorough from the photoexcitation ,i.e. from the ground state to the excited state and then the formation of the hydrated electron up to localization including readical formation.\nNew experimental techniques have been developed incluing e.g. ultrafast electron diffraction which allows to monitor this process experimentally. Thus a theoretical study which involes all steps is of relevance and has not been reported before.\nThe data analysis and the choice of the theoretical methods should be checked in detail by theoretical referees.\nBeing outside of this field I can only state that the results are interesting and noteworthy.", "reviewer_5": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "Comments for: Simulating the\nPhotochemical Birth of the Hydrated\nElectron in Liquid Water\nJanuary 2026\nDear Reviewers,\nWe sincerely appreciate the time and effort that the Reviewers have dedicated to eval-\nuating our manuscript and providing thoughtful and constructive feedback. In the\nfollowing sections, we provide detailed responses to each comment raised by the Re-\nviewers. The corresponding revisions are indicated in the file “manuscript marked”,\nwhere page and line numbers are referenced for clarity.\nResponses to Reviewer #1:\n1.1: I thank the authors for thoroughly addressing my comments. I rec-\nommend publication of the manuscript in Nature Communications. In\naddition, it would be helpful if the authors could clarify the following\npoint: In Figures 1B and 1C, the y-axis is labeled ”Probability counts”,\nwhich is not normalized and therefore does not sum to 1. However, the\ncaption states that these panels show the probability distribution func-\ntion, which implies a normalized quantity. The figure would be easier to\ninterpret if this inconsistency were resolved.\nWesincerelythanktheReviewerforthiscommentandforrecommendingourmanuscript\nfor publication. We apologize for this oversight and have now addressed the issue in\nthe caption of Figure 1. The following modification has been made in the revised\nmanuscript (page 8, Figure 1):\n1\n... Panel B: Probability Distribution Function (PDF) Counts of water\nmolecules involved in the excitation using the Inverse Participation Ratio of\nthe spin densities for all the conformations. The insets show two examples\ninvolving 1 and 5 water molecules. Panel C: PDF Probability Counts of the\ndifferent defects in the Hydrogen Bond (HB) Network of all the water molecules\ninvolved in the initial excitation. ...\nResponses to Reviewer #2:\n2.1: In the revised version of the manuscript ”The Photochemical Birth of\nthe Hydrated Electron in Liquid Water” the authors have addressed main\nissues pointed out in the first round of peer review: provided validation of\nROKS for the HAT mechanism, discussed and evaluated the importance\nnon-adiabatic effects and performed simple non-adiabatic (Landau-Zener\ntheory). In addition, they referred to more relevant publications. I will\nsupport the publication of this work in Nature Communications if the dis-\ncussion of the benchmark of ROKS against CASPT2 is done properly and\nfairly: their evaluation of the methods is way too optimistic, contradicting\nthe data. Unfortunately, I have reasons to believe it is not the case. Pan-\nels D and E in Figure 7 show what I was afraid such a benchmark would\nshow: very significant qualitative differences. Yes, the barrier height is\nindeed similar, but the geometry of TS is embarrassingly different (which\nis ingeniously, and reasonably, downplayed). However, the values in be-\ntween the reactant and the conical intersection are nowhere near. Say at\nthe coordinate value of 2 (units are not given at the axis, please add it)\nthe the CASPT2 energy is ca. 2 times larger than that of ROKS, the\ndifference being immense 2 eV. The shape of conical intersection region\nis dramatically different: it is clear from the figure that non-adiabatic ef-\nfects at the multireference level of theory will be significantly larger since\nthe gap squared appears in the exponential in the Landau-Zener theory.\nI believe it’s non proper to write about ”close agreement with TDDFT\nand CASPT2” and so on in the last paragraph of Section 3. Instead, the\nauthors should come up with realistic evaluation of their approach as at\nmost semiqualitative.\nWe thank the Reviewer for bringing this important point to our attention. We would\nfirst like to clarify how the validation against CASPT2 was carried out. In the\nCASPT2 study by Segarra et al.[1], the authors reported a minimum energy path\n(MEP) analysis but provided optimized geometries only for the ground state mini-\n2\nmumandtheconicalintersection. Startingfromthesetwogeometries, weconstructed\na linear interpolated path (LIP) and evaluated the ROKS and TDDFT energies along\nthis coordinate. As a consequence, a direct, point-by-point comparison between the\npublished MEP with CASPT2 and LIP pathways (that we construct with ROKS and\nTDDFT) is not possible. In particular, only the initial geometry shown in Figure 7\npanels D-F corresponds to the calculation with the same nuclear configuration at all\nlevels of theory. We recognize that this limitation was not sufficiently clear in the\noriginal version of the manuscript, and we have therefore revised the text to explic-\nitly state this point. The corresponding modification is reported below (page 25, lines\n458-462):\nUsing both the ground-state and conical intersection optimized geometries at\nthis high level of theory, we computed the potential energy surfaces along the\nLinear Interpolated Path (LIP) with ROKS (Figure 7E) and TDDFT (Figure\n7F) employing the CAM-B3LYP functional. It is worth emphasizing that the\nnumerical values along horizontal axes for the MEP and LIP coordinates do not\ncorrespond to identical nuclear configurations, except for the initial point shown\nin Figure 7D-F. We thus focus on the overall trends of ground and excited state\npotential energy surfaces.\nWe agree with the Reviewer that our original statement regarding the quantitative\nvalidation was perhaps overly optimistic. In particular, the difference of approxi-\nmately 0.8 ˚A in the HO•–H• distance at the conical intersection between ROKS and\nCASPT2 highlights a significant discrepancy in the predicted intersection topology.\nTo address this point more carefully, we have revised the discussion to better reflect\nthequalitativeratherthanquantitativenatureofthiscomparison. Thecorresponding\nmodification is reported below (pages 25-26, lines 465-470):\nAlthough this difference is not insignificant, a 2 ˚A separation already indicates\na fully dissociated configuration, confirming that ROKS and CASPT2 capture\nthe same underlying photodissociation physics Although an O-H distance of 2 ˚A\ncould be considered as dissociated geometry, the topology of the conical intersec-\ntion differs between the two methods. Despite this difference, ROKS reproduces\nthe same essential physical process, namely the hydrogen atom dissociation.\nFinally, we have revised and expanded the concluding discussion of the validation\nsection to provide a more balanced assessment of the strengths and limitations of the\nROKS approach. These changes are reported below (pages 26–27, lines 476–494).\nROKS approach. We demonstrate that ROKS captures the essential physics\nof both the HAT and PCET photochemical pathways in liquid water, in close\nagreement with TDDFT and CASPT2 benchmarks. This consistency across\ndifferent mechanisms and details of the underlying electronic structure supports\nthe robustness of our methodology and reinforces confidence in its predictive\ncapability for modeling the excited-state dynamics of liquid water. Based on the\ngas phase water dimer model, ROKS does not reproduce the conical intersection\ngeometry for the HAT mechanism. In particular, a difference of approximately\nof 0.8 ˚A is observed in the HO•-H• distance at the conical intersection when\ncompared with CASPT2. Nevertheless, ROKS qualitatively describes the same\nphotodissociation pathway, reproducing both the slope of the excited-state poten-\ntial energy surface and the presence of the relevant chemical species along the\ndissociation coordinate. While such differences may lead to deviations in the\npredicted HAT lifetimes in vacuum, we expect their impact to be reduced in the\ncondensed phase. Indeed, in bulk liquid water, ROKS yields lifetimes and quan-\ntum yields that are in good agreement with experimental observations. Over-\nall, the extensive validation presented here demonstrates the robustness of our\nmethodology and reinforces confidence in its predictive capability for modeling\nthe excited-state dynamics and photochemistry of liquid water. We nonetheless\ncaucion that more systematic studies with multireference methods for condensed\nphase aqueous systems needs to be explored in the future.\nWe believe that these revisions more accurately convey the scope of the validation\npresented in our manuscript.\nResponses to Reviewer #3:\n3.1: The authors satisfactorily addressed my comments, so I believe the\nms. may be published.\nWe sincerely thank the Reviewer for their positive assessment of our responses and\nfor supporting the publication of our manuscript.\nResponses to Reviewer #4:\n4.1: This is a theoretical simulation study on the birth of a solvated\nelectron following all the steps thorough from the photoexcitation ,i.e. from\n4\nthe ground state to the excited state and then the formation of the hy-\ndrated electron up to localization including readical formation. New ex-\nperimental techniques have been developed incluing e.g. ultrafast electron\ndiffraction which allows to monitor this process experimentally. Thus a\ntheoretical study which involes all steps is of relevance and has not been\nreported before. The data analysis and the choice of the theoretical meth-\nods should be checked in detail by theoretical referees. Being outside of\nthis field I can only state that the results are interesting and noteworthy.\nWe sincerely thank the Reviewer for their thoughtful and positive assessment of our\nwork. We believe that our study provides complementary insight into these experi-\nmentsbyofferingamolecularlevelinterpretationoftheprocessesinvolvedinthebirth\nand localization of the hydrated electron, which we hope will be useful for guiding\nand interpreting future experimental investigations.\nResponses to Reviewer #5:\n5.1: I co-reviewed this manuscript with one of the reviewers who provided\nthe listed reports. This is part of the Nature Communications initiative\nto facilitate training in peer review and to provide appropriate recognition\nfor Early Career Researchers who co-review manuscripts.\nWe sincerely thank the Reviewer for their time and effort in evaluating our work. We\ntruly appreciate all the feedback which helped us to improve the quality and clarity\nof our manuscript.\nReferences\n[1] Javier Segarra-Marti, Daniel Roca-Sanjuan, Manuela Merchan, and Roland\nLindh. On the photophysics and photochemistry of the water dimer. The Journal\nof Chemical Physics, 137(24), 2012.\n5"}]

Decision letter