Tracking Cavity Formation in Electron Solvation: Insights from X-ray Spectroscopy and Theory

We present time-resolved X-ray absorption spectra of ionized liquid water and demonstrate that OH radicals, H3O+ ions, and solvated electrons all leave distinct X-ray-spectroscopic signatures. Particularly, this allows us to characterize the electron solvation process through a tool that focuses on the electronic response of oxygen atoms in the immediate vicinity of a solvated electron. Our experimental results, supported by ab initio calculations, confirm the formation of a cavity in which the solvated electron is trapped. We show that the solvation dynamics are governed by the magnitude of the random structural fluctuations present in water. As a consequence, the solvation time is highly sensitive to temperature and to the specific way the electron is injected into water.

The optical ionization pump-soft x-ray probe measurements were performed at the Chem-RIXS endstation of the Linac Coherent Light Source (LCLS) x-ray free-electron laser (XFEL).Fig. S1 shows a schematic of the experimental setup.The 800 nm output of a Ti:sapphire amplifier was focused onto a 700 nm-thick converging liquid sheet jet target 1,2 to ionize liquid water by strong-field ionization.The surface of the liquid jet was oriented normal to the incident 800 nm ionization pump and soft x-ray probe beams, with a crossing angle of ∼1°S between the two beams.The pump laser was combined with the x-rays by a holey mirror upstream of the target and skimmed out by another holey mirror downstream of the target.
With a pulse duration of 55 fs, a focal spot diameter (FWHM) of 75 µm and a pulse energy of 0.5 mJ, the peak intensity of the pump laser was 1.3 × 10 14 W/cm 2 at the interaction region.
The intensity and polarization of the 800 nm pump laser could be controlled by an ND filter and a λ/2 waveplate in the 800 nm delivery path.The 20 fs XFEL pulse becomes stretched to ∼ 30 fs after being monochromatized to a spectral width of 0.1 eV.It was then focused by a pair of Kirkpatrick-Baez (KB) mirrors to a spot size (FWHM) of 32(H) × 130(V) µm 2 on the target.The focal point was located between the arrival time monitor (ATM) 3 and the target interaction point (IP), to provide sufficient x-ray intensity on the ATM screen.The ATM, as shown in Fig. S2, positioned 1.5 m upstream of the target, cut off about half of the x-ray beam before the focus and recorded the relative arrival times of the optical laser and the XFEL pulses.For the experimental results shown in Fig. 1 in the main manuscript, the optical-XFEL timing jitter limited the effective time resolution to ∼65 fs.The time delay between the optical pump and x-ray probe pulses at the interaction region was controlled by the upstream global delay on the optical arm.The local optical delay stage upstream of the ATM (shown in Fig. S1) is employed to fulfill the operational requirements of the ATM.

S-I.1 Incoming energy calibration and normalization
The incident energy (I 0 ) of the monochromatic x-ray pulse was recorded by fluorescence intensity monitors (FIM) installed above the KB mirrors. 4The transmitted energy (I) of x-ray pulse was recorded by a camera (Andor Newton SO DO940P) in full-vertical-binning mode, where the 2D image was projected to a vertical 1D array.The 1D array maps the transmitted x-ray intensity from various points of the jet leaf from upstream to downstream in the interaction region, which makes it possible to extract the x-ray absorption spectrum (XAS) for various thicknesses and various pump laser intensities.Based on the Beer-Lambert law ∆A = log(I 0 /I), the XAS of the water sheet, shown in Fig. S4, can be easily obtained by subtracting the no-sample XAS from the with-sample XAS.The correlation between I 0 and I without sample is shown in Fig. S3.The XAS is recorded by scanning the x-ray photon energy with a step size of 0.2 eV, with more than 1000 shots per step.The x-ray photon energy scan is performed by scanning the undulator gap in coordination with the monochromator; this gave about 50% variation of the averaged x-ray pulse power over a S-3 relatively large photon-energy range, 520 to 550 eV.The energy of the monochromator readout was calibrated by the pre-edge feature of liquid water absorption.

S-I.2 XAS of liquid water
The water XAS shown in Fig. S3 was taken by a synchronized undulator-monochromator scan near the oxygen K-edge region.The background absorption was taken without the water jet.Fig. S4 shows the experimental results here reproduce the reference water XAS spectrum obtained by Meibohm et al.

