Possibility of realizing superionic ice VII in external electric fields of planetary bodies.

In a superionic (SI) ice phase, oxygen atoms remain crystallographically ordered while protons become fully diffusive as a result of intramolecular dissociation. Ice VII's importance as a potential candidate for a SI ice phase has been conjectured from anomalous proton diffusivity data. Theoretical studies indicate possible SI prevalence in large-planet mantles (e.g., Uranus and Neptune) and exoplanets. Here, we realize sustainable SI behavior in ice VII by means of externally applied electric fields, using state-of-the-art nonequilibrium ab initio molecular dynamics to witness at first hand the protons' fluid dance through a dipole-ordered ice VII lattice. We point out the possibility of SI ice VII on Venus, in its strong permanent electric field.

INTRODUCTION

The phase diagram of water is exceptionally rich and unusual, because of water molecules’ small size, relatively high molecular symmetry, and, especially, their ability to interact with nearby neighbors via hydrogen bonding. There are currently at least 18 known crystalline ice phases, spanning wide temperature and pressure conditions (–). Most of these represent solid ice structures with various hydrogen bond networks among the molecules sitting in regular crystal lattices. However, as pressure increases, the intramolecular hydrogen interactions with oxygen atoms are weakened, and they start to be shared by oxygen pairs of neighboring molecules, and, eventually, molecules can even dissociate, with protonic diffusions through the crystal lattice. In such a superionic (SI) ice phase, oxygen atoms remain crystallographically ordered while protons become fully diffusive as a result of intramolecular dissociation (–). Ice VII’s importance as a potential candidate for a SI ice phase has been conjectured from anomalous proton diffusivity data (–), although this speculation is controversial and clearcut evidence of proton hopping is elusive.

Molecular simulation would serve as a potentially productive direction to investigate suggested proton hopping and related properties, albeit without success thus far in ice VII (); simulation has found this in ice X at 450 GPa and 1500 to 2000 K (–), with experimental support suggestive at 50 to 60 GPa and at 4750 K (). Still, experimental studies of SI ice are challenging—and especially so in realising static compression using a diamond anvil cell for pressure exceeding 300 GPa (–), even given progress in (overheated) higher-pressure shockwave methods (). Wider theoretical studies indicate possible SI prevalence in large-planet mantles (e.g., Uranus and Neptune) and exoplanets (–); the existence of SI ice under such and similar conditions was proven also very recently by experiment via onset of SI behavior in subnanosecond, laser-ablation “electro-shockwave” propagation in ice (, ). In addition, of great interest for the present study is the observation of electric fields on Venus () and a good likelihood of water thereon (): This may allow for conditions ripe for field-induced SI in ice VII, as we discuss below.

Of the 18 ice polymorphs found thus far up to ~200 GPa, ice VII’s stability across a wide-ranging region above 2 GPa () renders it interesting for potential near-ambient-temperature SI-conducting applications. Ice VII’s simple structure of two interpenetrating, and effectively independent, cubic-ice (Ic) sublattices offers an intriguing “playground” to explore rather unusual electronic behavior of protons moving in the O···O double-minimum potential while under pressure (, ). Anomalous proton diffusivity data in ice VII (–) underpin the conjectured—and exciting—prospect of SI behavior at relatively low pressures (~10 to 15 GPa) compared to more established SI behavior in ice X at much higher (extraterrestrial mantle–like) pressures and temperatures (–, , ).

Achieving—and sustaining (with an intact lattice framework)—lower-pressure SI properties in ice VII in the ~10- to 15-GPa region by an external perturbing agent remains a tantalizing possibility. Given difficulties in confirming experimentally beyond question proton SI-hopping mechanisms as the root of anomalous ice VII proton diffusivity (–) and the failure (so far) of ab initio (AI) molecular simulation to boot (, ), using nonequilibrium AI molecular dynamics (NE-AIMD) to capture in vivid detail the effects of a perturbing agent in overcoming water’s intramolecular-dissociation/proton-transfer free-energy barriers to hop through the lattice and realize SI behavior represents an intriguing challenge. In the present study, we apply external electric fields to do exactly that.

