Metallic water: Transient state under ultrafast electronic excitation

Being omnipresent in biological materials, water is the prototypical sample for studies of living organisms. The response of water to various types of radiation is used as the standard in radiation safety research and radiation therapy.^1–5 The disintegration of water under irradiation into various fragments (different species)—water radiolysis—was intensively studied over decades.^1–5 It plays a crucial role in biochemistry, defining survival rates of bio-molecules embedded in water solutions—in particular, living cells. Response of water molecules to high-energy photons such as x-rays and gamma rays, electrons, and fast ions (in material science called swift heavy ions, SHI, also known in bioscience as HZE, which stands for high-Z and energy) was a subject of numerous research.^2,6–8

Typically, it is assumed that each incoming high-energy particle produces a set of well-defined fragments embedded in an unaffected water environment. A small number of water molecules disintegrate into typical products: free radicals H, OH, O, and free electrons, and further chemical kinetics may lead to damage to other water molecules and create secondary products via interaction of the primary ones (such as H₂O₂, H₂, O₂, etc.).⁹ Since different radicals damage biological materials at different rates, it is important to trace the evolution of concentrations of them: e.g., the most damaging for human cells are considered reactive oxygen species (ROS), such as hydroxyl (OH⁻) and superoxide (⁠$O2−$⁠).^10–12 

The processes of formation and evolution of various kinds of radicals are usually described with a set of chemical rate equations, or corresponding kinetic Monte Carlo procedures.^9,13 Both methods are based on the assumption that the chemical kinetics of created species takes place in normal water media. These conditions presume low fluxes (dose rates) of incoming radiation, such that the bulk of the water is either unaffected or has sufficient time to relax between impacts of successive particles.

These conditions are not satisfied in many irradiation scenarios. Modern femtosecond lasers, especially x-ray free-electron lasers (XFEL) produce extremely short and intense pulses, which ionize the target to high degrees, such that the densities of excited electrons approach the densities of atoms.^14,15 Comparable levels of excitation may transiently be produced in the nanometric proximity of the swift heavy ion trajectories.^16,17

The excitation of a massive amount of electrons leads to the formation of unusual phases in materials.^18,19 When electrons are excited to high temperatures, whereas ions of the target transiently stay relatively cold, the matter departs from the equilibrium phase diagram and forms new states, unachievable otherwise. The physical and chemical characteristics of such states may be drastically different from the equilibrium phases and require dedicated studies. Recently, it has been reported that water irradiated with a femtosecond XFEL pulse becomes highly electronically conducting.²⁰ It was attributed to the transition of water into a liquid metal state, however, the exact mechanisms were unclear since results of density-functional theory (DFT) simulations did not agree with the experiment.²⁰ 

Here, we demonstrate that under ultrafast excitation of the electronic system, bulk water transiently forms a metallic state via the collapse of the bandgap, making it electronically conducting. This transition is induced by an interplay of the nonthermal modification of the interatomic potential and enhanced electron–ion coupling. It is accompanied by the softening of the interatomic potential, eventually leading to the disintegration of the water molecules into various species, depending on the deposited dose (electronic temperature reached): at low doses, the common reactive species such as OH are produced, whereas, with the increase of the dose, they fragment into atomic H and O species.

Upon ultrafast irradiation, such as with XFEL or SHI, primarily the energy is deposited into the electronic system of the target, whereas its ionic system stays cold.^14,21,22 The material response started with the excitation of electrons, proceeds via secondary cascades, exciting new electrons from the valence band or core shells to the conduction band of the material. At the same time, created core–shell holes decay via the Auger channel, also exciting secondary electrons.²³ Secondary electron cascades typically finish within a few (or a few tens) of femtoseconds, when all electrons lose their energy below a certain threshold and cannot ionize new electrons.²⁴ The majority of electrons are in local thermodynamic equilibrium even during the cascading stage, with only a minority of electrons performing the cascades.¹⁴ 

With femtosecond pulse duration and fluences typically produced at XFELs, the electronic system of the target can reach temperatures of a few electron-Volts.^14,15,25 At such densities and energies of the electronic ensemble, it quickly thermalizes and forms the two-temperature state: with an electronic temperature high above the atomic one.¹⁴ 

