Probing aqueous ions with non-local Auger relaxation

Abstract

Non-local analogues of Auger decay are increasingly recognized as important relaxation processes in the condensed phase. Here, we explore non-local autoionization, specifically Intermolecular Coulombic Decay (ICD), of a series of aqueous-phase isoelectronic cations following 1s core-level ionization. In particular, we focus on Na⁺, Mg²⁺, and Al³⁺ ions. We unambiguously identify the ICD contribution to the K-edge Auger spectrum. The different strength of the ion–water interactions is manifested by varying intensities of the respective signals: the ICD signal intensity is greatest for the Al³⁺ case, weaker for Mg²⁺, and absent for weakly-solvent-bound Na⁺. With the assistance of ab initio calculations and molecular dynamics simulations, we provide a microscopic understanding of the non-local decay processes. We assign the ICD signals to decay processes ending in two-hole states, delocalized between the central ion and neighbouring water. Importantly, these processes are shown to be highly selective with respect to the promoted water solvent ionization channels. Furthermore, using a core-hole-clock analysis, the associated ICD timescales are estimated to be around 76 fs for Mg²⁺ and 34 fs for Al³⁺. Building on these results, we argue that Auger and ICD spectroscopy represents a unique tool for the exploration of intra- and inter-molecular structure in the liquid phase, simultaneously providing both structural and electronic information.

---

1 Introduction

The recent exploration of the molecular structure in soft condensed matter and liquids forced a revision of a number of previously helpful paradigms and revealed novel, hitherto unobserved effects.^1–3 In the context of the interaction of radiation and liquids, the discovery of Intermolecular Coulombic Decay^4–6 (ICD) is probably one of the most exciting phenomena. Upon first glance, ICD shows similarities with the well-known Auger decay, however, it is a decay mechanism characteristic of condensed systems only. It takes place in an atom or molecule after an inner-valence or inner-shell vacancy is created. The vacancy is refilled with a valence electron and the released energy ionizes a neighbouring atom or molecule. The whole process leads to an energetically favourable delocalized two-hole state {Center⁺⋯Surrounding⁺} in contrast to an Auger decay, where a higher energy {Center²⁺} state is formed.

ICD and other non-local Auger-type processes have entered the field of X-ray science relatively recently and have primarily been discussed in the context of their role in radiation chemistry or X-ray photochemistry.^7–9 The simultaneous ionization of neighbouring species leads to completely new reaction channels that should be considered within the overall radiolysis mechanisms.¹⁰ Non-local Auger phenomena also offer new possibilities to probe liquid-phase molecular structure, as the decay processes are highly dependent on the intermolecular distances. Some of us have already suggested that the analysis of electron spectra after autoionization of an ionic metal centre could be interpreted with respect to the environment of the respective ion.^11,12 Moreover, in the liquid state, ICD spectra of metals involving surrounding water molecules have already been observed.^6,13,14 In fact, the sensitivity of electronic relaxation spectra to the chemical environment of the emitter started to be discussed immediately after the discovery of ICD.^15,16 Much like Förster energy transfer, ICD also exhibits a 1/R⁶ dependence on the intermolecular distance, R (ref. 6). The impact of ICD, however, is much broader than Förster transfer, since every atom or molecule may act as the receiving end of an inter-centre energy transfer, leading to the ejection of a secondary electron. Corresponding X-ray spectroscopies thus represent a relatively new tool to reveal liquid structure, complementing more conventional techniques such as dielectric spectroscopy or neutron scattering. Previously, the latter methods have been used in conjunction with molecular dynamics (MD) simulations to infer the arrangement of anions and cations in an electrolyte solution.¹⁷ X-ray-based spectroscopies and associated theoretical treatments represent an alternative to these well-established techniques due to their ability to directly and atom-specifically probe the electronic structure of the sample, while simultaneously delivering comparable information about the geometrical microstructure. With the combination of bright high-energy light sources, accurate electron analyzers, and liquid-microjet approaches, we can now relate the liquid structure to its signature in electron spectroscopy.¹⁸ Hand-in-hand with advances in X-ray-based experimental techniques for probing liquids' electronic structure, also the methods to simulate such structures have advanced immensely.^19–23

Despite the enormous interest in the non-local Auger-like decay processes in the last decade, their application to elucidate liquid-phase molecular structure remains scarce. In the cases probed so far, these non-local signatures were weak, and since LVV Auger decay was considered, they consisted of a convolution of the water and metal valence shells.^13,24 This complicated the interpretation of the spectral signatures. In the present work, we focus on the ICD signatures in the K-shell spectra. Such measurements provide a clear-cut case: the signal should be dominated by a valence vacancy of the surrounding and an L-shell vacancy in the metal centre, with a fairly well-defined energy.

We present experimental K-shell Auger and ICD spectra of aqueous-phase Al, Mg, and Na ions, associated with AlCl₃, MgCl₂, and NaCl solutions. These cations were selected because they are isoelectronic but represent different types of interaction with neighbouring molecules. While the sodium cation is only weakly bound to neighbouring water molecules, dicationic magnesium is more strongly coordinated and the aluminium ion forms a regular coordination-covalent bond; the ICD signal is shown to reflect this diversity. The measurements are accompanied by molecular dynamics (MD) simulations and ab initio calculations, which aid in the interpretation of the experimental data and allow for an in-depth analysis of the spectra.

2 Methods

2.1 Experimental methods

The Auger and ICD electron spectra were measured using the EASI photoemission setup,²⁵ equipped with a liquid microjet and a hemispherical electron energy analyzer, at the P04 beamline of the synchrotron radiation facility PETRA III, DESY, Hamburg.²⁶ The beamline has a high on-target photon flux of about 2 × 10¹² photons s⁻¹ at a resolving power of 10 000, with a photon-energy tuning range spanning 250–3000 eV, and variable circular polarization.²⁷ The beamline's 1200 lines per mm grating yields a photon-energy resolution of 250 meV at 1200 eV photon energy and 350 meV at 1500 eV, using an exit-slit opening of 100 μm. At these settings, the vertical spot size amounts to approximately 35 μm. The general properties of the liquid-microjet system are described elsewhere.^25,28 The liquid microjet of the sample solutions was introduced into the vacuum chamber using an HPLC pump at a flow rate of 0.8 ml min⁻¹ with a backing pressure of ≤ 12 bar, and was directed horizontally. The glass capillary nozzle used to introduce the sample into the chamber had an inner diameter of 28 μm. The synchrotron radiation was incident perpendicular to the flow of the solution. A near-ambient-pressure hemispherical electron analyzer (Scienta Omicron HiPP-3), mounted at a 50° backward-scattering angle with respect to the beamline (near magic angle), was used to measure the electron kinetic energy.²⁵ The solutions were prepared by dissolving commercially purchased AlCl₃, MgCl₂, and NaCl (Sigma-Aldrich with purity > 98%) salts in MilliQ (18.2 MΩ cm⁻¹) water. Aqueous solutions of MgCl₂ and NaCl had a concentration of 1 M, while that for AlCl₃ was 2 M. At pH ≤ 4, the aluminium is predominantly present as the aluminium hexahydrate cation ([Al(H₂O)₆]³⁺), i.e., the Al³⁺ surrounded by six water molecules.²⁹ The utilised aluminium chloride solutions were highly acidic (pH < 2), which results in the exclusion of large amounts of chloride from direct contact with Al³⁺.