Recently, the Berry phase approach to NE-AIMD in external electric fields has been extended to long-time second-generation Car-Parrinello MD (CPMD); this allows external electric fields to be applied, e.g., to water systems (). Intrinsic electric fields in condensed water phases have been estimated by MD at ~1.5 to 2.5 V/Å (). In liquid water, Saitta et al. () performed nondeterministic CPMD coupled with transition path sampling, to study liquid water and its breakup in (Berry phase) fields, finding a “threshold” dissociation intensity of ~0.35 V/Å. Cassone et al. () studied, using similar NE-CPMD–based methods, ~0.2 V/Å (and higher-intensity) field effects on ice Ih and ice XI mechanical/electric properties. Here, water dissociation and steady-state proton conduction were established above 0.25 and 0.36 V/Å, respectively (ice Ih), with 0.22 V/Å for both in ice XI. However, it has been established by NEMD that for external-field intensities of this order (up to ~15% of intrinsic-field strengths) that high levels of dipolar ordering/alignment is achieved in ices (), so modeling proton-diffusing SI behavior in these ice phases with deterministic NE-AIMD is relevant and timely to enhance our understanding in superionicity’s underlying physicochemical contours. Naturally, ice ionization (no matter how induced) is a quantum process and cannot be described by classical force fields. Therefore, we performed NE-AIMD based on density functional theory (DFT) with state-of-the-art handling of nonlocal dispersion, applied to dipole-ordered ice VII with molecular dipoles oriented along the applied field direction—cf. Fig. 1 (see Materials and Methods for a fuller description). External static electric fields with intensities of 0.02, 0.05, 0.10, 0.25, 0.33, 0.40, and 0.50 V/Å were run for 100 ps of second-generation CPMD at ~10 GPa and 410 K in the canonical (NVT) ensemble (). We computed site-site radial distribution functions (RDFs), infrared (IR), vibrational Raman spectra, and vibrational density of states (VDOS) (cf. Materials and Methods), as well as characterizing proton-hopping events (in terms of hydrogen bond distortion and fraction of molecules involved).

Fig. 1

Trajectory visualization of representative proton-transfer pathways.

Initial structure of simulation supercell (A) containing 128 water molecules with dipoles aligned in direction of the applied static electric field (with a magnitude of 0.33 V/Å). Snapshots from the beginning (B) and from the end (C) of the NEMD run depict the highlighted proton-transfer pathway. The ice structure is enlarged in the blue-circle region and is shown at times (D) 50 ps, (E) 60 ps, and (F) 70 ps of the NEMD simulations, as indicated in (G), showing the system’s potential energy evolution during NEMD. a.u. atomic units.

RESULTS

We did not observe any intramolecular water dissociation or proton transfer in fields with strengths below 0.33 V/Å. However, we did uncover some evidence of a perturbation in structural features, such as a “blurring” in site-site RDFs (fig. S1), increasing levels of strain in intramolecular geometry (fig. S2), widening of dipole-moment distributions (fig. S3), and a broadening distribution of hydrogen bond A-H lengths and A-D-H angles (fig. S4). In short, this is similar to previously observed effects under zero-field conditions of greater pressurization from 5 to 20 GPa, via equilibrium AIMD (). However, at and above 0.33 V/Å, we witness the initiation of Grotthuss-type proton-jump “cascades” as shown in Fig. 1. Although greater amplitude fluctuations in O-H lengths are evident in time-averaged O-H RDFs, with a “spread” by almost 0.1 Å of the intermolecular peak’s “lower tail” to shorter distances in 0.5 V/Å fields relative to the zero-field case (cf. movement from ~1.5 to ~1.4 Å in fig. S1), allowing scope for a certain degree of proton extra-molecular “trespassing” and more substantial intermolecular electron-cloud overlap, the use of NE-AIMD allows the wait-time kinetics for these processes to be gauged efficiently. In Fig. 1, a rather remarkable and concerted “rare-event” proton-jump event occurs after ~60 ps, which would have been previously inaccessible in time scale, vividly illustrating the utility of long-time second-generation CPMD. The proton stretches to greater O-H distances from the “parent” oxygen atom via the electric field’s translational force () pulling the proton away toward the next hydrogen-bonded water molecule along the field direction.

