Velocity autocorrelations across the molecular—atomic fluid transformation in hydrogen under pressure

Non-monotonous changes in velocity autocorrelations across the transformation from molecular to atomic fluid in hydrogen under pressure are studied by ab initio molecular dynamics simulations at the temperature 2500 K. We report diffusion coefficients in a wide range of densities from purely molecular fluid up to metallic atomic fluid phase. An analysis of contributions to the velocity autocorrelation functions from the motion of molecular centers-of-mass, rotational and intramolecular vibrational modes is performed, and a crossover in the vibrational density of intramolecular modes across the transition is discussed.


Introduction
Exotic properties of condensed matter at high pressures like a transition to non-metallic state in Li [1] and Na crystals [2] or emergence of ring-like jump diffusion of atoms in solid iron at conditions of the Earth's inner core [3,4] caused huge interest in simulations and theory of matter at extreme conditions. Liquids at high pressures have their own specific features defined by topological disorder and partial localization of electrons [5,6]. Especially in the field of dynamic properties of liquids at high pressures, there is a lack of understanding the behaviour of collective modes as well as of single-particle dynamics, that is connected with possible transformations in atomic structure at high pressures [7]. The situation with understanding the dynamic properties in molecular liquids which under high pressure undergo dissociation of molecules and even transformation to purely atomic fluid is even more fascinating and less studied. According to the theory of collective dynamics [8] such systems should be considered at least as a mixture of molecular and atomic components. One of the most simple and mysterious (because of the long-time search for its metallic state) molecular systems, which transform under pressure into atomic ones, are solid and liquid hydrogen [9].
Recently we observed in ab initio simulations of hydrogen fluid an interesting behaviour of the long-wavelength charge fluctuations, which made evidence, that right in the transition region from purely molecular to atomic hydrogen fluid there can appear unscreened ions with long-range Coulomb interaction [10]. This is in line with the recent results obtained with Coupled Electron-ion Monte Carlo simulations focused on detection of the possible formation of stable molecular ions at the liquid-liquid phase transition in hydrogen [11]. Hence, the molecular-atomic fluid transformation region can be much more sophisticated than just simply a mixture of neutral molecules and atoms, that must have some fingerprints in the dynamic properties of these fluids.
Actually, the dynamic properties of hydrogen fluid at high pressures were not studied enough. Mainly the old experimental [12] and simulation [13] studies of the sound propagation in hydrogen fluid are well known. Having in mind the study of collective dynamics, which is extremely sensitive to the existence of correlations between different species as well as to the presence of ions with long-range Coulomb interaction [14,15], we will first study the single-particle dynamics across the transition region. In particular, we aimed at the study of velocity autocorrelation functions for hydrogen fluids in a wide range of densities: from molecular liquid up to the metallic atomic liquid. The remaining paper is organized as follows: in the next section we report the details of our ab initio simulations, then will report our results for the velocity autocorrelation functions, their Fourier-spectra and diffusion coefficients, and the conclusion of this study will be given in the last section.

Ab initio molecular dynamics simulations
First-principles simulations for seven densities of hydrogen fluid were performed using a model system of 1000 particles at T = 2500 K by the VASP package. The chosen temperature is high enough for application of classical equation of motions to protons. The starting configuration for the smallest density was taken from the classical simulations of molecular dimers. The standard NVT ensamble was used for each density. Upon an equilibration of up to 8 000 timesteps, we checked out the temperature and pressure fluctuations and then switched to production runs, each of 20 000 timesteps. The time step in ab initio simulations was 0.2 fs. The range of simulated densities corresponded to the pressures from 2.5 GPa (purely molecular fluid) up to 166.5 GPa (metallic atomic hydrogen fluid [10]).
The electron-ion interaction was the projector-augmented wave (PAW) potential [16,17], and the plane-wave expansion cut-off was 250 eV. The generalized gradient approximation in PBE version [18] was used in account for exchange-correlation effects. In the sampling of the Brillouin zone we used only the Γ point because of the large size of our simulated systems.