The kinetic energy of the electrons produced in the decay processes, both local (Auger) and non-local (ICD), are independent of photon energy (hν). Therefore, to distinguish the decay features from the photoelectron peaks, electron spectra were measured for two photon energies differing by 3 eV in all cases. Hemispherical-electron-analyzer pass energies of 100 or 200 eV were used together with an analyzer slit width of 800 μm, leading to an estimated analyzer resolution between 0.2 and 0.4 eV.

Photon energies and kinetic energies were calibrated as detailed in the ESI,† in part making use of additional measurements carried out at the U49-2_PGM-1 beamline of the BESSY II synchrotron-radiation source at the Helmholtz–Zentrum Berlin für Materialien und Energie, using the SOL³PES setup for liquid-jet photoemission spectroscopy.^28,30

The observed ICD features were interpreted by a line-shape analysis using Voigt profiles in the SPANCF (Spectrum Analysis by Curve Fitting) macro package³¹ for Igor Pro (Wavemetrics, Inc., Lake Oswego, USA). In our case, either a 2p or 2s electron of the cation recombines with the 1s hole following direct photoemission, and the released energy leads to the secondary emission of electrons from the water valence states. The multipeak structure of each individual ICD feature (see Fig. 2) thus arises from the release of electrons from different valence states. During the data fitting process, the kinetic energy, the intensity of the peaks, and the Gaussian width were free to vary and the Lorentzian width was considered to be that of the cation 2p orbital (see Fig. 4). The relative kinetic energy differences of these peaks are almost the same as the relative binding energies of the water valence-band peaks.³² The main KLL Auger peak notably has a higher-kinetic-energy asymmetric tail that could be modelled by a PCI (post-collision interaction) profile, as originally developed for gas-phase work.

Fig. 1 Scheme of the two QM/MMPol models employed in the present study. The QM part contained one metal cation and six water molecules or one cation and one water molecule (minimal model), the remainder of the 20 Å-radius sphere was treated at the polarizable embedding MM level.

Fig. 2 Photoemission spectra of aqueous solutions of AlCl₃, MgCl₂, and NaCl (panels a–c). Two photon energies were selected and implemented slightly above the respective K-edge of the metal and the resulting spectra are compared. In each panel, the left-hand side shows the main (KLL) Auger peak, and the right-hand side photoemission peaks due to direct L-shell ionization and due to K-shell ICD involving an L-shell electron (see labels). In the case of Na, no ICD peak can be observed, and the water valence band extends into the observed spectral region. Right-hand side spectra were scaled to improve visibility. See text for details.

Fig. 3 The schematic energy-level diagram of the relaxation of a 1s core hole in aqueous metal (M) ions. The local decay channel (KLL Auger decay) and non-local decay (ICD) are shown. The neighbouring water molecules are indicated by ‘w’.

Fig. 4 Magnified view of the L-shell photoelectron and ICD peaks of AlCl₃ (upper panel) and MgCl₂ (lower panel) aqueous solutions shown in Fig. 2 (symbols). Spectra were measured with hν = 1569.8 eV and 1315.25 eV, respectively. A decomposition of the ICD features into components pertaining to final states of different energy, modelled by Voigt profiles and found by least squares curve fitting is shown by coloured traces (see text for details). Subcomponents of the ICD structure correspond to different water valence-band orbitals, namely 1b₁ (yellow), 3a₁ (green), 1b₂ (blue), and 2a₁ (purple) and appear for both ICD_2p and ICD_2s. The peak fits to the two spectra measured with higher hν are given in Fig. S2 of the ESI.†

2.2 Theoretical methods

Molecular dynamics. Classical molecular dynamics (MD) simulations were performed to sample the configurational space. MD simulations were performed for 1 M NaCl and MgCl₂, and 2 M AlCl₃ solutions to match the experimental conditions. The classical non-polarizable force fields were employed to generate large-scale structural snapshots for subsequent QM/MMPol calculations. The details of the classical MD simulations are summarized in the ESI.†

Ab initio calculations. For the ab initio calculations, the systems were divided into relatively small quantum (QM) and extensive molecular mechanics (MM) parts. In this work, we employed the QM/MMPol embedding scheme^33,34 with the AMOEBA polarizable force field.^35,36 In the polarizable force fields, each atom is described by a static point charge and an isotropic polarizability; more sophisticated force fields such as AMOEBA also contain higher multipoles. Recently, Lipparini et al.^37,38 have introduced an efficient implementation for polarizable QM/MMPol based on the Fast Multipole Method,³⁹ which makes the calculation of even very large systems affordable. In our case, the MMPol part was a 20 Å-radius sphere of the respective solution surrounding the QM part. All structures were taken from classical molecular dynamics calculations. The binding energies of the aluminium cation and of a chloride anion were also checked for a 30 Å-radius sphere; the calculated values were within the error bars for smaller systems, i.e., the values for a 20 Å sphere can be considered converged. The present model acknowledges the granularity of the solvent around the solute, takes into account the effect of varying ionic strength of the solutions, and allows us to control the convergence of the energetics with the increasing size of the simulation box. Importantly, the model provides an optimal treatment of electronic polarization which is instrumental for a correct description of the ionization energetics.

For the QM part, we considered two models – (1) a minimal model containing a single cation and a single water molecule (or chloride anion), (2) the first solvation-shell model containing one cation and six water molecules, which corresponds to the water coordination numbers for the investigated cations.^40–46 The considered QM/MMPol systems are shown in Fig. 1.

The core-level energies, binding energies (BEs) and energies of the lowest two-hole final states (of a triplet multiplicity) were calculated by the Maximum Overlap Method (MOM).⁴⁷ This approach allows the variational convergence of states with specifically localized hole(s) with any ground-state method. In this work, we performed the MOM calculations at the LC-ω PBE level with the range-separated parameter set to a default value of 0.4 bohr⁻¹ with the cc-pCVTZ basis set for a cation and the cc-PVTZ basis set for all other atoms. We did not tune the ω paratemer for particular systems or geometry since we mainly aim at pointing out the relative differences in spectra and energies.

Yet, the method has some limits, in clusters with a high density of energetically close-lying electronic states, the convergence of the MOM method can be poor. The BEs and energies of the final two-hole states were calculated for a set of 20 geometries selected from the classical molecular-dynamics simulations. The BEs were calculated as an energy difference of the ground and singly ionized states, the energies of the final two-hole states were evaluated relative to the ground state.

The valence photoemission spectra for cations and their first solvation shell were also calculated with a recently introduced ionization-as-an-excitation-into-a-distant-center (IEDC) approach.^23,48,49 The method is based on modelling the ionization from a selected orbital space as an excitation into a continuum using time-dependent density functional theory (TDDFT), similarly as was previously suggested by Stanton and Gauß⁵⁰ as well as Coriani and Koch.⁵¹ As DFT is in principle an exact many-body theory, we can obtain correlated orbital energies. Similar to the MOM method, the IEDC approach was performed in the QM/MMPol arrangement at the same level of theory with a sodium cation as a distant centre (placed at a distance of 1000 Å from the system). Because in the QM/MMPol model the distant centre was not solvated, we shifted the excitation energies so that the lowest-energy TDDFT transition agrees with the first ionization energy of the cluster in the MMPol embedding scheme. In this setting, the choice of a sodium cation as a distant centre is arbitrary. A 200-frame set of 20 Å spheres was cut from the classical molecular dynamics, the QM part contained one cation and six neighbouring water molecules. Excitation energies were calculated at the LC-ω PBE/cc-pVTZ level; the ω parameter was set to 0.4 bohr⁻¹.