It is energetically less favorable for OH− and H3O+ moieties to persist in these “electrified ices,” as found also by Cassone et al. () (although this plays some role in ice Ih). The comparatively larger water-dissociation threshold here with respect to those in ice Ih and ice XI in () arises from directly using deterministic NE-AIMD, as opposed to biased-sampling (and pragmatic) methods designed explicitly to estimate minimal possible field-intensity thresholds for inducing water breakup. Like in (), we observe a Grotthus-like proton-hopping mechanism, from neighboring oxygen to oxygen atoms, facilitated by dipole-aligned structures. The protons, largely screened by electronic density and bearing total (Hirshfeld) charge +0.25 e in our simulations, move rapidly from one oxygen site to another on a singlet-spin-state potential energy surface, and are stabilized for tens of picoseconds before the next such event occurs (see potential energy plot in Fig. 1). This is not unlike Grotthuss-type proton diffusion in ice XI (), albeit the lesser degree of dipolar orientations due perpendicular alignment of oxygen planes in ice XI renders ice VII protonic “electro-diffusion” more facile. In practice, defects present in real ice samples (such as ice VII or other polymorphs) would be expected to facilitate proton diffusion a great deal, perhaps allowing greater scope for greater H3O+ involvement, although in ice VII, Grotthuss electro-hopping mechanisms dominate.

To underscore the marked change in proton-diffusion behavior, especially at and above 0.33 V/Å field strengths, we present in Figs. 2 and 3 vibrational spectra, time-averaged over the simulations’ durations. The broadening of the librational VDOS mode (cf. Fig. 2A) in more intense fields shows how the greater suppression of librational motion with less scope for rotation-oscillation leads to a more “rigid” hydrogen-bonded hop pathway along which protons conduct superionically via Grotthus jumps. Equally, if not more, marked is the field-induced promotion and suppression of asymmetric (ν3) and symmetric (ν1) O-H stretch bands, respectively (cf. Fig. 2 and especially the Raman spectra in Fig. 3): The reverse-direction “momentum kick” of a departing proton tugged away by the external field in the same direction, by way of an opposite-direction mechanical reaction force in its wake, leaks/resonates into the asymmetric stretch (ν3) mode, which serves to hinder its symmetric (ν1) counterpart; such marked comparative effects are not seen in static-field ν1-ν3 effects on liquid water via NE-AIMD ().

Fig. 2

Vibrational spectra of ice VII under the influence of a static electric field.

(A) VDOS. (B) IR spectrum obtained from the dipole autocorrelation function.

Fig. 3

Raman vibrational spectrum obtained from autocorrelation function of velocities projected onto O─H bond vectors.

The effect of applying an external static electric field is shown in different colors.

DISCUSSION

The essential results, outlined in Figs. 1 to 3, serve to elucidate the subtle intricacies of how electric fields perturb the intermolecular energy landscape in ice VII to realize sustainable SI behavior via a Grotthus-like concerted proton-jump mechanism. Given the recent experimental interest from the quasi-elastic neutron scattering community in inducing potential SI behavior in previously unknown salty ice VII phases (featuring an 18% volume expansion compared to conventional ice VII) (), together with confirmed SI phenomena from electro-shockwaves from nanosecond laser ablation (, ), the present work serves as a powerful “proof of concept,” in that the external electric fields here are similar in intensity to local intrinsic fields in the immediate molecular milieux of salt ions (), which would be expected to promote further water dissociation and proton jumping. In this way, cutting-edge nonequilibrium molecular simulation serves as a convenient test bed to prototype innovative experimental methodologies for realizing lower-pressure SI ice phases with substantially lower field amplitudes, closer to those used in typical experiments (). Given previously mentioned conjectures of SI behavior in (exo-) planets and the confirmed presence of a strong, permanent electric field on Venus with the possibility of water in its atmosphere and more extreme temperature/pressure conditions present that would stabilize ice VII, there seems a good possibility of field-induced superionicity in some planetary ices. We note that ice VII has been suggested to exist in appreciable quantity especially in high-pressure environments, such as Ganymede (, ). At Ganymede, we note very recent reports of electric and magnetic field amplitude enhancements of the order of a million-fold for chorus-mode waves (); with such strong electric fields thereat, coupled with the prevalence of ice VII, Ganymede may also be poised as a ripe location for marked SI behavior (e.g., in ice VII). The present NEMD study suggests the need for a permanent, sustained electric field (such as documented on Venus) to ensure continued forces/momentum for Grotthus-like jumping, in the face of lattice friction and constant energy landscape barriers.