S-I.3 Transient XAS of ionized water at 2 ps delay
The transient x-ray absorption spectrum ∆A of strong-field ionized water was recorded with a pump-probe delay at 2 ps.A trigger kicks out the pump laser 30 ns after the x-ray probe pulse and was used to turn the pump laser 'off'.The ratio of the number of shots with pumpon and with pump-off was 1:1.Fig. S5 shows the x-ray/optical laser overlap -together with the averaged x-ray transmission signal as a function of X-ray camera pixel number.The center position of the laser projected to the camera is determined by the measured ∆A signal.The relative widths of the x-ray and laser were determined by knife-edge scans.By selecting the pixel number of the transmitted x-ray signal on the camera, we select the ∆A for specific pump laser intensities.The purple region corresponds to a pump laser intensity of ∼ 10 13 W /cm 2 , which is ∼ 7% of the peak pump intensity and similar to that used in earlier study of strong-field-ionized liquid water. 6This yields an ionization fraction of ∼ 1%.We note that the intensity required for a given ionization fraction is inversely proportional to the target thickness.
Figure S5: The spatial overlap between X-ray and laser.Green dash-dot line shows the averaged X-ray transmission intensity on the X-ray camera, from which the X-ray profile was fitted as shown by the green solid line.The integral of the absolute ∆A over 525 − 540 eV region on each X-ray camera pixel was used to fit the center position of the pump laser.The spot sizes were measured by knife-edge scans and the Gaussian laser profile was reconstructed as shown by the red solid line.The transient absorption spectrum at 2 ps shown in Fig. 1B in the main manuscript is from the purple region.

S-I.4 Liquid Jet
The thin water jet avoids saturation of absorption by the target.The averaged thickness over x-ray spot is estimated to be ∼ 700 nm as deduced by comparison with the absolute absorbance of Ref. 7 at the water main edge.As shown in Fig. S6(a), the thickness of the water jet is also measured through the transmission thin-film interferometry 8 at Argonne National Laboratory.The white light from a white light lamp (Thorlabs, OSL2 Fiber Illuminator) is focused on the water sheet and the transmitted spectrum is measured with a spectrometer (Ocean optics, USB2000).Here the flow rate of the HPLC pump was set to 1.75 ml/min, the pressure was 275 psi and the thickness of water sheet as a function of the distance from the nozzle is shown in Fig. S6(c).Due to the node at the bottom of the first water sheet leaf, the thickness there shows a rapid increase.Fig. S6(b) shows an example of fitting the water-sheet thickness from the transmission of the white light spectrum, in the S-6 following fitting function 8 t and r are the transmission and reflection of water, respectively, n is the refractive index of water, d is the thickness of water sheet.The parameters of a, c, k are added to get better fitting which are not sensitive to the oscillation period and therefore do not affect the measured thickness.

S-II.1 Dynamics and characterization of electron solvation
The modelling of the solvation dynamics is based on a force field implemented via a neural network (through the n2p2 package 9 ).This neural-network force field was trained at the restricted/unrestricted MP2 level of the electronic ground state of neutral/anionic water with a triple-zeta quality correlation-consistent basis set as implemented in CP2K. 10 For more information on the training of the machine learning potential (MLP), see the SI from Ref. 11.The RPMD simulations were performed using 47 water molecules, with 32 beads for each atom, in a box with an edge length of 11.3 Å and periodic boundary conditions.The propagation was performed using the i-PI code 12 with a 0.25 fs time step.With a similar approach as the one used by Lan et al., we begin by launching a simulation of neutral water for 5 ps at a given temperature (300 K or 340 K).Once these calculations with neutral water finished, the water box was considered to be equilibrated.The last snapshot of each of these trajectories was used as the starting point for a trajectory of the electron solvation dynamics.
To introduce the excess electron in the simulation, we simply had to switch to the neuralnetwork force field trained for the anionic water box and run the calculation for another 5 ps.
We launched three sets of trajectories of the electron solvation dynamics, one at 300 K with a thermostat, one at 340 K with a thermostat, and another at 300 K with no thermostat.
For the thermostated calculations, we employed an SVR thermostat (all parameters for the thermostat can be found in the aforementioned Ref. 11).We launched between 34 and 37 trajectories for the three different sets of calculations.During the dynamics, some of these trajectories failed because they explored configurations for which the force field was not properly trained, so we had to discard them.In the end, we obtained 27 viable trajectories for the calculations at 300 K with a thermostat, 25 for the calculations at 340 K with a thermostat, and 26 for the calculations at 300 K without a thermostat.All files necessary to replicate these calculations, such as the i-PI input files and the files containing the parametrization of the neural-network force field, can be also found in the SI of Ref. 11.