The population analysis was performed for [Al(H₂O)₆]³⁺ and [Mg(H₂O)₆]²⁺ complexes optimized at the BH&HLYP 6-31 + g* level in the polarizable continuum. The Löwdin reduced orbital population per molecular orbital was performed in ORCA 4.2.0.⁵² All other ab initio calculations were performed using the locally modified current development version of Gaussian 16.⁵³ All classical MD simulations were performed with the GROMACS 5.1.2 package.^54,55

3 Results and discussions

Experimental manifestation of ICD

The K-shell Auger and ICD spectra investigated in this work are produced by 1s photoionization of a metal cation in the respective electrolyte. Inner-shell photoelectron spectra resulting from this primary process were measured to determine the binding energies of the respective core levels, and are shown in Fig. S3 of the ESI.† The secondary-electron spectra measured above the metal ion 1s ionization thresholds of the NaCl, MgCl₂, and AlCl₃ aqueous solutions are the main subject of this article. In order to maximise their intensity, photon energies near the respective K-shell thresholds were chosen to measure the spectra shown in Fig. 2. The secondary-electron spectra associated with each sample were measured with two different photon energies to distinguish peaks from different processes: The kinetic energies of the primary photoelectron peaks disperse with photon energy, whereas the kinetic energies of the Auger and ICD features remain constant. For all three ions, the KLL Auger spectra are shown in the left panels and 2s⁻¹ and 2p⁻¹ photoelectron peaks are seen in the right panels of Fig. 2. In the case of Mg and Na solutions, a different baseline height (resulting from scattered electrons) was found at the two photon energies probed, probably because spectra were recorded in the proximity of the respective K-shell thresholds. To enable comparison of the ICD features, spectra were scaled and shifted on the intensity scale appropriately, which leads to most or all of the apparent difference in Mg 2p and 2s intensity for the two photon energies shown in Fig. 2b. For Na⁺, the water valence-band peaks (w⁻¹) are also in the selected kinetic energy range. Note that the Na⁺ 2p peak coincides in energy with one of the water valence-band peaks, 2a₁⁻¹. The peaks with constant kinetic energy are due to different 1s⁻¹ core-hole decay processes, which we outline in Fig. 3.

Among the peaks with constant kinetic energy, the most intense ones for all three ions are due to the local KLL Auger decay, shown in the leftmost panels of Fig. 2. The three ions are isoelectronic with configuration 1s²2s²2p⁶, which implies that their KLL Auger spectra resulting from the 1s⁻¹ → 2p⁻² + e⁻_Aug decay process should be (and are) relatively similar in terms of multiplet pattern, with the two peaks in this energy range representing the ¹S and ¹D states of the 2p⁻² configuration. Similar decays ending up in 2s⁻² and 2s⁻¹ 2p⁻¹ states occur at lower kinetic energies, although they are not shown in Fig. 2. In addition to the main KLL Auger decay, the intermediate M 1s⁻¹ state can also relax by filling the core hole by the electrons from the 2s or 2p orbitals and emitting an electron from one of the neighbouring molecules. This ICD process can be described as M 1s⁻¹ → M 2s⁻¹X⁻¹ + e⁻_ICD (ICD_2s) or M 1s⁻¹ → M 2p⁻¹X⁻¹ + e⁻_ICD (ICD_2p), with M and X designating the metal ion and a neighbouring molecule, respectively. These ICD electrons form the non-dispersing, broad, multipeak structures observed in the right-hand panels of Fig. 2 for the Mg²⁺ and Al³⁺ ions, but not for the Na⁺ ions. The energetics of the former features are further discussed in the Binding energies and ICD energies Section of the ESI† and compiled in Table S1 (ESI†).

ICD takes place between the ion and water

We will start by examining the species involved in the ICD process. The structure of the observed ICD signal should reflect the states corresponding to charge delocalization between the central ion and a neighbouring molecule (either a water molecule, as shown in Fig. 3, or a chloride anion). We assessed the contribution of the chloride anion in the ICD spectra of Al³⁺ and Mg²⁺ with the help of MD simulations. In the case of Al³⁺, our simulations showed only a very limited number of [Al(H₂O)₅Cl]²⁺ complexes at 2 M concentration: Only 1.5% of the observed structures featured a chloride anion in the first coordination shell. The majority of these structures also contained a second chloride anion, resulting in an average coordination number for chloride in the first hydration shell of only 0.04. For Mg²⁺, the MD simulations revealed that the fraction of contact ion pairs in solution is practically zero, in agreement with previous simulations as well as available neutron scattering and X-ray diffraction data.^40,56–58 The sodium cation forms a limited number of contact ion pairs, however, no ICD signal was observed, preventing the study of states in which the generated charge is delocalized between the central ion and the chloride anion. In the following discussion, in accord with the above, we assume that the ICD signal observed for Mg²⁺ and Al³⁺ is dominated by water molecules in the first solvation shell.

ICD intensity depends on the metal–ligand distance

We will compare the Na⁺, Mg²⁺, and Al³⁺ ICD spectra and focus on peak intensities next. Fig. 2 shows that the ratio of ICD intensities to the Auger signal is very different for each of the cations studied. Al³⁺ exhibits the highest relative ICD intensity, for Mg²⁺ the ICD intensity is significantly smaller, and for Na⁺ the ICD peaks effectively disappear. To quantify the differences, we compare the integrated intensities of the KLL Auger peaks (final state 2p⁻²) to the ICD_2p peaks (final state 2p⁻¹w⁻¹), i.e., we determine the intensity ratio of ICD to local Auger decay. We find that the intensity ratio I(ICD_2p)/I(KLL) is ∼0% for Na⁺, ∼2.6% ± 0.5% for Mg²⁺, and ∼5.1% ± 1% for Al³⁺. The theory of ICD predicts that the intensity should be asymptotically proportional to 1/R⁶, with R being the distance between the centres. The local solvation patterns around the three metal ions can be conveniently described by MD simulations via radial distribution functions (see Fig. S4 in the ESI†). The peak positions for the cation–water oxygen distribution are summarized in Table 1 together with respective coordination numbers. The mean of the distance between a cation and the instantaneously closest water oxygen is also provided in Table 1 as r_min. These results highlight the tighter coordination of the solvent to the metal ion in going from Na⁺ to Mg²⁺ to Al³⁺. In other words, the increasing charge of the cation correlates with a decrease of the intermolecular distance between the cation and neighboring water molecules. Accordingly, based on the MD simulations, the ICD signatures should be 1.3 times more intense for Al³⁺ than for the Mg²⁺ and 3.2 times more intense than for Na⁺. These ratios are reasonably consistent with the experimentally observed values. We note that the binding of water to a sodium cation is rather loose. Hence, although a large number of water molecules surround the Na⁺ ions, they barely contribute to the solute secondary ionization, i.e. ICD, signals.