The effect of the high-temperature electronic system is twofold: it affects the interatomic potential—collective potential energy surface in the material (nonthermal effect), and it starts to exchange energy with ions via electron–ion coupling. The first one may lead to molecular destabilization and result in nonthermal phase transitions, for example, nonthermal melting.^26,27 These effects can be studied with ab initio methods, such as density-functional theory molecular dynamics (DFT-MD) within the Born–Oppenheimer (BO) approximation.²⁸ It assumes that electrons are much lighter and faster than the ions, and thus, electronic wave functions instantly adjust to atomic displacements. The BO approximation, thus, does not allow for non-adiabatic coupling between electronic and atomic movements (electron–ion scattering, including normal atomic oscillations), which triggers electronic transitions between their energy levels, thereby exchanging kinetic energy between the electrons and ions.²⁹ The latter process is responsible for the equilibration of the electronic and atomic temperatures in the two-temperature state.¹⁴ 

Here, we employ two different methods to study the above-mentioned effects. The first ab initio method is the DFT-MD, which relies on the BO approximation. It allows us to study precisely the effects of electronic excitation to high temperatures on the interatomic forces and the response of the ionic system to it. The second one is based on tight-binding molecular dynamics (TBMD). It is a more approximate method than DFT, but it allows us to independently validate our conclusions, treat larger systems, and, most importantly, study the nonadiabatic coupling effects. Below we briefly describe the two methods.

Both approaches enable calculations of the evolution of the electronic energy levels (molecular orbitals, or band structure), the electronic wave functions, and the potential energy surface—forces acting on atoms. Excitation of electrons directly affects all these, thereby allowing us to take into account nonthermal effects.

The DFT-MD tool is employed via the Quantum Espresso simulation package.³⁰ Within this approach, Kohn–Sham one-particle equations are solved at each MD step to obtain electronic energy levels and to compute interatomic forces.³¹ Perdew–Burk–Ernzerhof’s (PBE³²) approximation for the exchange–correlation energy functional—the key parameter in DFT—was chosen to simulate a 96-atom cell. Although it describes water at ambient conditions worse than more advanced functionals, such as SCAN,³³ and it underestimates band gaps at ambient conditions, it performs much better at elevated electronic temperatures and provides qualitatively the same results as more accurate but slower-performing meta-GGA functionals.³⁴ As was shown in Refs. 34, GW calculations modified with the Keldysh diagram technique produce very close results to local-density approximation (LDA) functional at electronic temperatures above T_e ∼ 1 eV, indicating that even the LDA (and, by extension, more advanced functionals) should be applicable at elevated electronic temperatures.

To remove core electrons from the simulation, we use fhi98PP norm-conserving pseudopotentials³⁵ from the ABINIT website library.³⁶ The size of the plane wave basis set is controlled by energy cutoff E_cutoff = 170 Ry (2313 eV). The molecular dynamics use the Verlet algorithm with the time step 0.5 fs.

The second tool used here is the earlier developed hybrid code XTANT-3.^14,37 In XTANT-3, the electronic band structure is calculated using the transferrable tight binding (TB) method.¹⁴ We applied mio-1-1 DFTB parameterization, which employs the sp³ basis set for the description of water, specifically developed to describe biological systems.³⁸ 

To trace the nonadiabatic energy exchange between the electrons of the valence and conduction band with atoms, XTANT-3 uses Boltzmann collision integrals (BCI) formalism.³⁷ The electron–ion collision integral depends on the transient matrix elements of the electron–ion interaction, and the electron populations (distribution function).

The atoms are traced with the molecular dynamics (MD) simulation. The interatomic potential, evolving together with the electronic system, is provided by the transient TB Hamiltonian. The energy transferred from electrons calculated with the BCI is fed to atoms via the velocity scaling algorithm at each time step of the simulation.³⁷ We apply the fourth-order Martyna–Tuckerman algorithm to propagate atomic coordinates with the timestep of 0.05 fs.³⁹ Typically, 192 atoms are used in the simulation box with periodic boundary conditions.

For both, DFT-MD and TBMD simulations, before productive runs, the water state was created with the following procedures: water molecules were placed randomly in the simulation box that was adjusted in size to produce normal water density for the chosen number of molecules. The starting rotational orientation of molecules was also chosen randomly. Then, the system was relaxed with the steepest descent algorithm to reach the minimum of the potential energy (the produced configurations can be found in the supplementary material). After that, the simulations were initialized with atoms in those positions and random velocities assigned according to the Maxwellian distribution (at room temperatures). Before the increase of the electronic temperature, mimicking the arrival of an x-ray pulse (with the chosen duration of 10 fs FWHM), the system was allowed to equilibrate for a few hundred femtoseconds.³⁷ 

