Simulating the ghost: quantum dynamics of the solvated electron

The solvated electron, e−(aq), is the smallest anion and a fundamental reducing agent with implications for electrochemistry, photochemistry, high-energy chemistry, and biology, where its non-equilibrium precursor contributes to radiation damage to DNA. Despite its apparent simplicity, modeling e−(aq) is challenging: standard density functional approximations suffer from delocalization error and struggle with radicals and even neat water unless carefully chosen. Many-body quantum chemistry methods such as MP2 can accurately describe water interactions but are prohibitively expensive for extended simulations. Previous spin-unrestricted MP2 under periodic boundary conditions enabled picosecond-scale MD of bulk e−(aq), supporting a cavity structure rather than non-cavity models. However, comprehensive characterization requires much longer timescales and inclusion of nuclear quantum effects (NQEs), which are known to be important for electron relaxation dynamics, redox properties, and can qualitatively alter aqueous behavior by sampling classically inaccessible configurations. This work addresses these challenges by developing a machine-learning interatomic potential trained on MP2 data to enable long-timescale quantum dynamics of the hydrated electron, capturing localization, structure, diffusion, and spectroscopy with MP2-level accuracy at much lower computational cost.

Prior studies highlighted limitations of standard DFT for radicals due to delocalization errors, while certain functionals can approximate high-level benchmarks and experiments for water. Many-body methods (e.g., MP2) accurately reproduce two- and three-body interactions in water and electron–water clusters, consistent with CCSD(T) benchmarks. Implementation of periodic-boundary UMP2 energies and forces enabled accurate ab initio MD for bulk e−(aq), with picosecond trajectories supporting a cavity model. NQEs have been shown to significantly influence electron relaxation in clusters and affect redox and attachment processes in condensed phase; they can enable sampling beyond classical configurations, leading to qualitative changes in aqueous systems. These insights motivate a method that combines high-level electronic accuracy with feasible quantum nuclear dynamics over longer timescales.

• Electronic-structure training data: Forces and energies computed at periodic spin-unrestricted MP2 level were used to train a Behler–Parrinello Neural Network (BPNN) interatomic potential. The training set included configurations from both classical and quantum MD to improve transferability. Benchmarking against UMP2 energetics is reported in Supplementary Methods.
• Machine-learning potential: The NN represents total energy as a sum of atom-centered contributions, encoding the effect of the excess electron implicitly via water–water structural correlations stabilized by the electron. The potential is trained for a fixed system size/concentration of one excess electron and should not be used for different sizes without retraining.
• Simulation protocol: Systems were equilibrated using a water-only NN potential, then switched to the e−-trained potential to model electron injection. Nuclear quantum dynamics were performed using thermostatted ring polymer MD (TRPMD). Thirty independent trajectories were run to probe localization dynamics and diffusion over several hundred picoseconds in total.
• Electronic observables: Spin densities ρspin(r) and band gaps were not produced directly by BPNN; instead, for selected frames spin densities and gaps were evaluated with DFT using a hybrid functional. Band gaps and gyration radii of spin density distributions were analyzed, including e−–H coordination states.
• Structural analysis: Single- vs twin-cavity structures were identified; radial distribution functions g(e−–H) and g(e−–O) were computed relative to the spin-density center. Hydrogen-bond (HB) network changes were quantified using a probabilistic HB-state classifier.
• Spectroscopy: Electronic absorption spectra were computed via time-dependent DFT for single- and twin-cavity structures. Vibrational density of states (VDOS) were computed via TRPMD for molecules within 3.5 Å of the spin-density center for H2O, D2O, and HOD, and compared with classical MD. Isotopic effects and mixed-isotope features were analyzed, including predicted bending doublets for HOD near the cavity.
• Validation: Cross-validation showed errors for twin-cavity structures comparable to training data; retraining with extended sets including validation structures gave similar behavior, supporting the MP2-consistent stabilization of twin cavities by NQEs.