Table 1 Structural parameters of the hydration shell for sodium, magnesium, and aluminium cations, according to MD simulations (this work). +O R_max refers to the position of the first maximum in the cation–oxygen radial distribution function, +O R_min is the position of the first minimum in that function, +On_{+ O} is the coordination number of water molecules, +O r_min is the mean distance between the cation and the closest water molecule and ±R_max is the position of the maximum in the radial distribution function between the respective cation and a chloride anion. All distances are in Å. ±R_max for Mg²⁺ is not given because this species does not form ion pairs

	Na^+	Mg^2+	Al^3+
+O Rmax	2.28	2.04	1.86
+O Rmin	3.12	2.78	2.16
+O n+ O	5.58	6.00	6.00
+O rmin	2.21	1.90	1.82
±Rmax	2.86	—	2.30

---

Having discussed how the change in ion–water distance from Na⁺ to Al³⁺ affects the total ICD intensity, we will now proceed to discuss a second observation related to the intensity, namely that the ICD_2p feature has about four times higher intensity than ICD_2s for both Mg²⁺ and Al³⁺. The two processes can be schematically written as 1s⁻¹→2s⁻¹w⁻¹+ e⁻_ICD (ICD_2s) and 1s⁻¹→2p⁻¹w⁻¹+ e⁻_ICD (ICD_2p). On a qualitative level, there are two reasons for this intensity difference. First, the 2s and 2p orbitals have different numbers of electrons, six in 2p and two in 2s. The number of ICD channels involving 2p is, therefore, higher than that of 2s. Second, as we mentioned before, the ICD probability scales as 1/R⁶ and the 2s and 2p orbitals have different spatial extents which is then reflected in relative ICD probability. Moreover, different directionality of the 2s and 2p orbitals also plays a role, for example, directional 2p orbitals might be oriented more favourably for ICD than the spherically symmetric 2s orbital.

The ICD process takes tens of femtoseconds

The observed I(ICD_2p)/I(KLL) intensity ratios depend on the relative branching ratios, and thus the relative timescales of the different decay channels. The 1s⁻¹ state lifetime t_1s is ∼ 2.3 fs for Na⁺, ∼2.0 fs for Mg²⁺, and ∼1.7 fs for Al³⁺,⁵⁹ and we found the experimentally observed I(ICD_2p)/I(KLL) intensity ratios as ∼0% for Na⁺, ∼2.6% for Mg²⁺, and ∼5.1% for Al³⁺. This allows us to use the core-hole clock, in which the core-hole lifetime is used as an internal timescale for the secondary decay processes.^60–63 We can then estimate the timescale for ICD as t_1s divided by the I(ICD_2p)/I(KLL) ratio, which results in ∼76 fs for Mg²⁺ and ∼34 fs for Al³⁺. The non-observation of ICD for Na⁺ supports the aforementioned weaker solute–solvent interaction and suggests that the ICD process occurs much more slowly for Na⁺ compared to Mg²⁺ and Al³⁺. Similar ICD processes following 2s ionization have been reported for Na⁺, Mg²⁺, and Al³⁺.²⁴ By line-shape analysis of the photoemission peaks, timescales of the IC-decays following 2s ionization were found as 3.1, 1.5, and 0.98 fs for Na⁺, Mg²⁺, and Al³⁺, respectively. These results follow the same trend as the lifetimes in the present case of ICD following 1s ionization: The decay channel becomes more efficient with a decrease in the ion–water distance, and with increasing nuclear charge Z. On an absolute level, the timescales of ICD following 2s ionization, however, are all substantially shorter. We tentatively interpret it as a consequence of a Coster–Kronig-like process – one of the final state holes, 2p⁻¹, is in the same shell as the primary 2s⁻¹ hole which makes the decay much more rapid. Radial matrix elements tabulated for the normal Auger decay of 1s and 2s holes in the metals in question do show an order of magnitude higher rate for the 2s decay.⁶⁴ Another contributing factor may be that the overlap between the water orbitals is larger with the 2s⁻¹ hole than the 1s⁻¹ hole.

ICD peaks reveal electronic structure of neighboring molecules

For both Mg²⁺ and Al³⁺ ions, we can see that the ICD_2p and ICD_2s channels exhibit a substructure (see Fig. 4 containing enlarged views of the Al³⁺ and Mg²⁺ ICD peaks). The final states of the two channels are 2s⁻¹w⁻¹ and 2p⁻¹w⁻¹, respectively. The w⁻¹ hole can be produced in any of the water valence orbitals, i.e., 1b₁, 3a₁, 1b₂, or 2a₁. Their shape, for the [Al(H₂O)₆]³⁺ complex, is exemplified in Fig. 5. The figure shows how the orientation of the orbital relative to the metal ion differs between the states. The observed substructure of the ICD peaks then corresponds to localization of the final-state hole in these different orbitals. The relative positions of the peaks used for fitting the ICD structure were found practically identical to the water valence-band peaks observed in photoemission.

Fig. 5 Selected molecular orbitals for the [Al(H₂O)₆]³⁺ complex. The Löwdin reduced orbital population per molecular orbital was performed at the BH&HLYP 6-31 + g* level in the polarizable continuum, respective molecular orbitals are depicted with an isovalue of 0.05 e.

The observed relative intensities of the peaks associated with the different water orbitals, however, are different for ICD_2p, ICD_2s, and the valence photoemission spectrum of the respective salt solution. (A valence-band spectrum for NaCl solution at the photon energy of our measurements is shown in the ESI,† Fig. S1.) For both Mg²⁺ and Al³⁺, the strongest contribution in ICD_2p is from the water 3a₁, followed by 1b₁ and 1b₂, with the 2a₁ orbital contribution being relatively low. In ICD_2s, however, the 1b₂ orbital peak has a slightly higher intensity than the 3a₁ orbital for both ions. This shows that the ICD process selects electrons from specific molecular orbitals in a way that differs from photoemission and even depends on the orbital involved in the metal ion.

A straightforward theoretical assignment of the ICD peaks would require calculations of decay rates (e.g., by means of the Fano theory^65,66). This is unfortunately intractable for the condensed phase. However, a hint about the observed shape of the ICD spectra can be provided by a simple orbital analysis; in our case, we selected the Löwdin population analysis. The preference for ICD electrons from specific molecular orbitals can conveniently be demonstrated for various molecular orbitals of water. The strongest ICD signal among the water molecular orbitals should arise from the 3a₁-type orbitals because they are aligned along the connecting line between water oxygen and the central metal ion (typical contribution amounts to 10–11%). The percentage describes a contribution of the water molecular orbitals on a given atomic orbital of a metal atom.⁶⁷ On the contrary, the 1b₁-type and 1b₂-type orbitals are oriented perpendicular to the connecting line between water oxygen and the metal ion, which is reflected in a much smaller orbital overlap (up to 2% at maximum). Reviewing the data and fits shown in Fig. 4, the orbital-overlap analysis is consistent with the observed relative intensities of the ICD_2p features, but the different relative intensities of the ICD_2s peaks show that the nature of the metal orbital also plays an important role. This calls for further research, especially as the present results indicate that ICD spectroscopy can represent a sort of ‘orbital tomography’ – in principle, it is possible to reconstruct dominant features of orbital shapes by a careful analysis of the ICD spectra following ionization of core electrons from different atoms.