XTANT-3 has previously been validated by comparisons with experimental data for various materials, see, e.g., Refs. 14, 37, and 40, and thus, we assume that it can be applied here to water. We validate its performance at elevated electronic temperatures against the DFT-MD calculations in Sec. III.

DFT-MD simulations predict changes in water band structure starting from the electronic temperatures of T_e ∼ 2.5 eV (the deposited dose of 1.6 eV/atom, 5.1% of valence electrons excited to the conduction band), while the atomic system starts at room temperature. At this dose [see Fig. 1(a)], a part of the conduction band energy levels starts to descend at ∼380 fs forming impurity-like levels in the bandgap.

A further increase of the electronic temperature results in faster creation of such levels also making the energy spectrum inside of the former bandgap denser. At T_e ∼ 2.75 eV (2 eV/atom, 6.2% of valence electrons excited to the conduction band) first levels of the conduction band reach the top of the valence band as fast as ∼50 fs [Fig. 1(b)].

At last, at T_e ∼ 3 eV (2.5 eV/atom, 7.2% of valence electrons excited to the conduction band) the bandgap collapses completely within ∼25 fs merging the valence and the conduction bands into a quasi-continuous energy spectrum [Fig. 1(c)]. This indicates a formation of the electronically conducting state, a liquid metal.

Figure 2 demonstrates the radial pair distribution functions of water at T_e ∼ 3 eV. A significant decrease of the O–H correlation peak at ∼1 Å indicates that the vast majority of O–H bonds are broken already at ∼25 fs after electronic temperature elevation, which coincides with the timescale of the bandgap collapse. At the same time, O–O peak is almost unchanged, demonstrating that the relative positions of water molecules are the same as in the original liquid state. This shows that the bandgap stability depends on the stability of molecular bonds rather than on the mutual positioning of water molecules.

The rise of the H–H correlation peak [Fig. 2(b)] at ∼1 Å may indicate attempts to form H₂ molecules by detached hydrogens, while the O–O function remains almost intact suggesting that oxygen atoms freed from bonds with hydrogen are not trying to recombine into O₂ molecules at these timescales. Longer times are shown in the  Appendix, and videos in the supplementary material.

These conclusions are confirmed by an analysis of fragments appearing in the simulation box during the simulation. The fragments are identified by an analysis of neighboring atoms within cutoff distances 1.15 Å for H–O, 1.8 Å for O–O and 0.9 Å for H–H bonds (changing these values only trivially alters the fragmentation picture). As it is demonstrated in Fig. 3, the main constituents of the box are atomic hydrogen and oxygen with the moderate presence of OH⁻ and H₂ molecules and the occasional transient appearance of H₃ clusters. The short-lived appearance of such species indicates a liquid flow of hydrogens, which may transiently approach each other forming hydrogen clusters and then drifting apart.

Figure 4 shows that the atomic kinetic temperature increases in the simulations with electronic temperatures elevated above the threshold. As BO approximation does not produce nonadiabatic electron–ion exchange of kinetic energy, this increase in the atomic temperature is the result of the conversion of the excited potential energy to the kinetic energy of atoms. The potential energy released during the bandgap collapse accelerates atoms.⁴¹ Due to this nonthermal acceleration of atoms, the atomic temperature (kinetic temperature) may reach far above the water evaporation point (see Fig. 4). In the absence of the electron–ion kinetic energy exchange within BO approximation, the instantaneous atomic temperature stabilizes at ∼1900 K at electronic temperature T_e ∼ 3 eV.

Having evaluated the threshold electronic temperature within the BO approximation using DFT-MD simulations, we now proceed to the model including nonadiabatic electron–ion coupling. It is done using tight-binding simulation. The quality of the DFTB parameterization mio-1-1 for water under normal conditions was extensively studied in the literature.^42,43 Here, we evaluate its performance at high electronic temperatures. The calculated electronic temperature threshold of the bandgap collapse within the BO approximation is T_e ∼ 3.7 eV (corresponding to 6.1 eV/atom deposited dose, or ∼15% of valence electrons excited to the conduction band). Similar to the described-above DFT results, the collapse of the bandgap starts with the shrinkage of the entire gap and the appearance of the first levels inside of it, induced by detachment of hydrogens from the water molecules (see  Appendix). As the interatomic potential binding hydrogen and oxygen atoms become shallower with an increase of the electronic temperature, at the threshold dose it becomes barely binding. Collisions of the water molecules in the liquid phase at room atomic temperature knock out the first “seed” hydrogen atom, which then knocks out the next one, forming new defect energy levels inside the gap.