• Localization dynamics: Upon switching to the e− potential, the initially delocalized excess electron undergoes a pre-solvation stage and localizes into a cavity within ~1 ps in most of 30 TRPMD trajectories, consistent with experiment and prior ab initio MD.
• Cavity structure: The localized electron cavity involves 4–5 water molecules, each donating one OH toward the cavity center. The mean e−–O coordination number is 4.5 (close to prior DFT estimate 4.7). Spin-density distributions exhibit characteristic negative density near H atoms.
• Radial distribution functions: NQEs broaden first-shell peaks in g(e−–O) and g(e−–H) due to zero-point motion and reduce longer-range order (blurred second-neighbor peaks). Oxygen RDFs have non-zero values within 1.5 Å of the spin-density center, supporting a fluxional picture of e−(aq). Heavy water (D2O) RDFs are similar to classical simulations, consistent with compensating quantum effects in D2O.
• Hydrogen-bond network: Cavity formation increases undercoordinated species; per cavity, singly-donating molecules increase by ~2 and singly-accepting by ~3.
• Electronic structure correlations: For single cavities, band gap inversely correlates with spin-density gyration radius and directly with e−–H coordination. Average band gap ≈ 2.6 eV. Anisotropy (including single and double cavities) is 0.0423 ± 0.0349, indicating twin cavities are rare and contribute little to averages.
• Twin-cavity states and NQEs: Quantum simulations (H2O) reveal twin-cavity configurations with gyration radius 2–3 Å and anomalously low gaps ≤1.5 eV—absent in classical simulations. Twin cavities enable a transient diffusion mechanism involving two adjacent voids and electron shuttling.
• Diffusion mechanisms and rates: A non-transient diffusion mechanism involves HB reorganization of neighboring waters (molecules A, B, C). The transient mechanism (unique to quantum H2O) features double cavities and electron transfer between them. Despite their presence, transient events contribute little to overall diffusion. Diffusion coefficients: quantum 0.40 ± 0.03 Ų/ps; classical 0.36 ± 0.04 Ų/ps; both agree reasonably with experiment 0.475 ± 0.048 Ų/ps.
• Absorption spectra: Single-cavity absorption maximum (s→p) at ~2.2 eV, in reasonable agreement with experimental ~1.7 eV. Twin-cavity spectra show a broad feature with maximum <1.0 eV (s-to-s between neighboring cavities). Given <10% twin-cavity population, their spectral contribution is marginal; they may appear in transient bleaching at lower frequencies.
• Vibrational spectroscopy (VDOS proxy for RR): For bending modes, TRPMD reproduces peak positions and electron-induced downshifts within ~10 cm−1 for H2O and D2O (exp. downshifts ~30 cm−1 for H2O and ~20 cm−1 for D2O). Classical MD shows >50 cm−1 blue shift and negligible e−-induced shift. In stretches, TRPMD shows a red shift ~150 cm−1 (exp. ~200 cm−1); line shapes less reliable, but trends consistent with increased HB defects (singly-donating OH red shift ~200 cm−1).
• Isotopic mixtures: For HOD near the cavity, predicted bending doublet at 1338 and 1399 cm−1 (splitting ~61 cm−1), consistent with experimental splitting (~60 cm−1), attributed to H–O–D–e vs D–O–H–e bending. Equilibrium isotopic segregation favors OH over OD at the cavity; experimental equilibrium constant K≈1.6 supports preferential H2O enrichment around e−.

By training a BPNN potential on MP2 reference data and performing TRPMD simulations, the study captures the essential physics of the hydrated electron without explicitly modeling the electron in the force field. The model reproduces known localization dynamics from a delocalized to pre-solvated to cavity state, substantiating the cavity picture. Inclusion of NQEs explains broadened structural distributions, reduced long-range order, and reveals rare twin-cavity configurations that enable a transient diffusion pathway observable only in quantum H2O simulations. Quantitative agreement with experiment is achieved for diffusion coefficients and for resonance-Raman-relevant bending shifts; stretching shifts are reproduced semi-quantitatively, consistent with known limitations of approximate quantum dynamics and the VDOS proxy. Electronic absorption features for single cavities align reasonably with experimental maxima; twin-cavity features are predicted at lower energies with negligible macroscopic impact due to their low population but could be probed by transient techniques at low frequencies. Heavy water results show near-compensation of competing quantum effects, rationalizing similarities to classical behavior. Overall, the findings validate that an MP2-accurate ML potential, combined with quantum nuclear dynamics, can faithfully describe structure, energetics, dynamics, and spectroscopic signatures of e−(aq), while uncovering NQE-enabled transient structural motifs with limited macroscopic impact.

This work delivers an MP2-accurate machine-learning interatomic potential that implicitly encodes the effect of the hydrated electron, enabling long-timescale quantum dynamics simulations at affordable cost. The approach reproduces cavity formation and localization dynamics, quantifies structure (4–5 water coordination, mean e−–O coordination ~4.5), diffusion mechanisms and coefficients consistent with experiment, and vibrational and electronic spectroscopies with good to semi-quantitative accuracy. NQEs are shown to broaden structural features, diminish long-range order, and stabilize rare twin-cavity states that mediate a transient diffusion pathway unique to quantum H2O. The low population of twin cavities limits their macroscopic impact but suggests experimental detection via low-frequency transient bleaching. Future work could expand training data to larger system sizes and varying electron concentrations, interfaces, and different thermodynamic conditions; refine quantum dynamics methods for improved stretching-region spectroscopy; and integrate direct electronic observables within ML frameworks to avoid post hoc electronic-structure evaluations.

• Fixed electron concentration: The ML potential is trained for a specific system size and one excess electron; changing size without retraining effectively alters electron concentration and can yield nonphysical results.
• Indirect electronic properties: Spin densities and band gaps are not produced by the BPNN and require separate hybrid-DFT calculations on selected frames.
• Training coverage: Twin-cavity configurations were underrepresented in original MP2 training; although cross-validation and retraining confirm trends, residual uncertainty remains in sparsely sampled regions.
• Approximate quantum dynamics: TRPMD, while effective, can incur errors (notably in stretching region up to ~100 cm−1), and VDOS may not perfectly represent resonance Raman line shapes.
• Spectroscopic comparisons: TDDFT-based absorption maxima and intensities carry functional-dependent uncertainties; low twin-cavity populations complicate experimental detection.
• Timescale and scope: Although several hundred picoseconds were sampled, rare-event statistics and broader conditions (e.g., varying concentrations, interfaces) were not explored.