It is interesting to compare this result to earlier work on core-level de-excitation processes involving several centres. A related decay mode of core-excited states is Electron Transfer Mediated Decay (ETMD), in which the energy transfer from centre to ligand, which is characteristic for ICD, is accompanied by an electron transfer to the centre ionized initially. This process has been studied, e.g., in the de-excitation of Li²⁺ core holes in aqueous Li⁺ electrolyte solutions by some of the authors.^11,12 Since solvated Li⁺ is devoid of any valence electrons, a 1s core hole (with a binding energy of 60.4 eV, ref. 12) can only decay when an electron is transferred from the solvation shell. Ionization by the excess energy released in this transfer has been experimentally observed, and similar to our current results showed a propensity for creating 3a₁ vacancies in the Li solvation shell, as found by accompanying calculations.¹¹ Since in that work a theoretical framework quite different from this article was used, a one-to-one comparison is not possible, but we consider it very plausible that in both cases the shape of the 3a₁ orbital is particularly suited to create overlap with the core hole which lends efficiency to this particular decay channel. As ETMD involves electron transfer, it is clear that it can only proceed if orbital overlap is given. In ICD, orbital overlap plays a minor role: In principle, the decay is possible for two completely separated entities, rather overlap may accelerate a decay channel which is open in any case.^6,68 Our work shows that this phenomenon can have an important impact on the shape of ICD spectra. This is reminiscent of results on the decay spectra of core holes in molecules featuring strongly electro-negative, e.g., fluorine, ligands. Here a strong, ICD-like involvement of the ligands in the core-hole decay was seen⁶⁹ and the orbital contributions of the ligands at the core-ionized centre were even found to determine the shape of the decay spectrum in some cases.⁷⁰

Water ligands differ from bulk water

Let us now discuss to what extent the electronic structure of the water molecules in the first solvation shell is affected by a neighbouring ion. It has been observed previously^23,71,72 that the presence of an electrolyte has only a minor effect on the photoemission spectrum of water, even for a very high concentration of the uni–univalent electrolyte. In the valence photoemission measurements, the signal is dominated by bulk water, and the electronic structure of water molecules in the first solvation shell is very difficult to identify. In contrast, the ICD process involves practically only water molecules in the first solvation shell of the core-ionized cation. This opens up the possibility to selectively probe the electronic structure of these bound water molecules, to see how they differ from bulk water. Here, we compare the valence photoemission spectrum of bulk water and photoemission spectra of the water molecules in the first solvation shell, as simulated by the IEDC technique.

The calculated valence photoemission spectra of pure water and water molecules directly coordinated to the cations are shown in Fig. 6, the respective peak positions are collected in Table S5 of the ESI.† It can be seen that the calculations reproduce the experimental binding energies for pure water within tenths of eV. We can observe that the cations significantly influence the binding energies of the closest water molecules. For the magnesium cation, the binding energies are shifted towards higher values (by about 0.3–0.5 eV), and the effect of the aluminium cation is even stronger, the values are shifted by about 0.6–1.3 eV. In both cases, the six water molecules forming the solvation shell are more strongly bound to the cation than to another water molecule. In the case of aluminium, the 3a₁ peak forms a double-peak structure, the interpretation of which is unclear. Such an increase of binding energies can be expected by taking into account the electric field of a (positively charged) metal ion. This claim, however, is too simplistic, as the sodium cation exhibits the opposite effect. For Na⁺, the energies are shifted towards lower energies by 0.5–0.7 eV because Na⁺ ‘breaks’ the structure of liquid water, e.g., it forms weaker bonds with water.⁷³ Note that the shift in binding energies of water in the solvation shell, ΔE_b,vi, should be taken into account in the theoretical interpretation of ICD spectra, e.g., experimental values for water valence-band peaks should not be used.

Fig. 6 Simulated photoemission spectrum for pure water, and for water molecules directly coordinated to a metal cation in NaCl, MgCl₂, and AlCl₃ solutions. Spectra were calculated by the QM/MMPol IEDC approach at the LC-ω PBE/cc-pVTZ level.

Energetics of the ICD process

We now further discuss the energetics of the ICD process, M 1s⁻¹ → M 2s/2p⁻¹w⁻¹ + e⁻_ICD. As can be seen in Fig. 3, the kinetic energy of the ICD electron upon 1s ionization (E_k) can be expressed as the difference between the intermediate M 1s⁻¹ (E_1s) and M 2s⁻¹/2p⁻¹w⁻¹ two-hole final state energies (E_2h)

E_k = E_1s − E_2h. (1)

We can obtain the energy of the final two-hole states from the experimental data as

E_2h = E_1s − E_k. (2)

The resulting two-hole energies for both the Auger and ICD final states of Mg²⁺ and Al³⁺ are shown in Fig. 7 and collected in Tables S2 and S3 in the ESI.† For Mg²⁺ (Al³⁺), the experimentally obtained two-hole energies are seen to decrease from 134 (187) eV for the localized 2p⁻² final states of the Auger decay to 70 (96) eV for 2p⁻¹ 1b₁⁻¹, the lowest delocalized two-hole state reached in the ICD process. This decrease in energy is an important driving force for the ICD process, and for the localized 2p⁻² final states this has been predicted to lead to secondary delocalization via the ETMD processes.⁷⁴

Fig. 7 The energies of the two-hole ICD final states for Al³⁺, Mg²⁺, and Na⁺. The experimental values, measured with 1569.8 eV and 1315.25 eV photon energies for Al³⁺ and Mg²⁺, respectively, were determined from the curve fitting of the measured spectra. MOM/LC-ω PBE/cc-pVTZ (with cc-cPVTZ basis for the cation) energies were obtained for a minimal model containing one cation and one water molecule; the dimer was solvated by a 20 Å sphere of molecules treated at the MMPol level. The values of experimental and calculated energies are collected in Tables S2, S3, and S6 in the ESI.† Experimental errors are approximately equal to or smaller than the symbol size. No experimental values for delocalized final states are given for Na⁺, as ICD leading to the respective two-hole states could not be observed.

To quantitatively understand this energy lowering, we can conceptually decompose the energy of the final two-hole states into three contributions

E_2h = E_b,v1 + E_b,v2 + E_Cp, (3)

where E_b,vi is the binding energy of the i^th electron participating in this process and E_Cp is the Coulomb penalty, which can be understood as the electrostatic potential energy of the positively charged metal and water ions in a close proximity. The values of E_b,vi and E_Cp are not directly experimentally accessible; the experiment provides bulk water valence-band peaks only. Ignoring this difference, the Coulomb penalty is generally found in the range of 3–6 eV, with a tendency towards somewhat lower values if states of s-character are involved. This is in contrast with the Coulomb penalty found for the ETMD process following the Li⁺ ionization where the Coulomb penalty was found to be very small.¹¹

Theoretically, we model the energetics with the MOM approach in the minimal model containing one cation and one water molecule for Na⁺, Mg²⁺, and Al³⁺ in a polarizable embedding. In this way, we obtain the binding energies of the i^th electron E_b,vi, the two-hole state energies, and the Coulomb penalty. The calculated energies are presented in Fig. 7 and collected in Tables S2 and S3 in the ESI.†