With the increase of the dose (electronic temperature), more and more water molecules disintegrate, and eventually, the entire bandgap disappears, forming merged valence and conduction band. Thus, in TBMD simulation, the state also becomes metallic (electronically conducting). The XTANT-3 calculations qualitatively agree with the DFT-MD calculations, even though quantitatively the damage threshold is higher. We thus conclude that the DFTB method can be applied for the qualitative analysis of water at high electronic temperatures. In Sec. III B, we employ it beyond the BO approximation to study nonadiabatic effects and their interplay with the nonthermal damage reported here.

Modeling water response accounting for the nonadiabatic energy exchange between electrons and ions (electron–ion coupling, e.g., scattering) in water noticeably reduces the damage threshold down to the dose of ∼1.1 eV/atom (peak electronic T_e ∼ 1.5 eV, peak density of excited electrons ∼2.8%), as calculated with XTANT-3. The damage takes place via an interplay of the nonthermal and nonadiabatic effects, similar to the reported one in solid oxides and polymers.^18,40 The initial kick to atoms due to nonthermal modification of the interatomic potential leads to a very fast increase in the atomic temperature (as described in Sec. III A). Since electron–ion coupling rate is proportional to the atomic temperature, an increased atomic temperature leads to an increase in the electron–ion coupling, triggering even faster energy exchange (for a detailed discussion, see Refs. 37 and 41). This self-amplifying process increases the atomic temperature at ultrashort, sub-picosecond, timescales, see the example in Fig. 5. The average atomic temperatures reach ∼5000 K after equilibration with the electronic one.

It is interesting to note that the oxygen temperature transiently increases higher than the hydrogen one, indicating that the oxygens experience stronger modifications of the interatomic forces due to electronic excitation to high temperatures. However, the deviations quickly disappear, and within a few hundred femtoseconds, the two species are equilibrated again.

The extremely fast heating leads to molecular fragmentation, which is accompanied by the formation of the defect levels within the bandgap, similar to the BO case in the last section. Here, the bandgap completely collapses at the dose of ∼1.5–1.7 eV/atom deposited to the electronic system (peak electronic temperatures of T_e ∼ 20 000 K = 1.7 eV), see Fig. 6. Water, thus, turns into a metallic liquid state at lower doses than those produced within BO approximation.

Note that this behavior and the electronic temperature coincide very well with the experimental report of the metallic water formation in Ref. 20. As was also noted in Ref. 20, the simulations within BO approximation do not produce conducting state at the experimental conditions (T_e ∼ 20 000 K), which agrees with our estimate in Sec. III A. We thus conclude that the liquid metallic water is produced via an interplay of the nonthermal and nonadiabatic effects, which require simulations beyond the BO approximation.

This process is accompanied by water fragmentation into various species, as shown in Fig. 7. The majority of the produced fragments are hydrogens (H, H₂, with little contribution of H₃), and O and OH radicals, with some contribution of O₂ molecules to which floating hydrogen may attach and detach. With the increase of the electronic dose (electronic temperature) above ∼1.5 eV/atom, the concentration of the OH radical reaches a plateau and new OH radicals stop forming, instead fragmenting into H and O (see  Appendix).

Having evaluated the damage threshold for the formation of metallic water, we now proceed with an estimate of the achieved doses after an SHI impact. For that, we apply the TREKIS-3 Monte Carlo code, which evaluates the ion energy loss in a given matter, secondary electron transport induced, Auger and radiative decays of core holes produces, transport of valence holes and emitted photons.⁴⁴ We estimate the energy density in the ensemble of the excited electrons, and the energy transferred to ions after the impact of the U ion with energy of 2000 MeV, see Fig. 8. Following our previous studies,^41,45 energy transfer to ions is assumed to take place via three channels: elastic scattering of excited electrons, elastic scattering of valence holes, and release of the potential energy of electron–hole pairs associated with the nonthermal acceleration described in Sec. III A above. This way, both nonthermal and nonadiabatic energy transfers may be approximately extracted from a Monte Carlo simulation, as discussed in Refs. 41 and 46.

The results in Fig. 8 suggest that even after the heaviest ion impact with the energy near its Bragg peak (maximal energy loss, S_e ∼ 16 keV/nm), the energy densities reached in the electronic and ionic systems in water are insufficient to produce a metallic state. We thus conclude that the effects of impacts of swift heavy ions in water may be described within the standard models, assuming the production of radicals in a regular insulating water environment.⁹ 

