ordering in liquid gallium under extreme

The atomic-scale structure, melting curve, and equation of state of liquid gallium has been measured to high pressure ( p ) and high temperature ( T ) up to 26 GPa and 900 K by in situ synchrotron x-ray diffraction. Ab initio molecular dynamics simulations up to 33.4 GPa and 1000 K are in excellent agreement with the experimental measurements, providing detailed insight at the level of pair distribution functions. The results reveal an absence of dimeric bonding in the liquid state and a continuous increase in average coordination number ¯ n GaGa from 10.4(2) at 0.1 GPa approaching ∼ 12 by 25 GPa. Topological cluster analysis of the simulation trajectories finds increasing fractions of fivefold symmetric and crystalline motifs at high p - T . Although the liquid progressively resembles a hard-sphere structure towards the melting curve, the deviation from this simple description remains large ( ≥ 40% ) across all p - T space, with specific motifs of different geometries strongly correlating with low local two-body excess entropy at high p - T . DOI:

Liquid metals and alloys have exceptional properties that make them particularly attractive for applications: potential uses include electrical energy storage and generation as, e.g., electrodes for all-liquid high capacity batteries [1] and efficient heat exchange fluids in concentrated solar power systems [2]. By virtue of their nontoxicity, low viscosity, and high thermal and electrical conductivity, low-melting point gallium-based liquid metals have applications from cooling integrated electronics to manufacturing flexible and reconfigurable electronic devices and soft robotics [3][4][5][6]. Such optimal thermophysical properties are governed by the atomic-scale structure of these liquids. Knowledge of structural changes and solidification pathways in liquid metals at nonambient pressure (p) and temperature (T) is essential for the development of new materials with novel physical properties and for operating under extreme conditions. Structural information of liquid metals is also key to understanding processes within deep terrestrial and exoplanetary interiors, including metallic core formation [7] and magnetic field generation [8]. While challenging, measuring liquid structure at elevated p-T conditions is a rapidly developing field [9][10][11][12][13][14][15][16].
Gallium is a remarkable metal, exhibiting a rich array of crystal structures at nonambient p-T [17]. At ambient p-T, gallium exhibits an orthorhombic structure (Ga-I) with mixed metallic-covalent bonding featuring Ga 2 dimers [18,19]. This mixed bonding gives rise to unusual characteristics including an anomalously low melting point (T m ¼ 302.9 K [20]) and consequently one of the largest liquid ranges of any element, a 3.2% volume contraction on melting, and a strong tendency for undercooling [20,21]. At elevated p the melting curve exhibits negative dT m =dp up to the I-II-liquid triple point at 1.2 GPa [22] (Fig. 1). The existence of a first-order liquid-liquid phase transition (LLPT) has been postulated [23][24][25][26] on the basis that other candidate polyamorphic liquids exhibit similar anomalous behavior [27], notably water [28], silicon [29], sulphur [30], and phosphorous [9]. Previous in situ structural measurements of liquid gallium at high p are limited to ∼6 GPa [26,[31][32][33][34][35][36][37] by synchrotron x-ray diffraction (SXRD) and 9 GPa by x-ray spectroscopy [38]. At ambient p and T m the average coordination numbern Ga Ga increases from 7 in the solid Ga-I phase to ∼10 [34], compared to a typical value of 11-12 in most liquid metals. A gradual increase inn Ga Ga is observed with increasing p with a simple close-packed liquid predicted by ∼15 GPa [34]. A similar evolution from a complex low-coordinated liquid to a simple liquid at high p has been reported recently in shock compressed tin [39].
In this Letter, we report the atomic-scale structure, melting curve, and equation of state of liquid gallium as measured by in situ SXRD up to 26 GPa, representing a > fourfold increase in the p range compared to previous experimental surveys. Complementary ab initio molecular dynamics (MD) simulations of the liquid atomistic structure were made to 33.4 GPa and 1000 K.
p-T conditions of up to 26 GPa and 900 K were achieved using a membrane driven diamond anvil cell (DAC) with Boehler-Almax anvils (∅500 μm culet) surrounded by a Watlow coiled resistive heater within a vacuum vessel [40,41]. Temperature was measured using a K-type thermocouple attached to one anvil, close to the gasket. To prevent alloying with the gasket, the liquid gallium droplet was loaded into an annulus of dry NaCl (130 μm inner diameter) [42] within the ∅165 μm sample chamber drilled in a pre-indented rhenium gasket. SXRD measurements were made at beam line I15, Diamond Light Source, UK using a Perkin-Elmer 1621 EN flat panel detector and a monochromatic 56 keV x-ray beam collimated by a 20 μm tungsten pinhole to ensure a clean signal from the wholly liquid sample. This offers a significant advantage over laser-heated DAC or shock compression experiments which both suffer from contamination of the liquid signal by diffraction peaks arising from, e.g., thermal insulation media or nonmelted solid [13,39]. Pressure was determined from additional SXRD measurements of the NaCl annulus at each step and the known p-T equation of state of NaCl [52]. We constrained the melting curve by observing either liquid diffuse scattering or Bragg peaks in regular 20 s SXRD acquisitions along a stepped p-T path (Fig. 1). The midpoints of the melting brackets were fitted to a Simon-Glatzel [53] equation modified to force the fit through the I-II-liquid triple point [22] using orthogonal distance regression and yielding a ¼ 7.6ð16Þ and c ¼ 1.39ð18Þ. Longer, 20 min acquisitions were made at ∼2-3 GPa steps in the liquid field just above the melting curve. At the end of the experiment an SXRD measurement of the DAC containing the recovered empty gasket was made to characterize the dominant background component originating from Compton scattering from the diamond anvils. The single crystal diamond reflections were masked prior to integration of the twodimensional diffraction images. Density functional theory (DFT) electronic structure calculations were performed using the Vienna ab initio software package (VASP) [54,55]. The electronic interactions were described by the projector-augmented wave (PAW) [56,57] pseudopotentials with an [Ar] core and 3d 10 4s 2 4p 1 valence electrons. The Perdew-Burke-Ernzerhof (PBE) formulation of the generalized gradient approximation (GGA) exchange correlation functional [58] was used with an energy cutoff value of 500 eV, sampling the Brillouin zone at the Γ point. Molecular dynamics trajectories were calculated in the canonical (NVT) ensemble with N ¼ 600 Ga atoms. Smaller boxes did not accurately reproduce the low-q features observed by experiment. Simulations at various initial volumes (V), were heated at 6000 K for 2 ps, cooled to target temperatures of T ¼ 400, 600, 800, or 1000 K over 2.5 ps, and equilibrated for 15 ps with a simulation time step of 1 fs. An additional simulation at the ambient-p density at T m of 0.0526 Å −3 (6.095 g cm 3 [59]) yields p ¼ 24.38 GPa via the computed stress tensor. This value was subtracted from the computed p at each V-T to correct for the inherent underbinding of the GGA functional. The compute time for a 5 ps simulation interval was approximately 7 days using the University of Bristol BlueCrystal Phase 4 supercomputer on 10 nodes with 28 central processing units (CPUs) per node.
The measured p-dependent structure factors S GaGa ðQÞ and pair distribution functions shown in Figs. 2(a) and 2(b), were obtained by normalizing the background-corrected diffraction patterns using the formalism of Eggert et al. [60], as implemented in our code LiquidDiffract [61]. This iterative procedure [42] exploits the simple behavior of the reduced pair distribution function GðrÞ ¼ −4πn 0 r½g GaGa ðrÞ − 1Þ prior to the first interatomic distance, to eliminate the Q-space manifestations of the unphysical low-r contributions and provide a converged solution for the liquid atomic number density n 0 (Å −3 ) [ Fig. 3(a)]. The ab initio MD g GaGa ðrÞ and corresponding S GaGa ðQÞ functions were computed from the final 5 ps of the simulation trajectories using the R.I.N.G.S. code [62]. The agreement between the S GaGa ðQÞ and g GaGa ðrÞ functions as measured in the DAC at 0.1 GPa (the lowest p) and computed from the ab initio MD simulations at ambient p and 303 K (this study) with previous ambient-p results from neutron and x-ray diffraction measurements [63], is excellent [ Fig. 2. (a)]. We find no evidence from the measured g GaGa ðrÞ or ab initio MD trajectories for short (<2.5 Å) Ga-Ga bonds under any p-T condition, indicating that dimeric bonding characteristic of the Ga-I structure does not persist in the liquid state. The p evolution of the computed S GaGa ðQÞ and g GaGa ðrÞ (Fig. 2) is in good general agreement with the experiment, although with increasing p the experimental reciprocal-space data suffer from increasing statistical uncertainty leading to poorer resolution and stronger Fourier transform artifacts in real space.
The first peak in S GaGa ðQÞ features a pronounced high-Q shoulder that becomes less pronounced by ∼15 GPa in the experimental measurements, matching the p at which liquid Ga is predicted to transform to a hard-sphere-like liquid [34]. However, we note the first peak remains asymmetric in the experimental measurements and a distinct shoulder can be resolved in the simulation results.
A Mie-Gruniesen-Debye thermal equation of state [64] determined from a fit to the high-p-T ab initio MD results, with parameters V 0 ¼ 19.043ð13Þ Å 3 , K ¼ 50.3ð6Þ GPa, K 0 ¼ 4.75ð4Þ, q ¼ −0.09ð21Þ, D ¼ 325 K, g 0 ¼ 2.07ð4Þ, and g ∞ ¼ 0, agrees with the experimentally derived density within the limits of uncertainty [ Fig. 3(a)]. The first peak in g GaGa ðrÞ shifts to smaller radii with increasing p, from r GaGa ¼ 2.79ð2Þ Å measured in the DAC at 0.1 GPa to 2.58(2) Å at 25.9 GPa and 891 K. The development of the average coordination numbern Ga Ga with increasing p, as obtained by integrating over the measured g GaGa ðrÞ or directly from the ab initio MD trajectories with a cutoff value r cut ¼ 3.5 Å, are shown in Fig. 3(b). The experimental and simulation results are in good agreement within the limits of uncertainty revealing a continuous increase in n Ga Ga on densification from ∼10 at ambient p towards closepacked liquid values of ∼12 by 26 GPa.
Recent studies of the local structure of gallium at elevated-p have employed reverse Monte Carlo (RMC) modeling [33,34,36]. RMC is a fitting strategy to generate an atomistic model by minimizing the difference between experimental data (e.g., the pair distribution function) and an input configuration. However, we show that this naive RMC approach can be misleading by comparison with direct analysis of the local structure of the liquid ab initio MD trajectories. This is illustrated in Fig. 4 which compares the local structure obtained by analysis of the ab initio MD result at T ¼ 1000 K, p ¼ 3.2 GPa with new RMC results obtained using two different initial guesses: a disordered configuration (obtained via a linear compression algorithm) and an ordered fcc configuration at the relevant density, constrained by the ab initio MD g GaGa ðrÞ function [42]. In our analysis of the simulation trajectories, we employ two descriptors of local structure: Voronoi indices and the topological cluster classification (TCC) algorithm [65]. In the former case, we count the fraction of particles in polyhedra (nonuniquely) identified by a vector of integers  representing the histogram of the number of edges on the faces. In the latter case, the local environment of the particles is compared to a predefined library of elementary motifs that are important in simple classical liquids, as they minimize the local energy.
For consistency with a recent RMC study of liquid gallium to 1.9 GPa interpreted by Voronoi tessellation [36], we first consider a selection of Voronoi motifs [ Fig. 4(a)]. These results show that naive RMC constrained solely to an input g GaGa ðrÞ not only fails to reproduce the ab initio MD structure but also produces different results depending on the starting configuration. In particular, icosahedral motifs [0,0,12,0], which are marginal in the ab initio MD, are more highly represented in the disordered RMC but entirely absent in the ordered RMC. In fact, the ordered RMC has preserved strong signatures of crystalline order [42]. Similar behavior is observed in the TCC analysis [ Fig. 4(b)]. For simplicity, we present the results for a set fSg of seven specific motifs which describe different types of local environments: tetrahedal ordering (6Z) which is a precursor to crystallization, motifs with four-membered rings (6A,10A,11F), and fivefold symmetric ordering (7A, 10B, 12D). Pentagonal structures such as 7A and 10B are overrepresented in the disordered RMC and absent in the ordered RMC, while crystal-like patterns such as 11F, which have small fractions in the ab intio MD, are overrepresented in the ordered RMC. All these differences emerge despite the naive impression of a good convergence of pair correlations [ Fig. 4(c)].
We continue the analysis focusing on the TCC motifs, due to their relative simplicity of interpretation compared to Voronoi indices. Figure 5 shows how the structural features change as we move from low to high p-T along the melting curve. With increasing T and p the abundance of larger motifs increases. Among these, the fivefold symmetric 10B and the crystalline 11F units stand out, as their abundance almost doubles from low to high p-T. In order to understand to what extent such structural changes result from effective excluded-volume effects, we performed a mapping onto a system of hard spheres in the Percus-Yevick approximation and compare the TCC spectra of ab initio MD with event-driven molecular dynamics [66] for hard spheres [42]. We define a scoring function as a weighted average of the relative deviations of the fractions n i , where the weight w i ¼ s i = P j∈fSg s j is proportional to the number of particles s i in the TCC motif i. Performing this calculation on all the ab initio MD results delivers a contour map of the deviation from hard-sphere behavior (Fig. 6). Although hard-sphere features become progressively more important towards the melting curve, the deviation is always Δ ≥ 40% such that they do not model a sizable part of the emerging structural correlations. Metastable hard-sphere liquids show the formation of low configurational entropy regions [67]. These can be related to the socalled local two-body excess entropy which accounts for fluctuations in the (smoothed) local pair correlationsg i ðrÞ of particle i [42]. This measurement is not an entropy difference, but connects local structural variations to entropic contributions [68][69][70]. Utilizing the ab initio MD configurations, we measure the local excess entropy and compare its distributions for atoms in different local environments. Figure 6 shows that despite broad fluctuations, the 10B and 11F motifs, as well as fivefold symmetric 12D motifs, have significantly lower values of s 2 , suggesting that gallium at high p-T forms regions of exceptionally low configurational entropy, which may help to stabilize the glassy phase beyond the melting curve. From our combined SXRD experimental and ab initio MD simulation approach we have considerably extended the p-T conditions at which the melting curve, equation of state, and nature of local structural ordering in liquid gallium is known. The results reveal an increase in local coordination numbern Ga Ga approaching ∼12 with increasing densification. Analysis of the ab initio MD trajectories reveals the concomitant increase of the number of fivefold symmetric and crystalline motifs at high p-T. Both form regions of low local entropy, a behavior that contrasts with purely repulsive hard spheres, which are dominated by fivefold symmetry. Previous studies predict the development of hard-sphere like behavior in liquid gallium [34] and tin [39]. However, we find that although the local structure progressively resembles that of hard spheres when approaching the melting curve, the deviation from this simple description is always ≥40% across all p-T space. The presence of low configurational entropy motifs in the liquid provides a mechanism for the promotion of metastable polyamorphic phases beyond the high-p melting curve. The emergence of novel amorphous phases from supercooled regimes may be explored in future work using effective potentials checked by liquid structure measurements using heating elements inside the DAC [71] for rapid T quenching at high p.