As seen in Fig. 7, the calculations reproduce the experimental values for the final two-hole energies for Mg²⁺ and Al³⁺ very well. It is important to note here that the calculated E_2h values correspond to the low-energy onset of the fitted ICD peaks (the MOM approach provides only the energies of the lowest triplet two-hole states). For Mg²⁺ and Al³⁺, the Coulomb penalty is found in the range 4–6 eV, which is reasonably consistent with the experimental estimate.

For Na⁺, there are no experimental values since we did not observe any ICD-related spectral features, but based on the good agreement between the calculation including solvation and experiment for Mg²⁺ and Al³⁺, we can use the calculated values for Na⁺. One interesting observation is that the energy difference between the localized 2p⁻² final state to the delocalized 2p⁻¹w⁻¹ states is much smaller for Na⁺ than for Mg²⁺ and Al³⁺. If we regard this lowering of the two-hole state energy as a driving force for the ICD process, this may be another reason why ICD is less efficient for Na⁺ than for Mg²⁺ and Al³⁺.

As can be inferred from Table S6 in the ESI,† the theoretically predicted Coulomb penalty E_Cp is higher for Al³⁺ and Mg²⁺ than for Na⁺, which is due to the formation of highly charged pairs of Al⁴⁺ H₂O⁺ and Mg³⁺H₂O⁺ compared to Na²⁺H₂O⁺ as the final ICD states. The energy difference between the ions is however rather minor, most probably due to a large amount of screening by the surrounding water.

Screening effects of water

The two-hole energies for the minimal system in the gas phase are collected in Table S8 in the ESI,† respective one-electron binding energies are provided in Table S9 (ESI†). These energies exhibit the same trend for energy lowering upon delocalization as the minimal system with solvation, but the absolute two-hole energies are much higher, ≈20 eV for Na⁺, ≈40 eV for Mg²⁺, and ≈60 eV for Al³⁺. These large differences observed in the gas phase further highlight the huge screening effect of water, which minimizes the Coulomb penalty and the binding-energy shifts. The Coulomb penalty in the gas phase is a factor 1.5–2 larger than in solution; the binding-energy shifts are around 7 eV for Na⁺, 15 eV for Mg²⁺, and more than 20 eV for Al³⁺ in the gas phase compared to less than 1.3 eV in the aqueous phase. The large reduction of the two-hole energies illustrates the importance of including the solvent effects when considering the energetics of ICD.

4 Conclusions

X-ray photoemission spectroscopy, molecular dynamics, and ab initio calculations have been combined to investigate ICD channels during the decay of 1s core holes in aqueous-phase Al³⁺, Mg²⁺, and Na⁺ in chloride solutions. We observed two individual ICD features originating from the relaxation involving either the 2s or 2p levels of the cations and mostly water valence states. Importantly, each of these features displays a substructure, which has been shown to originate from the participation of different water valence orbitals. According to our calculations and the analysis of our data, transitions including the water 3a₁ level contribute the largest share of the overall ICD signal. This is caused by the orbital overlap between the water valence states and the cation. The orientation of the water molecules toward the cation enhances the orbital overlap of the water 3a₁ orbital, whereas the water 1b₁ and 1b₂ orbital overlaps are disfavoured. However, ICD channels that lead to 1b₂ orbital emission appear to be favoured in the specific cases of Mg²⁺ and Al³⁺ ICD_2s processes relative to ICD_2p, indicating that also the electronic structure of the core-hole-excited species plays a role. If chloride is present in the first solvation shell (forming a contact ion pair), we expect a significant overlap of one of the Cl⁻ 2p orbitals, too. However, according to our calculations and the experimental data, this is rarely the case. Contributions from Cl⁻ states were not discernable in the experiments.

The charge of the cation has been found to have a profound impact on the overall strength of the ICD transitions. For example, the higher the charge on the cation, the stronger the ICD signal. ICD is pronounced in the decay of 1s holes in Al³⁺ and Mg²⁺, while only Auger decay is observable in Na⁺. The enhancement of the ICD signal with the charge of the cation reflects the intermolecular distances between the cation and the water molecules in the first solvation shell, as ICD scales with 1/R⁶, with R being the distance between the interacting entities (see entries r_min in Tab. 1 for calculated, averaged values). Accordingly, the strongest ICD_2p/Auger branching ratio could be determined for Al³⁺ with about ∼5% relative ICD efficiency. This number is about a factor of 10 higher than ICD rates, e.g., in Ar clusters.⁷⁵ A ∼5% intensity ratio also allows us to estimate the timescale of the ICD_2p transition to about 34 fs, based on the core-hole lifetime in Al 1s of about 1.7 fs. Note that this is about 20 times slower than the ICD decay of much shallower Al 2s holes.¹⁴

We have shown that the substructure of ICD spectra hinges on the electronic structure of the valence levels of the partner species during the de-excitation. One may regard the ICD spectrum as a projection of the valence states, but with the relative intensities and peak widths modified by the orbital overlap and the geometry of the reaction partners. If used under this premise, ICD can serve as a quasi-tomographic tool to explore the combined electronic structure and local geometry of aqueous solutes and their first hydration shell simultaneously.

From a radiation-chemistry perspective it is interesting to follow the further relaxation of the ICD final state. We provide a speculative outlook here: A Coulomb explosion will immediately follow the ICD in the gas phase⁶ and in aggregates of inert species.⁷⁶ However, we suggest that the process is more complicated in the liquid phase. An ionized water molecule, H₂O⁺, is an extremely short-lived species; in less than 50 fs, H₂O⁺ gives up a proton to another water molecule, forming H₃O⁺ and a hydroxyl radical.⁷⁷ Besides, we can speculate that a charge transfer from another water molecule to the metal cation might follow;⁷⁴ this process would give rise to another ionized water molecule and another proton transfer would follow. The overall result would then be formation of two OH radicals and two H₃O⁺, which would later slowly diffuse.

Data Availability

Data relevant for this study are available at DOI: 10.5281/zenodo.6372662.

Author contributions

G. G. and O. B. conceived the project, G. G., S. M., F. T., I. W., O. B., and B. W. performed the experiments, G. G., I. U., G. Ö. and U. H. analyzed the data, E. M. and P. S. performed all the ab initio calculations, B. M. and F. L. adjusted and implemented MOM into the polarizable force field, G. G., E. M., P. S., O. B., and U. H. wrote the manuscript that was discussed by all authors.

Conflicts of interest

The authors declare no competing financial interest.

Acknowledgements