We studied water under ultrafast deposition of energy into its electronic system using two methods: density functional and tight binding molecular dynamics. Both methods applied within the Born–Oppenheimer (BO) approximation predicted the formation of metallic liquid water above a certain dose deposited to the electronic system (electronic temperature ∼3–4 eV, the atomic system starting at room temperature). Including nonadiabatic (non-BO) energy exchange between electrons and ions lowers the threshold dose to ∼1.1 eV/atom. In this case, the collapse of the bandgap and formation of the metallic state is induced by an interplay of the nonthermal and nonadiabatic effects: acceleration (heating) of atoms due to modification of the interatomic potential by the electronic excitation enhances electron–ion coupling. Reaching high ionic temperatures at sub-picosecond timescales, together with high electronic temperatures, also leads to the disintegration of water molecules into various fragments. We show that the formation of OH radicals reaches a plateau at the doses of ∼1.5 eV/atom deposited to the electronic system, above that mostly H and O atomic species are formed. Our results on ultrafast formation of the metallic liquid state in water at electronic temperatures ∼20 000 K are supported by recent experiments.²⁰ 

The supplementary material for this work includes two datasets with the initial atomic positions of water for DFT-MD (32 water molecules) and TBMD (64 water molecules) simulations in the extended XYZ format. The three video files show the time evolution of the pair correlation functions for H–H, H–O, and O–O atoms in water at the electronic temperature of Te ∼ 3 eV up to 500 fs simulated with DFT-MD.

N.M. gratefully acknowledges financial support from the Czech Ministry of Education, Youth and Sports (Grant No. LM2018114). R.V. and A.E.V. were funded by a grant from Russian Science Foundation Grant No. 22-22-00676. Computational resources for XTANT-3 calculations were supplied by the project “e-Infrastruktura CZ” (e-INFRA LM2018140) provided within the program Projects of Large Research, Development and Innovations Infrastructures. DFT calculations have been carried out using computing resources of the federal collective usage center Complex for Simulation and Data Processing for Mega-science Facilities at NRC “Kurchatov Institute,” http://ckp.nrcki.ru/.

The authors have no conflicts to disclose.

Nikita Medvedev: Conceptualization (lead); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (lead); Writing – review & editing (equal). Roman Voronkov: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (supporting). Alexander E. Volkov: Conceptualization (supporting); Funding acquisition (equal); Methodology (supporting); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The distribution of O–H bond lengths in the created initial states, used to set the DFT and TB simulations, are shown in Fig. 9. The peaks of the first nearest O–H and H–H neighbors coincide between the two methods of simulations, DFT and TB, and with the standard values in water.³³ The angular distribution of H–O–H bonds is shown in Fig. 10, where the initial distributions (for DFT and TB) and the one averaged over 0.5 ps of the DFT-MD simulation are demonstrated. The mean angles are close to those simulated with more advanced methods and measured in experiments. Although it is known that PBE functional underperforms in comparison with more advanced ones such as SCAN, we consider that the quality of the produced simulated water is sufficient for this work, which focuses on the extremely high electronic temperatures rather than the fine effects SCAN can reproduce.³³ 

The pair correlation functions for H–H, H–O, and O–O atoms extended to 500 fs simulated with DFT-MD are shown in Fig. 11 (extended Fig. 2); see also movies in the supplementary material. The figure demonstrates that the O–H correlation is essentially lost, indicating the decomposition of H₂O molecules. At the same time, the O–O correlations remain essentially unchanged, as supposed to be in a liquid state. This supports our claim in the main text that the decomposition of water molecules is responsible for the bandgap collapse, and not merely the structure (relative positions of the molecules in the particular state).

The electronic energy levels in water after energy deposition, corresponding to the electronic temperatures of T_e ∼ 3.7 and 4 eV, calculated with TBMD (XTANT-3) in BO approximation are shown in Fig. 12. A qualitative comparison with DFT-MD calculations, shown in Sec. III A (cf. Fig. 1), demonstrates that the mio-1-1 DFTB approximation can be used to describe water at high electronic temperatures. However, overestimation by the absolute value should be kept in mind for comparisons with experiments.

The relative concentrations of various atomic and molecular species in water after irradiation with doses from 0.5 eV/atom to 2 eV/atom calculated with nonadiabatic simulation in XTANT-3 are shown in Fig. 13. The concentration of intact water molecules decreases with the dose increase. The fraction of OH radicals first increases with the increase of the deposited dose, then reaches its maximum at the dose of ∼1.5 eV/atom deposited to the electronic system. At high doses, OH radicals become unstable and dissociate into atomic H and O. Other molecular fragments almost do not form at the timescales studied here.