In closing, given the very recent and exciting experimental confirmation of subnanosecond, laser-driven electro-shockwave induction of SI behavior in ice (, ), the present study provides much microscopic insight into the evolution of electric field–driven SI behavior at short time scales and under strong electric field intensities.

MATERIALS AND METHODS

For nonlocal dispersion, the van der Waals density functional (vdW-DF) was used. We used the opt-B88 exchange-correlation part (, ), based on generalized gradient approximation, together with DRSLL nonlocal correlation correction of Dion et al. () and Román-Pérez and Soler () describing the van der Waals interactions. This type of functional was shown to perform exceptionally well for liquid water and ice Ih (–). Recently, we demonstrated its good performance on aqueous system under high pressures on studies of clathrate hydrates and ice VII (, ). Following these works, we used Goedecker-Teter-Hutter pseudopotentials () and the triple-zeta polarized (TZV2P) basis set to describe the wave function. All calculations were done in CP2K software package ().

Analyses of electric field effects on ice VII structure and dynamics requires relatively long simulation times to allow the system relax. Therefore, we applied second-generation CPMD simulation protocol (, ), when the system is propagated in time with a 1-fs time step and the numerical errors caused by insufficient convergence of electronic wave function at each step was corrected by Langevin thermostatting with the γ factor set to 0.005. The system is studied in canonical (NVT) ensemble with 0.1-ns-long productive MD simulations. While the temperature of 410 K is imposed by the thermostat, lattice constant of a cubic 2 × 2 × 2 supercell containing 128 water molecules (see Fig. 1A) was fixed at 12.62 Å, which corresponds to a pressure of 10 GPa with the above-described density functional. Each system was equilibrated carefully by ~15-ps MD runs to ensure that stably propagated MD runs without any considerable energy drift (see fig. S5). Since water-dipole alignment with the applied strong external field is expected, bringing the disordered ice VII structure to its ordered ferroelectric state, before the ionization proceeds, we used the already aligned structure as our computational model to avoid the long-time ordering process.

The homogeneous static electric field was applied, uniform across the supercell, within periodic boundary conditions applied in all three directions, using the Berry phase formulation (). As analogous simulations performed by Saitta and Saija () on liquid water indicated a field intensity of ~0.35 V/Å as the dissociation threshold of the water O─H bond, here, we compare ice VII behavior under field intensities of 0.02, 0.05, 0.10, 0.25, 0.33, 0.40, and 0.50 V/Å, as well as at zero-field conditions. The external-field vector’s direction was aligned with the ordered water molecular dipoles in the ice VII structure. In addition, NEMD simulations, with fields applied perpendicular to this original laboratory direction, were performed as well, confirming the already-found threshold field-strength magnitude.

RDFs g(r) were calculated for all O-O, O-H, and H-H pair interactions in the system to detect structural changes. Further, histograms of water O─H bond lengths, H-O-H valence angles, and individual-water dipoles were collected from the calculated trajectories, and their mean values were calculated from the corresponding normalized distributions 〈x〉 = ∫ xP(x)dx. Hydrogen bonds selected on the basis of geometry criteria (maximum allowed distance between H-bond donor D and acceptor A of 3.5 Å, while the A-D-H angle cannot exceed 35°) were analyzed in analogous way.

The IR spectrum was computed as the power spectrum of the system-collective-dipole autocorrelation functionwhere (t) can be written as sum of individual water dipoles . The individual dipoles were calculated semiclassically from Hirshfeld point charges obtained by DFT. The vibrational density of states (VDOS) was calculated by Fourier transforming the respective velocity autocorrelation function

The autocorrelation function itself was computed on the basis of Wiener-Khinchin relations using the fast Fourier transform algorithm with efficient O(n log n) scaling. To obtain the vibrational Raman spectra, we Fourier-transformed the O─H bond length autocorrelation function. In this way, the vibrational spectrum is decoupled from the translational and rotational modes and is therefore directly comparable with experimentally measured data. All analyses were performed on the last 10 ps of the calculated trajectories, where the system was equilibrated at given external field.