The authors would like to thank Chin Lee, Dana Bloß, Karen Mudryk, Wilson Quevedo, Tillmann Buttersack, and Clemens Richter for their help in data acquisition, and Robert Seidel for his help with the SOL³PES setup. We acknowledge DESY (Hamburg, Germany), a member of the Helmholtz Association HGF, for the provision of experimental facilities. Parts of this research were carried out at PETRA III and we would like to thank Moritz Hoesch and his team for assistance in using beamline P04. We also thank HZB for the allocation of synchrotron-radiation beamtime. F. T. and B. W. acknowledge support by the MaxWater initiative of the Max-Planck-Gesellschaft. B. W. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No. 883759. P. S. and E. M. thank the Czech Science Foundation for the support via project number 21-26601X (EXPRO project). E. M. acknowleges that the research leading to the results was supported by the European Structural and Investment Funds, OP RDE-funded project ‘CHEMFELLS III’ (CZ.02.2.69/0.0/0.0/19_074/0014006) funding her MCSA-IF proposal. O. B. acknowledges funding from the Swedish Research Council (VR) for the project VR 2017-04162. C. C. acknowledges the Swedish Research Council, project 2018-00740, and the Helmholtz Association through the Center for Free-Electron Laser Science at DESY. Open Access funding provided by the Max Planck Society.

Notes and references