Results and discussion
The main features in the behaviour of the structural and dynamic properties of hydrogen fluid in a wide range of pressures are defined by the region of transformation from purely molecular to purely atomic fluid, in which the hydrogen molecules break up due to the pressure increase with simultaneous increase of a contribution from purely atomic subsystem. In figure 1 (a) we show the density dependence of the velocity autocorrelation functions (VACF) which estimates the time correlation in proton velocity along the trajectory, averaged over all the particles and over all the possible origins of the calculated time correlations t 0 . In order to analyze the changes in VACFs, we plotted in figure 1 (b) the pair distribution functions g(r) for four densities. They make evidence of the pure molecular fluid at the lowest density 0.2847 g/cm 3 with the intramolecular peak of g(r) centered at 0.75 Å, while at the highest studied density of 0.9610 g/cm 3 , the g(r) looks typical as for purely low-density atomic fluids. For intermediate densities, one observes a gradual reduction of the intramolecular peak and a transformation (in amplitude and with a shift towards the smaller distances) of the intermolecular peak into the first interatomic peak (at high densities). Here, we should mention that although the classical equations of motion of protons were considered, much more precise pair distribution functions and in general a picture of the structural features across the transformation from molecular to atomic hydrogen fluid are observed in the Quantum Monte Carlo simulations [9]. The typical time dependence of VACFs for dense fluids lies in the existence of the so-called cageeffect, which is reflected by the negative region of ψ(t) i.e., when the direction of particle velocity at time t is opposite (or rather almost opposite in order to result in the negative scalar product) with the initial velocity v i (0) due to back-scattering on the nearest neighbours. In a wide range of densities, we do observe in figure 1 (a) the cage-effect. As a consequence of the molecular structure, we observed a combination of the high-frequency oscillations and the cage effect in the velocity autocorrelations, which are very pronounced, especially for the two lowest studied densities [see figure 1 (a)]. For density 0.4920 g/cm 3 , the high-frequency oscillations are smeared out, while for the higher densities they are not observed. For the densities up to 0.6052 g/cm 3 , the back-scattering effect is clearly observed, though for the higher densities studied here the negative region of ψ(t) vanishes and the VACFs show the time decay typical of gases. For the highest studied density 0.9610 g/cm 3 , the VACF again becomes nonmonotonously decaying in time showing a weak oscillation at time ∼ 0.015 ps, which with the further increase of density (and pressure) would perhaps increase in amplitude and should definitely again lead to the back-scattering effect at ultra-high pressures.
We calculated the Fourier-spectra of VACFs (figure 2) in order to look how the low-frequency (mainly acoustic excitations and diffusion) processes and high-frequency modes change across the transformation from molecular to atomic hydrogen fluid. For the low densities up to 0.4920 g/cm 3 , the high-frequency intramolecular modes are well-defined in the spectra. Note that due to the rotational-vibrational coupling, one observes a splitting, i.e., two peaks in the region of frequencies 600-900 ps −1 . For higher densities, the Fourier-spectra of VACFs look monotonously decaying functions in the high-frequency region, and the well-defined peak in the region of frequencies ∼ 180 ps −1 disappeares too. Herein below we will try to decompose the VACFs into contributions from the motion of center-of-mass, rotations and intramolecular modes in order to assign the different frequency regions to these processes. However, we previously note that the Fourier-spectrum at zero frequency, which actually is directly connected to the diffusivity via the Kubo integral, changes non-monotonously with density.
In figure 3, we show the density dependence of the diffusion coefficient D calculated via Kubo integral (plus symbols with error bars), and its check via calculations of the long-time asymptotes of mean-square displacements (MSD) R 2 (t) [19,20] (cross symbols in figure 3). Both methods of the estimation of diffusivity result in perfect agreement (figure 3). As it was seen from the Fourier-spectra at zero-frequency in figure 2, the density dependence of the diffusion coefficient shows three regions: decay in the region of densities up to 0.4920 g/cm 3 due to reduction of the free volume in mainly molecular fluids, further increase of D due to the rapid increase of the free volume because of the break-up of hydrogen molecules, while for densities higher than 0.7558 g/cm 3 , again there is observed a decay of D due to the reduction of the free volume in mainly atomic fluid. A relation between the free volume and diffusion coefficient was proposed by Bylander and Kleinman [21] and was used in [7] to explane the pressure dependence of diffusivity in liquid Rb across the structural transformation under pressure. To assign different peaks in the Fourier-spectrum of VACFs, we decomposed it into the contributions due to translational motion of the centers-of-mass, rotational motion around the centers-of-mass, and the intra-molecular stretching vibrational normal mode. Such a decomposition is quite straightforward in the molecular frame adopting the vibrational normal mode to oscillations along the molecular axis, while the rotation corresponds to the projection of the instantaneous proton velocity onto an orthogonal axis to the molecular one. In figure 4, we show this decomposition for the case of the purely molecular fluid at the lowest density. The two high-frequency peaks are due to the intra-molecular stretching normal mode split into two as a result of the coupling to the rotational mode of the molecule. This can be proved by removing the rotational motion that results in a single peak at 780 ps −1 , as is shown in the inset of figure 4. The low-frequency peak consists of contributions both from the rotation of the molecule and the translational motion, the latter extending to very low frequencies and being responsible for diffusivity. For higher densities, we performed similar decompositions when the two neighbour protons were separated by distances not larger than 1.1 Å.
In figure 5, we show the power spectrum of the intra-molecular stretching normal mode at different densities. An increase of density first leads to a small reduction of the peak position, then to a rapid  decrease in the frequency, and eventually to a transformation into a wide peak at 470 ps −1 , corresponding to the short-wavelength mode in the atomic fluid. Actually, these three regimes are consistent with the non-monotonous change of the diffusion (in figure 3) across the transformation from purely molecular to atomic hydrogen fluid. The maximum position of the normal mode distribution shows a decay from the state of molecular fluid with well-defined intramolecular normal modes towards very smeared-out distributions of short-wavelength modes in the atomic fluid with practically the same amplitude. This disappearance of the normal mode can be considered as some kind of the order parameter for the observed transition from molecular to atomic fluid in hydrogen under pressure.

Conclusions
Ab initio molecular dynamics simulations were performed for seven densities of the hydrogen fluid at 2500 K with the purpose to estimate the features in velocity autocorrelations across the transformation from purely molecular to atomic fluid. We observed a non-monotonous density dependence of diffusion coefficients across this transformation, which shows a region of densities with an increase of the diffusivity due to the break-up of hydrogen molecules. The velocity autocorrelation functions reveal changes from the time dependence with high-frequency oscillations due to intramolecular normal modes to the typical monotonously decaying one for low-density atomic fluids, which, however, for the highest studied density (pressure) again shows a non-monotonous time decay. The Fourier spectra of velocity autocorrelation functions were analyzed using a decomposition into contributions from translational and rotational motions as well as from normal intramolecular modes. The density dependence of contribution from normal modes to vibrational spectra shows their gradual disappearance and transformation of the corresponding high-frequency peak at high density into the one coming from the short-wavelength modes involving the most nearest neighbours in the purely atomic fluid.