S-II.1.1 Detection of Cavities
To identify the formation of cavities, we followed the same analysis procedure as Lan et al. 11 That is, from selected snapshots of the RPMD trajectories, we calculated the spin density using density functional theory with the CP2K code. 10 We employed the PBE(α) exchange-correlation functional with 40% Hartree-Fock exchange, which has been shown to provide a satisfactory description of the solvated electron, particularly in terms of producing a spin density comparable to that obtained using MP2, see SM from Ref. 11.Because these calculations are very demanding, we calculated the spin density only for snapshots of the centroid at every 100 fs.We observe that the spin density (and thereby the density of the excess electron) is initially completely delocalized throughout our simulation box.Only at some later time, the spin density localizes somewhere, which is accompanied by the formation of a cavity around the solvated electron.An example of such a formed cavity can be seen in Fig. S7.
Once we obtained the spin density from the DFT calculation, we shifted the periodic boundary conditions such that the spin density is minimal on the surfaces of the simulation box.We then computed the center of the spin density ρ s (r), Figure S7: Spin density of the solvated electron for an isovalue of 0.0008.Water molecules represented with spheres belong to the first solvation layer.As can be seen, spin density can be found not only in the central excluded volume (the cavity), but also on the oxygen atoms of the first and even second solvation layer.
and the gyration radius by diagonalizing the gyration tensor where m and n refer to the different x, y, z Cartesian components.We consider that a cavity has formed when the gyration radius is less than 3 Å.From this point, we can identify all water molecules that belong to the cavity of the solvated electron.To do so, we construct the isosurface of the spin density for an isovalue of ρ s = 0.001.Generally, as soon as the electron has localized, such an isosurface encloses one large volume and several smaller volumes, i.e., spin bubbles.We focus on the isosurface associated with the largest volume and identify water molecules within a certain distance from this isosurface as belonging to the cavity.

S-9
Figure S7 illustrates a typical example, where one can see five water molecules that are within 1.4 Å from the isosurface of the largest spin bubble.For illustration purposes the depicted isosurface has a somewhat smaller isovalue of 0.0008.
The results for the gyration radius from all trajectories as a function of time are shown in Fig. S8.We have included the results of a global fit to the three sets of calculations to an exponential law where the baseline, R g (∞), represents the asymptotic gyration radius of each set of trajectories.For the calculation at 300 K and without thermostat, a simple monoexponential fit yielded considerably worse results at early times compared to the fits for the thermostated calculations (see bottom panel of Fig. S8).However, using a biexponential fit produced better results.Both fits are included in the lower panel of Fig. S8 S-II.1.2Characteristic Formation Times: Fitting Procedures Given that the initial delocalization of the electron in our simulations depends on the box size, potentially influencing the solvation timescale obtained through gyration radius fitting, we performed an additional analysis for a more in-depth analysis of the characteristic times of the solvation process.For each trajectory, we adopted a binary criterion and assigned a value for each of the two possible states of the cavity at any given time, formed or not formed.
A cavity is considered to be formed if the gyration radius at a given time is less than 3 Å.We assign a value of 1 to trajectories without a cavity (R g > 3 Å), and a value of 0 to those with a cavity (R g < 3 Å).The average of all trajectories corresponds to the fraction of trajectories in which, at a given time, cavity formation has not yet taken place.The results, together with a global fit to an exponential function, are shown in Fig. S9.As in the previous case, for the nonthermostated calculation, we obtained a better fit using a biexponential function rather than a monoexponential function.In this case, the improvement of the fitting is even more apparent.All global fits shown in Figure S8 and Figure S9 were performed using the lmfit library 13 together with bootstrapping to estimate standard error of the fit parameters.Figure S8: Gyration radius of the spin density as a function of time for the RPMD trajectories (colored thin lines) calculated at 300 K with a thermostat (top-left panel), 340 K with a thermostat (top-right panel), and 300 K with no thermostat (bottom panel).Black and red lines correspond to the mean gyration radius and its fit to a monoexponential decay law, respectively.The shaded area represents the standard error due to sampling.The results for the case of 300 K with no thermostat (bottom panel) also include a biexponential fit in blue.

