Evolution of Aqueous Electron with Varying Temperature
The temperature-dependent properties of the aqueous electron have been extensively studied using mixed quantum-classical simulations in a wide range of thermodynamic conditions based on one-electron pseudopotentials. While the cavity model appears to explain most of the physical properties of the aqueous electron, only a non-cavity model has so far been successful in accounting for the temperature dependence of the absorption spectrum. Here, we present an accurate and efficient description of the aqueous electron under various thermodynamic conditions by combining hybrid functional-based molecular dynamics, machine learning techniques, and multiple time-step methods. Our advanced simulations accurately describe the temperature dependence of the absorption maximum in the presence of cavity formation. Specifically, our work reveals that the red shift of the absorption maximum results from an increasing gyration radius with temperature, rather than from global density variations as previously suggested.
Introduction<!>Conclusions<!>Methods
<p>The aqueous electron (e − aq ) is the simplest condensed-phase anion in nature, consisting of an excess electron solvated in liquid water. [1][2][3] In the early 1960s, Hart and Boag discovered e − aq through its absorption spectrum, 4 which is characterized by a main peak resulting from transitions between the s ground state and the p excited states. 2 Since then, efforts in both theory and experiment have led to a better understanding of the structural, dynamical, and electronic properties of e − aq . [5][6][7][8][9][10][11][12][13][14] However, despite nearly a half-century of interest in e − aq , many questions remain controversial. 2 In particular, the question of whether e − aq occupies a cavity has been highly debated. [15][16][17][18] This lack of consensus mostly originated from the use of mixed quantum-classical simulations based on the socalled one-electron pseudopotential method, where the electron is treated quantum mechanically and the other interactions classically. [19][20][21][22][23][24] In more recent years, cavity formation has drawn strong support, both from molecular dynamics simulations based on more accurate electronic structure methods [25][26][27][28] and from the measurement of partial molar volumes. 29 Despite this progress, the variation of the optical absorption spectrum of e − aq with temperature has remained elusive. 2 More specifically, the absorption maximum E max (T ) shows a red shift of 2.4-2.8 meV/K at constant pressure, 30,31 which has not yet found a consistent explanation. A similar shift with temperature is observed at constant density, 30,31 suggesting that density variations only play a minor role. Within the one-electron pseudopotential scheme, 21 many theoretical attempts have been undertaken to understand the dependence of E max upon varying thermodynamic conditions. Turi and Borgis (TB) developed a potential leading to cavity formation and giving a fair absorption spectrum. 19,20 The TB model reproduces the experimental temperature dependence of E max , 32 but this agreement was demonstrated to originate from a density effect, 32,33 unlike in the experiment. 30,31 A more elaborate polarizable electron-water pseudopotential model also failed in reproducing the temperature dependence. 2,34 Subsequently, Larsen, Glover, and Schwartz (LGS) developed a non-cavity model 22 , which was found to describe E max (T ) at fixed density in a qualitatively correct fashion, apart from the slope with temperature being overestimated by a factor larger than two. 33 A more accurate description of e − aq can be achieved by describing the electronic structure within density-functional theory. In the early 2000s, ab initio molecular dynamics (AIMD) simulations based on the generalized gradient approximation found that decreasing the density resulted in a red shift of the absorption spectrum, in qualitative agreement with the TB model. 35,36 However, it is by now well known that an accurate description of e − aq requires more advanced methods to overcome delocalization errors and to capture subtle correlations effects. 11,25,26,28,37 To limit the computational cost of such methods, machine learning schemes 38 offer a viable strategy to simulate the quantum dynamics of e − aq with inexpensive force fields. 28 Here, we combine hybrid-functional molecular dynamics and machine learning techniques to elucidate the evolution of the adsorption spectrum of the aqueous electron with temperature. Using the energies and forces determined at the hybrid functional level, we train an inexpensive force field through a neural network. 38 The molecular dynamics are accelerated by a multiple timestep method. 39,40 The long time scale enabled by such a scheme allows us to accurately determine the dependence of the structural and electronic properties of e − aq on temperature. The present results lead to cavity formation and to a temperature dependence of the absorption maximum in agreement with experiment, thereby demonstrating the strength of the applied methodology.  in a quantum-mechanics/molecular-mechanics (QM/MM) set-up to calculate absorption spectra through time-dependent density functional theory (TDDFT) (see Methods for more details).</p><p>The electronic structure scheme used in our work correctly captures the temperature dependent density of water, whereas the simple point charge potential used in the one-electron pseudopotential model is unable to accurately reproduce such a property (see Fig. S2). Our MLP performs well for temperatures up to about 600 K (see Fig. S2), but the description deteriorates at higher temperatures, where the low density regime (< 0.6 g/cm 3 ) might be affected by the lack of longrange interactions in the MLP. Therefore, we here focus on the temperature range from 300 to 573 5 K.</p><p>In Fig. 2(a,b), we give radial distribution functions of oxygen (g eO ) and hydrogen atoms (g eH )</p><p>with respect to the electron for various temperatures. We find that the excess electron localizes at all temperatures, forming a cavity structure. Increasing the temperature appears to accelerate the localization process, as can be seen from the evolution of the band gap and the gyration radius (Fig. S6). At all temperatures, the hydrated electron is characterized by a similar geometrical structure, with the hydrogen atoms of the nearby water molecules pointing to the center of the electron. In response to a temperature increase, the vibrational motions globally broaden the features in g eO To shed light onto the temperature dependence of E max (T ), we carry out TDDFT calculations of the absorption spectrum in a QM/MM set-up and average over the structural configurations achieved with the MTS-AIMD simulations (see Methods). For each configuration, we consider the transitions between the ground state and the lowest five excited states. E max is then obtained as the peak position of the convoluted absorption spectrum and mainly results from the s-p transitions (Fig. S5). We find a close to linear dependence of E max on temperature. The calculated E max (T )</p><p>shows a substantial redshift of 2.2 ± 0.2 meV/K with temperature, in good agreement with experimental data (2.4 meV/K) 30 obtained in the same thermodynamic conditions (see Fig. 3). The calculated E max lies higher by about 0.4 eV than the experimental values. This offset should be assigned to the energy levels of the unoccupied p orbitals, which are known to sensitively depend on the adopted electronic structure method. 42 Indeed, at 300 K, the transition from the ground state to the vacuum level, i.e. the vertical binding energy (VBE), is found at 3.26 ± 0.17 eV, in excellent agreement with its experimental counterpart at 3.2 eV. 43 The dependence of the VBE on temperature is given in Fig. S8.</p><p>Since temperature and density vary concurrently in our simulations, the effect originating the variation of E max remains undetermined. To clarify this issue, we carry out two additional series of simulations. In the first series, we vary the temperature (423, 473 and 523 K), while keeping the density fixed at 0.977 g/cm 3 . In the second series, various densities (1.029, 0.977 and 0.907 g/cm 3 ) are considered at the same temperature of 473 K. In Fig. 3(b), we collect E max obtained in these various thermodynamics conditions. It appears clearly that the temperature affects E max more significantly than the density. At constant density, E max redshifts with a slope of 1.9 ± 0.4 meV/K, comparable to the slope of 2.2 ± 0.2 meV/K found in simulated constant-pressure conditions (red circles in Fig. 3(a)). In contrast, E max redshifts only slightly with density at constant temperature.</p><p>We estimate rates of 0.34, 0.36, and 0.15 eV per g/cm 3 at 423, 473, and 523 K, respectively, in accord with the measured value of 0.2 eV per g/cm 3 . 44 By the same token, the structure of e − aq is more affected by the temperature than by the density, as evidenced by the radial distribution functions, the gyration radii, and the distributions of R eH and θ (see Figs. S7 and S8). In particular, the gyration radii remain almost constant upon density variations at fixed temperature, but they are highly sensitive to temperature. Overall, the gyration radius emerges as the key factor for determining E max .</p><p>To investigate the origin of the temperature dependence of E max , we consider a simple onedimensional square-well model with finite barriers V . Here, we assume the separation between the ground-state level and the onset of the continuum to correspond to the vertical binding energy of e − aq , and the half-width of the well to the electron gyration radius r gyr (see Fig. 3(c)). We numerically solve the Schrödinger equation and obtain the excitation energy ∆E to the first excited state as a function of both VBE and r gyr (see Fig. 3(d)). Using the values for VBE and r gyr as obtained from our ab-initio simulations (see black dots in Fig. 3(d)), we obtain ∆E redshifting with temperature by 1.6 meV/K (red triangles in Fig. 3(a)), in reasonably good agreement with the trends seen in the experiment and in our ab initio simulations. This model implies that VBE and ∆E decrease with increasing gyration radius. Since our ab initio simulations reveal a clear increase of the gyration radius with temperature, this provides an interpretation for the temperature dependence of E max . As a further outcome, this model puts forward the prediction that VBE decreases with temperature (Fig. S8), which can be validated by future experiments.</p><p>Using the one-electron pseudo-potential method, the TB (green circles in Fig. 3(a)) and the</p><p>LGS model (blue circles) have been able to reproduce the experimental trend of E max with tem-perature. However, these two models rely on very different underlying physical mechanisms to support the observed behavior. On the basis of the TB model, Nicolas et al. found that the variations of E max originate from a density rather than from a temperature effect, obtaining a redshift of 1.4 eV per g/cm 3 for increasing density and a vanishing shift for temperature variations at constant density (blue squares in Fig. 3). 32,33 These results are not supported by experimental observations, which measure redshifts of 0.2 eV per g/cm 3 at constant temperature 44 and of 2.2 meV/K at constant density. 31 The inaccuracy of the TB model could result from the underestimation of the interaction between the electron and its neighboring water molecules, as suggested by the height of the first peaks in the radial distribution functions 32,33 (cf. Fig. 2). At variance, the LGS model yields a temperature dependent structure for e − aq , which causes E max to redshift by 5.2 meV/K in constant-density simulations (blue circles in Fig. 3(a)), 33,45 in qualitative agreement with the measured slope of 2.2 meV/K. 31 However, increasing evidence supporting cavity formation [25][26][27][28][29] disfavors the LGS model, in which the non-cavity structure results from an overestimated electronwater interaction. 46</p><!><p>In this work, we provide a comprehensive theoretical study of the temperature dependent properties of e − aq by combining accurate electronic structure calculations at the hybrid-functional level, machine-learning potentials, and a multiple-time-step scheme to achieve long time-scale molecular dynamics simulations. The achieved configurations are then used in a QM/MM set-up to obtain the absorption spectrum through TDDFT. We develop two machine-learning potentials for describing neat water and e − aq in various thermodynamic conditions. In this way, we first simulate liquid water at a fixed pressure of 250 bar and at various temperatures, thereby providing a good description of the temperature-dependent density of water between 300 and 573 K. Starting from the equilibrium trajectory of liquid water at various temperatures, we then address e − aq within our accelerated machine learning scheme yielding equilibrated temperature-dependent properties. We find that the structure of e − aq is strongly affected by temperature but is less sensitive to density. Our results demonstrate that cavity-forming electronic structure theory predicts E max accurately under a variety of thermodynamic conditions. In particular, we reveal that higher temperature leads to larger gyration radii of e − aq , thereby causing the redshift of E max . At last, we expect that the computational framework used here can efficiently be used for modeling excess electrons and holes in condensed phases.</p><!><p>Our methodology sequentially combines various advanced simulations. In the first step, we separately simulate neat water and e − aq using AIMD under various thermodynamic conditions (see Fig. S2). We use heavy water to match the conditions in the pulse radiolysis experiments of Bartels et al. 30 We do not include nuclear quantum effects, because they have been shown to not affect the structure of e − aq in heavy water at room temperature 28 and because they are expected to be even less significant at higher temperatures. 47 The underlying electronic structure is described with the spin-polarized hybrid functional PBEh(0.40), in which the fraction of Fock exchange is set to 0.40. The functional is supplemented with van der Waals interactions through the nonlocal rVV10 functional, 48 where the parameter b is set to 5.3 to ensure that the structural and dynamical properties of liquid water are properly reproduced. 49 The PBEh(0.40)-rVV10 functional provides an accurate description of the electronic structure of e − aq , comparable to higher-level results from perturbation theory. 26,28 In particular, this set-up yields a Kohn-Sham band gap of 8.8 eV for liquid water, 49 in agreement with recent bounds set for the experimental value. 50 All the periodic AIMD calculations are performed with the CP2K code 51 employing TZV2P basis sets 52 and Goedecker-Teter-Hutter pseudopotentials. 53,54 The exchange integrals are calculated with the auxiliary density matrix method. 55 In the second step, we use the data obtained in the first step to train two MLPs, denoted MLP(water) for neat water and MLP(e − aq ) for e − aq . We adopt the neural network scheme introduced by Behler and Parrinello, 38 which is implemented in the n2p2 code. 56 The MLP is further improved by running MLP-based molecular dynamics in different thermodynamics conditions.</p><p>Hybrid-functional energies and forces are calculated for representative configurations selected on the basis of their atomic fingerprints 57 and used for extending the initial data set. Cross-validation of the MLPs achieved in this way is illustrated in Fig. S1 of the Supplementary Information (SI).</p><p>In the third step, we simulate liquid water with the MLP(water) under various thermodynamic conditions. 30 All simulations of e − aq start from configurations of neat water as equilibrated within the canonical ensemble at experimental density and temperature (NVT), relying on the fact that MLP(water) gives a satisfying density for neat water (Fig. S2). For e − aq , we applied a multiple time step (MTS) method. The full force based on the PBEh(0.40)-rVV10 functional is evalu-ated every 2.0 fs, corresponding to four time steps evolved with the reference force based on the MLP(e − aq ). We perform two independent simulations lasting 10 ps for each thermodynamic condition. The MTS-AIMD simulations are performed with the i-PI code 40,58 interfaced with CP2K for the electronic structure calculations and LAMMPS 59 for the MLP. The e − aq is characterized by the spin-density distributions, the center of which is calculated from the second moment tensor (see SI).</p><p>In the final step, we carry out TDDFT calculations for configurations of e − aq from the MTS-AIMD simulations using a QM/MM protocol as implemented in the ORCA code. 60 For each adsorption spectrum in the various thermodynamic conditions, we use 40 configurations homogeneously distributed in time between 2 and 10 ps after the addition of the electron. The QM region comprises 70 water molecules closest to the center of e − aq , corresponding to the cluster with the largest radius fully contained in the simulation cells. The MM region is described through a point charge model and constructed on the basis of 560 surrounding water molecules in the periodic supercells. The rest of the environment is described by the conductor-like polarizable continuum model 61 with the experimental dielectric constant of water at the given thermodynamic condition. 62 In the TDDFT calculations, we use the PBEh(0.40) functional with ma-def2-TZVP basis set 63 and consider the first five excited states. The QM/MM set-up is illustrated in Fig. S4.</p>