1. P. Jungwirth and D. J. Tobias, Chem. Rev., 2006, 106, 1259–1281.
2. N. F. A. van der Vegt, K. Haldrup, S. Roke, J. Zheng, M. Lund and H. J. Bakker, Chem. Rev., 2016, 116, 7626–7641.
3. O. Björneholm, M. H. Hansen, A. Hodgson, L.-M. Liu, D. T. Limmer, A. Michaelides, P. Pedevilla, J. Rossmeisl, H. Shen, G. Tocci, E. Tyrode, M.-M. Walz, J. Werner and H. Bluhm, Chem. Rev., 2016, 116, 7698–7726.
4. L. S. Cederbaum, J. Zobeley and F. Tarantelli, Phys. Rev. Lett., 1997, 79, 4778.
5. S. Marburger, O. Kugeler, U. Hergenhahn and T. Möller, Phys. Rev. Lett., 2003, 90, 203401.
6. T. Jahnke, U. Hergenhahn, B. Winter, R. Dörner, U. Frühling, P. V. Demekhin, K. Gokhberg, L. S. Cederbaum, A. Ehresmann, A. Knie and A. Dreuw, Chem. Rev., 2020, 120, 11295–11369.
7. K. Gokhberg, P. Kolorenč, A. I. Kuleff and L. S. Cederbaum, Nature, 2014, 505, 661–663.
8. F. Trinter, M. S. Schöffler, H.-K. Kim, F. P. Sturm, K. Cole, N. Neumann, A. Vredenborg, J. Williams, I. Bocharova, R. Guillemin, M. Simon, A. Belkacem, A. L. Landers, T. Weber, H. Schmidt-Böcking, R. Dörner and T. Jahnke, Nature, 2014, 505, 664–666.
9. P. Slavíček, N. V. Kryzhevoi, E. F. Aziz and B. Winter, J. Phys. Chem. Lett., 2016, 7, 234–243.
10. S. Thürmer, M. Ončák, N. Ottosson, R. Seidel, U. Hergenhahn, S. E. Bradforth, P. Slavíček and B. Winter, Nat. Chem., 2013, 5, 590–596.
11. I. Unger, R. Seidel, S. Thürmer, M. N. Pohl, E. F. Aziz, L. S. Cederbaum, E. Muchová, P. Slavíček, B. Winter and N. V. Kryzhevoi, Nat. Chem., 2017, 9, 708–714.
12. M. N. Pohl, C. Richter, E. Lugovoy, R. Seidel, P. Slavíček, E. F. Aziz, B. Abel, B. Winter and U. Hergenhahn, J. Phys. Chem. B, 2017, 121, 7709–7714.
13. W. Pokapanich, H. Bergersen, I. L. Bradeanu, R. R. T. Marinho, A. Lindblad, S. Legendre, A. Rosso, S. Svensson, O. Björneholm, M. Tchaplyguine, G. Öhrwall, N. V. Kryzhevoi and L. S. Cederbaum, J. Am. Chem. Soc., 2009, 131, 7264–7271.
14. N. Ottosson, G. Öhrwall and O. Björneholm, Chem. Phys. Lett., 2012, 543, 1–11.
15. S. Barth, S. Marburger, S. Joshi, V. Ulrich, O. Kugeler and U. Hergenhahn, Phys. Chem. Chem. Phys., 2006, 8, 3218–3222.
16. E. Fasshauer, M. Förstel, S. Pallmann, M. Pernpointner and U. Hergenhahn, New J. Phys., 2014, 16, 103026.
17. A. A. Chialvo and L. Vlcek, Fluid Phase Equilib., 2016, 407, 84–104.
18. B. Winter and M. Faubel, Chem. Rev., 2006, 106, 1176–1211.
19. B. Winter, R. Weber, I. V. Hertel, M. Faubel, P. Jungwirth, E. C. Brown and S. E. Bradforth, J. Am. Chem. Soc., 2005, 127, 7203–7214.
20. D. Ghosh, A. Roy, R. Seidel, B. Winter, S. Bradforth and A. I. Krylov, J. Phys. Chem. B, 2012, 116, 7269–7280.
21. A. P. Gaiduk, M. Govoni, R. Seidel, J. H. Skone, B. Winter and G. Galli, J. Am. Chem. Soc., 2016, 138, 6912–6915.
22. T. A. Pham, M. Govoni, R. Seidel, S. E. Bradforth, E. Schwegler and G. Galli, Sci. Adv., 2017, 3, e1603210.
23. M. N. Pohl, E. Muchová, R. Seidel, H. Ali, Š. Sršeň, I. Wilkinson, B. Winter and P. Slavíček, Chem. Sci., 2019, 10, 848–865.
24. G. Öhrwall, N. Ottosson, W. Pokapanich, S. Legendre, S. Svensson and O. Björneholm, J. Phys. Chem. B, 2010, 114, 17057–17061.
25. S. Malerz, H. Haak, F. Trinter, A. B. Stephansen, C. Kolbeck, M. Pohl, U. Hergenhahn, G. Meijer and B. Winter, Rev. Sci. Instrum., 2022, 93, 015101.
26. J. Viefhaus, F. Scholz, S. Deinert, L. Glaser, M. Ilchen, J. Seltmann, P. Walter and F. Siewert, Nucl. Instrum. Methods Phys. Res., Sect. A, 2013, 710, 151–154.
27. P04 beamline parameters, https://photon-science.desy.de/facilities/petra_iii/beamlines/p04_xuv_beamline/beamline_parameters/index_eng.html, Accessed: 2021-09-23.
28. R. Seidel, M. N. Pohl, H. Ali, B. Winter and E. F. Aziz, Rev. Sci. Instrum., 2017, 88, 073107.
29. J. Duan and J. Gregory, Adv. Colloid Interface Sci., 2003, 100–102, 475–502.
30. Helmholtz-Zentrum Berlin für Materialien und Energie, Journal of Large-Scale Research Facilities, 2016, 2, A72.
31. E. Kukk, Spectrum Analysis by Curve Fitting (SPANCF) macro package for Igor Pro, 2012.
32. B. Winter, R. Weber, W. Widdra, M. Dittmar, M. Faubel and I. V. Hertel, J. Phys. Chem. A, 2004, 108, 2625–2632.
33. D. Loco, É. Polack, S. Caprasecca, L. Lagardère, F. Lipparini, J.-P. Piquemal and B. Mennucci, J. Chem. Theory Comput., 2016, 12, 3654–3661.
34. M. Bondanza, M. Nottoli, L. Cupellini, F. Lipparini and B. Mennucci, Phys. Chem. Chem. Phys., 2020, 22, 14433–14448.
35. J. W. Ponder, C. Wu, P. Ren, V. S. Pande, J. D. Chodera, M. J. Schnieders, I. Haque, D. L. Mobley, D. S. Lambrecht, R. A. DiStasio, Jr., M. Head-Gordon, G. N. I. Clark, M. E. Johnson and T. Head-Gordon, J. Phys. Chem. B, 2010, 114, 2549–2564.
36. P. Ren and J. W. Ponder, J. Phys. Chem. B, 2003, 107, 5933–5947.
37. F. Lipparini, L. Lagardère, C. Raynaud, B. Stamm, E. Cancès, B. Mennucci, M. Schnieders, P. Ren, Y. Maday and J.-P. Piquemal, J. Chem. Theory Comput., 2015, 11, 623–634.
38. F. Lipparini, J. Chem. Theory Comput., 2019, 15, 4312–4317.
39. L. Greengard and V. Rokhlin, J. Comput. Phys., 1987, 73, 325–348.
40. E. Duboué-Dijon, P. E. Mason, H. E. Fischer and P. Jungwirth, J. Phys. Chem. B, 2018, 122, 3296–3306.
41. E. Cauët, S. A. Bogatko, E. J. Bylaska and J. H. Weare, Inorg. Chem., 2012, 51, 10856–10869.
42. T. S. Hofer, B. R. Randolf and B. M. Rode, J. Phys. Chem. B, 2008, 112, 11726–11733.
43. I. S. Joung and T. E. Cheatham, J. Phys. Chem. B, 2008, 112, 9020–9041.
44. G. Lamoureux and B. Roux, J. Phys. Chem. B, 2006, 110, 3308–3322.
45. H. Yu, T. W. Whitfield, E. Harder, G. Lamoureux, I. Vorobyov, V. M. Anisimov, A. D. MacKerell, Jr. and B. Roux, J. Chem. Theory Comput., 2010, 6, 774–786.
46. Y. Marcus, Chem. Rev., 1988, 88, 1475–1498.
47. A. T. B. Gilbert, N. A. Besley and P. M. W. Gill, J. Phys. Chem. A, 2008, 112, 13164–13171.
48. M. Rubešová, E. Muchová and P. Slavíček, J. Chem. Theory Comput., 2017, 13, 4972–4983.
49. E. Muchová and P. Slavíček, J. Phys.: Condens. Matter, 2019, 31, 043001.
50. J. F. Stanton and J. Gauss, J. Chem. Phys., 1999, 111, 8785–8788.
51. S. Coriani and H. Koch, J. Chem. Phys., 2015, 143, 181103.
52. F. Neese, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1327.
53. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, Gaussian Development Version Revision J.06+, Gaussian Inc., Wallingford CT, 2020.
54. M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess and E. Lindahl, SoftwareX, 2015, 1-2, 19–25.
55. S. Pronk, S. Páll, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, B. Hess and E. Lindahl, Bioinformatics, 2013, 29, 845–854.
56. S. Friesen, G. Hefter and R. Buchner, J. Phys. Chem. B, 2019, 123, 891–900.
57. K. M. Callahan, N. N. Casillas-Ituarte, M. Roeselová, H. C. Allen and D. J. Tobias, J. Phys. Chem. A, 2010, 114, 5141–5148.
58. F. Bruni, S. Imberti, R. Mancinelli and M. A. Ricci, J. Chem. Phys., 2012, 136, 064520.
59. J. Campbell and T. Papp, At. Data Nucl. Data Tables, 2001, 77, 1–56.
60. O. Björneholm, A. Nilsson, A. Sandell, B. Hernnäs and N. Mårtensson, Phys. Rev. Lett., 1992, 68, 1892–1895.
61. W. Wurth and D. Menzel, Chem. Phys., 2000, 251, 141–149.
62. P. A. Brühwiler, O. Karis and N. Mårtensson, Rev. Mod. Phys., 2002, 74, 703–740.
63. A. Föhlisch, P. Feulner, F. Hennies, A. Fink, D. Menzel, D. Sanchez-Portal, P. M. Echenique and W. Wurth, Nature, 2005, 436, 373–376.
64. M. H. Chen, F. P. Larkins and B. Crasemann, At. Data Nucl. Data Tables, 1990, 45, 1–205.
65. U. Fano, Phys. Rev., 1961, 124, 1866–1878.
66. H. Feshbach, Rev. Mod. Phys., 1964, 36, 1076–1078.
67. P.-O. Löwdin, J. Chem. Phys., 1950, 18, 365–375.
68. V. Averbukh, I. B. Müller and L. S. Cederbaum, Phys. Rev. Lett., 2004, 93, 263002.
69. T. D. Thomas, C. Miron, K. Wiesner, P. Morin, T. X. Carroll and L. J. Saethre, Phys. Rev. Lett., 2002, 89, 223001.
70. F. O. Gottfried, L. S. Cederbaum and F. Tarantelli, J. Chem. Phys., 1996, 104, 9754–9767.
71. S. Thürmer, S. Malerz, F. Trinter, U. Hergenhahn, C. Lee, D. M. Neumark, G. Meijer, B. Winter and I. Wilkinson, Chem. Sci., 2021, 12, 10558–10582.
72. B. Credidio, M. Pugini, S. Malerz, F. Trinter, U. Hergenhahn, I. Wilkinson, S. Thürmer and B. Winter, Phys. Chem. Chem. Phys., 2022, 24, 1310–1325.
73. I. Waluyo, C. Huang, D. Nordlund, U. Bergmann, T. M. Weiss, L. G. M. Pettersson and A. Nilsson, J. Chem. Phys., 2011, 134, 064513.
74. V. Stumpf, K. Gokhberg and L. S. Cederbaum, Nat. Chem., 2016, 8, 237–241.
75. A. Hans, C. Küstner-Wetekam, P. Schmidt, C. Ozga, X. Holzapfel, H. Otto, C. Zindel, C. Richter, L. S. Cederbaum, A. Ehresmann, U. Hergenhahn, N. V. Kryzhevoi and A. Knie, Phys. Rev. Res., 2020, 2, 012022(R).
76. F. Wiegandt, F. Trinter, K. Henrichs, D. Metz, M. Pitzer, M. Waitz, E. Jabbour al Maalouf, C. Janke, J. Rist, N. Wechselberger, T. Miteva, S. Kazandjian, M. Schöffler, N. Sisourat, T. Jahnke and R. Dörner, Phys. Rev. A, 2019, 100, 022707.
77. Z.-H. Loh, G. Doumy, C. Arnold, L. Kjellsson, S. H. Southworth, A. Al Haddad, Y. Kumagai, M.-F. Tu, P. J. Ho, A. M. March, R. D. Schaller, M. S. Bin Mohd Yusof, T. Debnath, M. Simon, R. Welsch, L. Inhester, K. Khalili, K. Nanda, A. I. Krylov, S. Moeller, G. Coslovich, J. Koralek, M. P. Minitti, W. F. Schlotter, J.-E. Rubensson, R. Santra and L. Young, Science, 2020, 367, 179–182.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d2cp00227b

‡ These authors contributed equally to this work.

---