S-10
S-11  Black and red lines correspond to the mean gyration radius and its fit to a monoexponential decay law, respectively.The shaded area represents the standard error due to sampling.The results for the case of 300 K with no thermostat (bottom panel) also include a biexponential fit in blue.
the respective oxygen atoms.The description of the final, core-excited electronic state |ψ j ⟩ having N electrons is based on molecular orbitals optimized for the core-ionized (N − 1)electron configuration.These orbitals are further projected into a larger basis set where additional 11 s, 11 p, and 11 d-type Kaufmann basis functions 28 have been added on the respective core-ionized oxygen atom.The respective Kaufmann-basis functions are optimized for the specific ionized final charge state.For example, for the absorption by the hydrated hydronium cation (H 3 O + ), Kaufmann-basis functions were taken that are generated for a dicationic continuum.Using this larger basis set, we describe the final electronic states |ψ j ⟩ by diagonalizing the electronic Hamiltonian in the space of N -electron configurations, in which an electron has been added to a virtual or singly occupied orbital of the core-ionized configuration.For the cavity water molecules around a solvated electron, the final-state expansion is augmented by further electronic configurations where the solvated electron can occupy the n lowest unoccupied orbitals, where n is the number of water molecules included in the electronic structure calculation.This provides more flexibility in the final configuration expansion, accounting for the fact that the rather diffuse solvated electron orbital adapts quite significantly in the presence of the core hole and shows rather small overlap in the two orbital sets.For the calculations involving OH • , configurational mixing in the two almost degenerate π hole orbitals is considered in the initial and final states as well.
The transition dipole operator in Eq. ( 5) is then calculated using Löwdin's generalized expression of transition matrices.The two spectral parts, the bound-bound and the bound-continuum (obtained via the Stieltjes imaging procedure), are then merged at an energy value where both spectral parts overlap.

S-II.3.1 Sampling of structures for XAS calculations
From snapshots of the PBE0+D3 MD simulations, we calculated the XAS, for x-ray absorption of the OH radical, the H 3 O + cation, water molecules neighboring H 3 O + , and water molecules in bulk water.To that end, the local environment of the absorbing species was considered in the calculation by incorporating the absorbing molecule and the closest water molecules in the electronic structure calculations for the x-ray absorption (QM region).The remaining water molecules up to a distance of 10.6 Å were incorporated as point charges, using values from the SPC/E water force field. 33The selection was based on the oxygen-oxygen distance.
The individual criterion on how the geometrical structures for the QM region were selected is described in the following.
• Bulk water (MD): We picked snapshots of the water structure at two time frames, namely at 7.5 ps and at 10 ps of simulation time.From both snapshots, individual x-ray absorption spectra calculations were performed considering every water molecule as the absorbing water molecule.For each calculations the 9 closest water molecules were incorporated into the QM region.
• Bulk water (RPMD): We picked a single snapshot of the water structure (centroid po-

Figure S1 :
Figure S1: Schematic layout of the experimental setup at the ChemRIXS beamline: 800 nm strong-field ionization pump, time-delayed x-ray probe, arrival time monitor and transmission x-ray camera.

Figure S2 :
Figure S2: Optical layout for arrival time monitor (ATM) and XAS measurement.The x-ray was partially cut off by the ATM screen and focused upstream of the water jet; the transmitted x rays were recorded by the x-ray camera.The x-ray intensity distribution from the x-ray camera is shown on the right.The relative sizes of the x-ray (purple) and laser spots (red) on the water jet is shown in the dashed circle.

Figure S3 :
Figure S3: The XAS with and without water sample.The arbitrary vertical scale is offset from zero in the water window because different types of detectors were used to measure the incident and transmitted x-ray intensity.

Figure S4 :
Figure S4: The XAS of water sample.Black dot is experimental data, red curve is the XAS spectrum at 23 • C.5

5
Figure S4: The XAS of water sample.Black dot is experimental data, red curve is the XAS spectrum at 23 • C.5

Figure S6 :
Figure S6: Measurement of the water sheet jet thickness.(a) Schematic for the whitelight interferometry, (b) Typical transmission recorded by the spectrometer, (c) Water sheet thickness fitted by the equation 1, as a function of the distance to the nozzle.

Figure S10 :
Figure S10: Evolution of the solvation process as a function of the relative formation time ∆τ = t − τ , where τ represents the cavity-formation time of each trajectory.Gyration radius of the electron spin density, R g (∆τ ), for the various sampled trajectories for the simulations without a thermostat (left) and the simulations at 340 K.
31By taking two different sets of orbitals, orbital relaxation effects are directly taken into account and there is no imbalance by choosing a common orbital set adequate for a specific state (e.g., either the initial or the final state).The resulting absorption resonance energies ϵ ij are thus considerably more accurate compared to resonance energies evaluated using a common orbital set.Compared to experimental data there remains an absolute energy shift of 1 eV-2 eV due to relativistic effects in the oxygen 1s shell (relativistic effects are about 0.8 eV for core ionization in neon31), limited basis set effects, and further many-electron correlation effects not considered in the calculation.To test our procedure, we inspected the summed oscillator strength for the considered absorption transitions and found values that are on average 0.93 (for H 2 O), 0.95 (for OH • ), 1.01 (for H 2 O neighboring H 3 O + ), 1.25 (for H 3 O + ), and 0.79 (for H 2 O in a solvated electron cavity). Th low value for a cavity water molecule is connected with the mentioned diffuse nature of the solvated electron that responds particularly sensitively to the creation of the core hole.