# Dynamical and transport properties of liquid gallium at high pressures

D. Sheppard, S. Mazevet, F. Cherne, R. Albers, K. Kadau, T. Germann, J.

Kress, L. Collins

► To cite this version:

D. Sheppard, S. Mazevet, F. Cherne, R. Albers, K. Kadau, et al. Dynamical and transport properties of liquid gallium at high pressures. Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, 2015, 91 (6), 10.1103/PhysRevE.91.063101. [hal-02483246](https://hal.science/hal-02483246)

HAL Id: hal-02483246

https://hal.science/hal-02483246

Submitted on 15 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Dynamical and transport properties of liquid gallium at high pressures

D. Sheppard,<sup>1,*</sup> S. Mazevet,<sup>2,3</sup> F. J. Cherne,<sup>1</sup> R. C. Albers,<sup>1</sup> K. Kadau,<sup>4</sup> T. C. Germann,<sup>1</sup> J. D. Kress,<sup>1</sup> and L. A. Collins<sup>1</sup>

<sup>1</sup>Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA

<sup>2</sup>CEA, DAM, DIF, 91287 Arapajon Cedex, France

<sup>3</sup>LUTH UMR 8102, Observatoire de Paris, CNRS Université Paris Diderot, 5 place Janssen, 92195 Meudon Cedex, France

<sup>4</sup>Siemens Energy Inc., Charlotte, North Carolina, USA

(Received 27 August 2014; published 1 June 2015)

Quantum molecular dynamics (QMD) simulations are used to calculate the equation of state, structure, and transport properties of liquid gallium along the principal shock Hugoniot. The calculated Hugoniot is in very good agreement with experimental data up to a pressure of 150 GPa as well as with our earlier classical molecular dynamics calculations using a modified embedded atom method (MEAM) potential. The self-diffusion and viscosity calculated using QMD agree with experimental measurements better than the MEAM results, which we attribute to capturing the complexity of the electronic structure at elevated temperatures. Calculations of the DC conductivity were performed around the Hugoniot. Above a density of 7.5 g/cm³, the temperature increases rapidly along the Hugoniot, and the optical conductivity decreases, indicating simple liquid metal behavior.

DOI: 10.1103/PhysRevE.91.063101

PACS number(s): 52.25.—b, 02.70.—c, 52.65.—y, 52.50.—b

## I. INTRODUCTION

Gallium is a complex element with four equilibrium and several metastable allotropic phases known to exist [1–3]. The low-pressure A11 equilibrium phase has a relatively complex, partially covalent molecular crystal structure that retains some character into the metallic liquid [4–7]. In marked contrast to more typical elements such as hydrogen, nitrogen, and oxygen, where the closure of the band gap with pressure has been associated with the disappearance of covalent bonds [8–10], gallium remains (with iodine) among the few elements where the existence of covalent bonds has been experimentally established in the metallic A11 phase.

While liquid Ga has been extensively studied both experimentally and theoretically, measurements at high pressures have thus far been limited due to its low melting temperature of 303 K. Experimental measurements of the principal shock Hugoniot (the locus of final states accessible by a single shock compression wave) were reported in 2001 by Fritz and Carter [11]. In this paper, we use quantum molecular dynamics (QMD) simulations to study the evolution of the dynamical, structural, and transport properties of liquid gallium along the principal shock Hugoniot.

In quantum molecular dynamics [12] (QMD), also called *ab initio* molecular dynamics, the electrons receive a quantum mechanical treatment within the context of finite temperature density functional theory (FT-DFT), while the ions are propagated classically on the corresponding Born-Oppenheimer surface. Albeit computationally expensive, this method has been extensively used in the past few years to predict the properties of various systems and for a broad range of extreme conditions characterized by large temperature, pressure, and/or radiation field [8,13,14]. QMD simulations can be combined with linear response theory, which is particularly attractive, as it provides a consistent set of dynamical, electrical, and optical properties from the same simulations and for a broad range of densities and temperatures [15,16].

For gallium, this method has previously been applied [17] to study the properties of the liquid phase at ambient density (6.09 g/cm³) for temperatures below 1000 K. Using a limited number of atoms (e.g., 64) in the simulation supercell and a limited *k*-point grid, the method reproduced the structure of the liquid and the experimental DC conductivity while the experimental diffusion coefficient was underestimated by a factor of two. The present study builds on this previous work by reporting simulations where both the convergence of supercell and *k*-point grid sizes were investigated, and by simulations performed along the principal Hugoniot up to a pressure of 150 GPa. We find that the calculated Hugoniot is in very good agreement with both the experimental data [11] and recent classical molecular dynamics calculations [18,19] performed using a modified embedded atom method (MEAM) potential [1]. Whereas the liquid diffusivity and viscosity obtained with the classical MEAM approach departs significantly from the experimental data, the present QMD simulations are in better agreement. Finally, the DC optical conductivity increases with increasing density below 7.5 g/cm³ (~20 GPa) due to increasing structural correlations, but then decreases at higher densities and pressures as a simple liquid metal behavior is recovered.

## II. METHOD AND COMPUTATIONAL DETAILS

### A. Quantum molecular dynamics in VASP

To calculate the Ga Hugoniot, we performed QMD simulations using the Vienna *ab initio* simulation package (VASP), a plane-wave pseudopotential code developed at the University of Vienna [20]. We used periodic supercells with up to 500 atoms initially arranged in an FCC crystal structure. Density functional theory calculations were preformed using the generalized gradient approximation with the PW91 functional [21]. We used a Vanderbilt ultrasoft pseudopotential [22,23] with a plane-wave cutoff energy of 300 eV to perform the molecular dynamics simulations. The value chosen was checked by performing simulations using a higher cutoff energy. With this

<footnote>*danielsheppard@lanl.gov</footnote>
plane-wave cutoff, the resulting internal energy and pressure converged to 0.01 eV/atom and 0.1 GPa, respectively.

To determine the principal shock Hugoniot, we performed isokinetic simulations where the temperature is fixed by a direct rescaling of the velocities. Within this ensemble, the volume and particle number are also constant while the temperature is chosen in order to span the Hugoniot points. We recall that the Rankine-Hugoniot equation [24] relates the specific volume $\bar{V} = 1/\rho$, internal energy $U$, and pressure $P$ at initial and final states (subscript 1 and 2, respectively) as

$$ (U_1 - U_2) + \frac{1}{2}(\bar{V}_1 - \bar{V}_2)(P_1 + P_2) = 0. \quad (1) $$

The shock Hugoniot, or adiabat, is the locus of final conditions {$\bar{V}_2, U_2$, and $P_2$} that satisfy this equation. The internal energy per atom $U = \frac{3}{2} k_B T + \frac{E}{N}$ is computed as the sum of the ion-kinetic energy and the time average of the DFT potential energy $E$, where $k_B$ is the Boltzmann constant and $N$ is the number of atoms in the simulation supercell. The pressure $P$ is computed in a similar fashion using the Hellmann-Feynman forces derived from the DFT potential energy. The simulations were advanced using time steps of 3 fs for a total simulation time of 15–25 ps. The resulting steady-state properties were least-square fitted to quadratic functions in $P$ and $T$ for each isochore to obtain the Hugoniot points that satisfy Eq. (1).

## B. Dynamical properties

The self-diffusion coefficient $D$ is computed from the trajectory after the system has equilibrated. This is calculated via two methods: by mean-square displacement (MSD),

$$ D = \frac{1}{6t} |\langle \mathbf{R}_i(t) - \mathbf{R}_i(0) \rangle|^2, \quad (2) $$

and by the velocity autocorrelation function,

$$ D = \frac{1}{3} \int_{0}^{\infty} \langle \mathbf{V}_{i}(t) \cdot \mathbf{V}_{i}(0) \rangle dt, \quad (3) $$

where $\mathbf{R}_i$ is the position of atom $i$, $\mathbf{V}_i$ is the velocity of atom $i$, and $t$ is the time. Brackets indicate statistical summations over all atoms and time origins. In this case the interval between time origins is $\Delta t$ so every step serves as a time origin.

The shear-viscosity ($\eta$) for the QMD runs is calculated directly from the autocorrelation function of the off-diagonal components of the pressure tensor [25],

$$ \eta = \frac{V}{k_B T} \int_0^\infty \langle \mathbf{P}_{\alpha\beta}(t) \cdot \mathbf{P}_{\alpha\beta}(0) \rangle dt. \quad (4) $$

To improve statistical precision we average over the five independent components of the pressure tensor, $\mathbf{P}_{xy}$, $\mathbf{P}_{yz}$, $\mathbf{P}_{zx}$, $(\mathbf{P}_{xx}-\mathbf{P}_{yy})/2$, and $(\mathbf{P}_{yy}-\mathbf{P}_{zz})/2$.

To compare $\eta$ and $D$ to experimental numbers where only the viscosity or diffusivity are given, we use the empirical formula established for several metals by Chisholm *et al.* [26],

$$ \frac{D\eta b}{k_B T} = c, \quad (5) $$

where $c = 0.18$ for Ga and $b = V_A^{1/3}$, with $V_A$ the atomic volume.

## C. Optical properties

Following each QMD simulation, a total of five uncorrelated configurations were selected from an equilibrated portion of the molecular dynamics trajectory. For each configuration, we calculate the optical conductivity using the Kubo-Greenwood formulation [27–29]. This gives the real part of the electrical conductivity, as a function of frequency $\omega$ and at a particular $k$ point, as

$$ \sigma_1(k, \omega) = \frac{2\pi}{3\omega\Omega} \sum_{j=1}^{n_b} \sum_{i=1}^{n_b} \sum_{\alpha=1}^{3} [F(\epsilon_{i,k}) - F(\epsilon_{j,k})] \times \langle \Psi_{j,k} | \nabla_{\alpha} | \Psi_{i,k} \rangle^2 \delta(\epsilon_{j,k} - \epsilon_{i,k} - \omega). \quad (6) $$

In Eq. (6), we use atomic units, with the electron charge $e$, Planck’s constant $\hbar$, and the electron mass $m_e$ all set to unity. The $i$ and $j$ summations range over the $n_b$ discrete bands (orbitals) included in the triply periodic calculation for the cubic supercell volume element $\Omega$. The $\alpha$ sum is over the three spatial directions and improves the statistics. $F(\epsilon_{i,k})$ is the Fermi weight corresponding to the energy $\epsilon_{i,k}$ for the $i$th band at $k$. In practice, because of the finite simulation volume and resulting discrete eigenvalues, the $\delta$ function must be broadened [15,16]. We use a Gaussian broadening that is as small as feasible without recovering the local oscillations in the optical conductivity resulting from the discrete band structure. We use a PAW (projected augmented wave) [30] pseudopotential to calculate the matrix elements appearing in Eq. (6). The PAW formulation, which is equivalent to an all-electron calculation, provides a simple formulation of the matrix elements without the nonlocal terms usually appearing in other pseudopotential formulations. The simulation parameters were otherwise kept constant, and the use of two different pseudopotentials for the dynamical and optical calculations was solely guided by performance considerations since the resulting dynamical properties were otherwise identical.

# III. RESULTS AND DISCUSSION

## A. Dynamical properties along the shock Hugoniot

As initial conditions for the cohesive energy $U_1$ and the atomic volume $\bar{V}_1 = 1/\rho_1$, we considered both the experimental and the calculated values of the A11 structure reported previously [1].

Figure 1 shows the effect of the QMD simulation parameters on the resulting principal Hugoniot. We first consider the influence of the reference state ($U_1$, $V_1$, $P_1$) in Eq. (1). The principal Hugoniot calculated is similar whether we use the experimental or calculated values for the A11 cohesive energy and density. The most noticeable influence is observed at the highest densities investigated in this work. When using the experimental reference state, $U_1 = -2.97$ eV/atom; $\rho_1 = 5.94$ g/cm$^3$, the Hugoniot point found at a density of 10.0 g/cm$^3$ is about 7% lower than when using the calculated GGA values, $U_1 = -2.94$ eV/atom; $\rho_1 = 5.7478$ g/cm$^3$. In either case, the 54-atom supercell gives a good description of the variation of experimental pressure along the Hugoniot. Consequently, we will use the GGA reference state as initial conditions for consistency, despite the slight improvement
<figure><img src="image_2.png" /><figcaption>FIG. 1. (Color online) Principal Hugoniot of liquid Ga obtained by experiments (Fritz [11]) and the present QMD simulations. The results are displayed as a function of (left) density and (right) particle velocity $u_p$. The QMD simulations involved 54, 108, or 256 particles in the simulation cell, and exp and GGA denote which reference state is used (see text).</figcaption></figure>

shown in Fig. 1(b), where the pressure is plotted as a function of particle velocity $u_p$. Here we see that the experimental reference state gives a slightly stiffer Hugoniot, which is in better agreement with the experimental data from Fritz and Carter [11]. We also note that increasing the supercell size to 108 or 256 atoms does not significantly improve the agreement with the experimental data, suggesting that a 54-atom unit cell is sufficient for calculating Hugoniot pressure-volume relationships.

However, simulation supercell size effects make a significant difference when considering the temperature variation predicted along the Hugoniot. Figure 2 compares the temperatures obtained for three supercells. The temperatures along the Hugoniot differ by several thousand Kelvin. We attribute this to long-wavelength correlations present in Ga, which require larger supercells to capture, and the internal energy per atom that results. The three calculations are performed at the $\Gamma$ point, and the convergence in the electronic structure is obtained by simply increasing the number of atoms. Increasing the number of $k$ points for the 108-atom supercell does not significantly alter the internal energies or the pressure calculated, suggesting that size effects only play a large role when predicting Hugoniot temperatures. Previous calculations of Ga liquid properties only used 64 atoms in the simulation cell [4,5]. Despite the slow convergence of the predicted temperature along the Hugoniot the experimentally observable properties (pressure, particle speed, and shock speed) are well converged by 256 atom simulations.

Figure 2 also shows a direct comparison with the results of classical molecular dynamics simulations performed using the MEAM potential developed by Baskes et al. for Ga [1]. The MEAM calculations were performed using 1372 atoms starting in a liquid configuration that had been equilibrated at T = 308 K and zero pressure. The calculated equilibrium den-

<figure><img src="image_3.png" /><figcaption>FIG. 2. (Color online) Variation of the temperature as a function of pressure along the principal Hugoniot of liquid Ga using either QMD or classical MEAM [18] molecular dynamics simulations.</figcaption></figure>

sity was 5.983 g/cm$^3$, slightly below the experimental density of 6.078 g/cm$^3$. From this configuration, the uniaxial Hugoniot technique developed by Maillet et al. [31] was used to reach the equilibrium Hugoniot state for a specified final volume $V_2$ by equilibrating for 50 ps using the Hugoniot equations of motion, and averaging the system properties (pressure and temperature) for the final 40 ps. The Hugoniot was determined by taking incremental compression steps of 2%. Calculations were performed to determine the effect of system size, with no difference detected between samples of 1372 and 10976 atoms. While this potential was adjusted with input from ab initio calculations similar to those performed here, it is still remarkable to notice the level of agreement in the prediction of the principal Hugoniot points for the 108-atom supercell, even at the higher densities where the electronic structure significantly departs from the ambient conditions where the MEAM was fit. At the higher Hugoniot points, temperature effects start to play a significant role, but the classical MEAM potential still provides a surprisingly accurate description.

This is further seen in Fig. 3, where we compare the particle and shock velocities, as well as pressure, along the Hugoniot as computed using the classical MEAM potential [18] and the current QMD calculations. The MEAM Hugoniot is somewhat softer (by 10%) than both experiment and QMD results at the highest pressures considered here; this is most clearly seen in Fig. 3(c) and may also be contributing to the temperature discrepancy seen in Fig. 2. As noted previously [18], this softening can be corrected by adjusting a cubic anharmonicity factor $\delta$, which represents the pressure derivative of the bulk modulus in the Rose universal equation of state [32]:

$$E^u(\mathbf{R}) = -E_c \left( 1 * a^* + \frac{r_e}{\mathbf{R}} \delta a^{*3} \right) e^{-a^*}, \quad (7)$$
<figure><img src="image_4.png" /><figcaption>FIG. 3. (Color online) Comparison between the QMD (circles) and classical MEAM [18] (solids lines) molecular dynamics simulations, and experimental measurements [11] (triangles), of the principal Hugoniot of liquid Ga.</figcaption></figure>

with

$$a^* = \alpha \left( \frac{\mathbf{R}}{r_e} - 1 \right) \quad \text{and} \quad \alpha^2 = \frac{9\bar{V}B}{E_c}, \qquad (8)$$

where $E_c$, $r_e$, $\bar{V}$, and $B$ are the cohesive energy, nearest-neighbor distance, atomic volume, and bulk modulus, respectively, all evaluated at equilibrium in the reference structure. $E^u(\mathbf{R})$ is the energy per atom of the reference structure as a function of nearest-neighbor distance $R$. For the case of Ga, the FCC structure was chosen as the reference structure and initial fit of the MEAM potential to ab initio calculations was performed [1]. This original Ga MEAM potential resulted in a cubic anharmonicity $\delta = 0.097$ to obtain better agreement with the cold compression curve and the pressure derivative of the bulk modulus predicted using ab initio molecular dynamics [1], although a modified $\delta = 0.127$ was subsequently suggested from the shock Hugoniot [18]. This simple modification extends the range of validity of the MEAM potential and enables calculation of the thermodynamical properties up to 150 GPa, as shown in Fig. 3(c) for both values of this parameter.

Figure 3(a) shows that for low particle velocities, the MEAM variation of the shock velocity as a function of particle velocity follows the ab initio results toward lower values than the experimental results. This can be directly associated to the different equilibrium volume obtained using GGA for the A11 structure. As reported in Ref. [18], the GGA equilibrium volume for the A11 structure is $\bar{V}^{\text{GGA}} = 20.15 \text{ Å}^3/\text{atom}$, about 3% higher than the experimental $\bar{V}^{\text{GGA}} = 19.47 \text{ Å}^3/\text{atom}$.

For the shock and particle velocities, Fig. 3(a) shows a very good agreement with the ab initio results for all the conditions along the principal Hugoniot studied here. Both values of the $\delta$ parameter tested lead to rather similar Hugoniot conditions in either $P$ or $u_s$ versus $u_p$ space, with good agreement between the MEAM calculations, the ab initio calculations, and the experimental data.

## B. Liquid structural and transport properties

The structural properties of the liquid are examined via the pair correlation function $g(r)$. Figure 4 shows the correlation function and the mean-square displacement at ambient density ($\rho = 6 \text{ g/cm}^3$) and a temperature $T = 1000 \text{ K}$ for which both experimental data [33] and previous ab initio calculations [17] are available. The pair correlation function agrees nicely with previous ab initio and MEAM calculations. Comparisons of $g(\mathbf{r})$ for MEAM and QMD along the Hugoniot will be addressed later. Figure 4(a) shows that $g(r)$ converges rather rapidly with respect to number of particles, but this is less so for the mean-square displacement, which can be used to calculate the diffusion coefficient from Eq. (2). A $\Gamma$-point calculation with a 54-atom supercell does not lead to a converged mean-square displacement. Also when calculating the trajectory with the $\Gamma$ point only and without scaling the density as performed in Ref. [17], the resulting pressure is negative by a few GPa. The MSD displacement converges slowly with system size compared to the convergence of the pressure. We attribute this slow convergence of the MSD to the fact that a limited simulation does not capture the complexity of the Ga electronic structure and thus does not lead to the correct metallic state for the system. The simulation performed with 256 atoms in the simulation cell leads to

<figure><img src="image_5.png" /><figcaption>FIG. 4. (Color online) (a) ab initio pair correlation function $g(r)$ (b) and mean-square displacement for liquid Ga at normal density and for a temperature of $T = 1000 \text{ K}$. The number of atoms in the simulation cell is indicated in the graph.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 5. (Color online) (a) Velocity auto-correlation and (b) Fourier transform for liquid Ga at normal density and for a temperature of $T = 1000$ K. The number of atoms in the simulation cell is indicated in the graph.</figcaption></figure>

a diffusion coefficient $D = 8.66 \times 10^{-5} \text{ cm}^2\text{s}^{-1}$ compared to $D = 8.35 \times 10^{-5} \text{ cm}^2\text{s}^{-1}$ for the 500-atom cell, which corresponds to a 4% difference. Given this difference we feel using the 256-atom cell is an acceptable cost-accuracy tradeoff for calculating diffusion constants. The much lower diffusivity obtained by Holender et al. [17], $D = 6.5 \times 10^{-5} \text{ cm}^2\text{s}^{-1}$ at a temperature of $T = 982$ K and an effective density of $\rho = 5.89 \text{ g/cm}^3$, was obtained using only 64 atoms in the simulation cell. The MSD shown in Fig. 4(b) suggests that these earlier results were insufficiently converged with respect to system size.

This convergence with respect to system size can also be seen in the velocity auto-correlation function $VAC(t) = \langle \mathbf{V}_i(t) \cdot \mathbf{V}_i(0) \rangle$ in Fig. 5(a), which appears to indicate convergence by 108 atoms. However, the Fourier transform of the VAC, Fig. 5(b), demonstrates that convergence of the low-frequency behavior requires 256 atoms, as with the MSD.

Figure 6 shows a comparison of the *ab initio* calculations with the experimental data [34,35] and the classical results obtained using the MEAM potential [1]. We performed calculations of diffusion and viscosity using Eqs. (2)–(4) for ($\rho = 5.76 \text{ g/cm}^3$, $T = 823$ K), and ($\rho = 5.55 \text{ g/cm}^3$, $T = 1163$ K). These two points correspond to conditions where experimental measurements of the viscosity have been performed. A third point at ($\rho = 6.0 \text{ g/cm}^3$, $T = 1000$ K) is included to illustrate the sensitivity of diffusion and viscosity to pressure and density. Equation (5) was used to convert the Tippelskirch’s viscosity measurements [35] to diffusivity data. The strong overlap of direct measurements of the diffusivity from Riedl [34] with those converted from Tippelskirch’s data give us confidence in using Eq. (5) as a vehicle for comparison. Figure 6 shows that Ga is more mobile than predicted from MEAM calculations, and these

<figure><img src="image_7.png" /><figcaption>FIG. 6. (Color online) Variation of the (a) density, (b) diffusivity, and (c) viscosity as a function of temperature for liquid Ga. Experimental results: Riedl (filled red triangles) [34]; Tippelskirch (empty blue triangles) [35]. Simulation: MEAM (green squares); QMD (black circles).</figcaption></figure>

present QMD results are in better agreement with experimental data. The MEAM calculations deviate from the experimental data as the temperature approaches the melting point. The low-temperature behavior is not completely surprising as one expects that the mismatch in the equilibrium volume would have a greater impact upon the dynamical properties as one approaches melting.

This discrepancy in diffusivity between quantum and classical MD simulations and experimental data carries over to the viscosity. The MEAM potential was developed based on *ab initio* calculations of *static* equilibrium properties such as lattice constant, cohesive energy, and bulk modulus. It is not surprising that it cannot accurately capture *ab initio* results for dynamical quantities, especially at low temperatures where small energy differences in the barrier for diffusion can lead to large differences in the rate.

We now turn to the evolution of the structural properties along the Hugoniot. All simulations were performed as close to Hugoniot points as possible. Figure 7 shows the variation of the pair correlation function, $g(r)$, for conditions along the Hugoniot. As the density and temperature increase along the Hugoniot, the correlation in the liquid decreases as shown by the first maximum in the pair correlation function varying from a value around 3.5 at $\rho = 6.5 \text{ g/cm}^3$ to around 3.2 at $\rho = 7.5 \text{ g/cm}^3$.

At a density of 6.5 and 7.5 g/cm³, the QMD calculations starting from an FCC structure remain in that structure for our limited simulation times. These simulations were started in a liquid state to circumvent overcoming the heat of fusion close to the melting point. It should be noted that calculating chemical properties near the melt transition is dangerous
<figure><img src="image_8.png" /><figcaption>FIG. 7. (Color online) 256-atom *ab initio* (solid) and MEAM pair correlation functions (dashed), $g(r)$, calculated at conditions surrounding the Hugoniot points.</figcaption></figure>

because determining equilibrium can be illusive. Thus, we do not predict an increase in the melting temperature as a function of density as commonly found for other simple metals. We conjecture that this behavior is more likely due to the usual overestimation of the melting temperature due to the limited number of atoms used in the simulation cell and limited simulation time.

The slow decrease in correlation along the low-pressure points of the Hugoniot is directly related to the small temperature variation taking place for densities up to 8.5 g/cm³. As the density further increases above this value, temperature significantly increases to reach up to 8430 K at the highest density. A rapid decrease in the correlation of the system is associated with this significant temperature elevation with the maximum of the pair correlation function now significantly broadened and reduced to a value close to one. A direct comparison with the correlation functions obtained using the MEAM potential shows an overall agreement with the QMD results with, however, some noticeable differences. The MEAM system appears less correlated between 7.5 and 9.2 g/cm³ with a pronounced shoulder around 3 Å. This shoulder is associated with the formation of dimers that remain up to very high pressures. We further note that to provide for the stability of the A11 structure, the MEAM parameters were selected to promote the formation of dimers. This may result in an overestimation of the dimerization in the liquid state, as the QMD results do not show any dimer formation for any of the conditions simulated along the Hugoniot.

### C. Electrical properties

Before studying the variation of the electrical properties near the Hugoniot, we first turn to the value obtained at ambient density and a temperature $T = 800$ K, where experimental measurements as well as previous *ab initio* results are available [17,33]. Figure 8(a) shows the *k*-point convergence performed on a single snapshot of 108 atoms using up to $10^3$ *k* points in the Brillouin zone as defined within the Monkhorst-

<figure><img src="image_9.png" /><figcaption>FIG. 8. (Color online) Optical conductivity for Ga at $\rho = 6.0$ g/cm³ and $T = 800$ K. (a) Convergence of the optical conductivity as a function of *k* points for a single snapshot of a 108-atom supercell. (b) Convergence of the optical conductivity as a function of *k* points for a single snapshot of a 256-atom supercell. (c) Average optical conductivity for multiple snapshots and Drude fit of conductivity.</figcaption></figure>

Pack scheme. The number of *k* points needed to converge the DC conductivity illustrates the fact that optical properties have a slower convergence rate than the thermodynamical properties with respect to *k*-point sampling. For a larger supercell of 256 atoms, shown in Fig. 8(b), the DC conductivity converges faster with the number of *k* points. For the 256-atom cell, $4^3$ *k* points was our computational upper limit because of memory limitations. Figure 8(c) shows the optical conductance averaged over five uncorrelated configurational snapshots from the MD runs. The optical conductance was calculated using $8^3$ and $4^3$ *k* points for the 108- and 256-atom cells, respectively. The conductance was then fit to a Drude form to extract the DC conductance at $\omega = 0$.

Once converged, we find that the DC conductivity slightly departs from a Drude-like form near zero frequency. However, this departure from Drude-like behavior decreases with the size of the supercell [Fig. 8(c)]. We attribute this to numerical inaccuracy and statistical error close to zero frequency. It is tempting to attribute this dip near zero frequency in $\sigma_1(\omega)$ to a small minimum in the density of states around the Fermi energy. A similar feature was previously identified for the A11 phase and was interpreted as the competing influence of the molecular and metallic nature of the solid phase [4,5]. However, upon examination of the density of states, we see no evidence to suggest this is the case. Holender *et al.* [17] suggested that the calculation performed by Gong *et al.* [5], at the $\Gamma$ point only, did not capture the complexity of the electronic structure with such a limited simulation cell. We find that the number of states around the fermi level does not change as the number of *k* points are increased. This suggests that the sensitivity of the conductivity, Eq. (6), with respect to *k*-point
<figure><img src="image_10.png" /></figure>

FIG. 9. (Color online) Variation of the DC conductivity, $\sigma_{DC}(\omega)$ near the Hugoniot for the 108-atom cell (blue diamonds), and $\sigma_{DC}(\omega)$ from Fig. 8(c) for the 256-atom cell (red circle). Experimental measurements of DC conductivity in the solid [36] and liquid [33] phase are shown as green squares. Solid and solid-like phases shown with open shapes.

sampling is due to some dispersion in the bands. We believe that the inaccurate description of the electronic structure is also the origin of the poor description of the thermodynamical properties found in the previous section. While the result found here confirms, to first order, the calculation performed by Holender et al. [17], we emphasize that the better converged calculations presented here still indicate a small deviation from a Drude form in the liquid. Given the trend with cell size we expect that if we could afford to perform this calculation with more atoms this deviation would disappear. This result confirms that this system behaves very closely to a simple metal as suggested earlier [17].

The value of the DC conductivity obtained [$\sigma_{DC}(\omega) = 2.97 \times 10^6 (\Omega m)^{-1}$] compares well with the experimental result [33], [$\sigma_{DC}(\omega) = 2.79 \times 10^6 (\Omega m)^{-1}$] near this condition. This result is also in better agreement with the experimental value than the previous *ab initio* calculations performed using 64 atoms, eight *k* points, a norm conserving potential and within the LDA approximation [17] [$2.0-2.5 \times 10^6 (\Omega m)^{-1}$]. We see in Fig. 8(a) that the calculation using only $2^3$ *k* points leads a DC conductivity in better agreement with their reported value and much closer to a Drude form, which shows that previous calculations were not fully converged in either *k*-points sampling and/or particle number.

We now turn to the behavior of the optical conductivity as both the density and temperature increase along the Hugoniot.

Due to computational limitations we use a 108-atom cell and $8^3$ *k* points. Figure 9 shows QMD calculations of $\sigma_{DC}$ near the Hugoniot in blue. Experimental conductivity measurements in liquid-Ga [33] versus temperature are shown in green with the density values extracted from the density-temperature relations in Ref. [35]. The first two Hugoniot calculations are very close to the melt transition, making it very difficult to determine a true equilibrium state via QMD. Liquid-like snapshots are shown in solid blue diamonds. Solid-like configurations are shown with open diamonds. There is a large difference in correlation between solid- and liquid-like structures, which is reflected in the calculated values of $\sigma_{DC}$. This difference is readily seen in experimental measurements of $\sigma_{DC}$ in the liquid [33] and solid [36]. Given the difficulty in reaching equilibrium via QMD close to the melt transition, experimental measurements are needed to ascertain the conductivity for temperatures below 1000 K along the Hugoniot. Above a density of 7.5 g/cm$^3$ the temperature increases rapidly along the Hugoniot. This increase in temperature corresponds to an optical conductivity recovering simple liquid-metal behavior where the DC conductivity as a function of temperature decreases.

## IV. CONCLUSION

Using *ab initio* molecular dynamics simulations, we calculate the thermodynamic and transport properties of liquid gallium up to 150 GPa. These simulations show that the MEAM potential accurately reproduces the *ab initio* dynamic properties, although it was fit only to equilibrium *ab initio* data. The discrepancies remaining with the experimental data for the diffusivity and viscosity can thus be attributed to either a shortcoming of the DFT at describing liquid gallium or a problem with the accuracy of experimental measurements. There are several factors contributing to the DC conductivity along the Hugoniot. The DC conductivity is dominated by density increasing correlation and causing an increase in conductivity at densities below 7.5 g/cm$^3$ along the Hugoniot trajectory. Above 7.5 g/cm$^3$, at a shock pressure of approximately 20 GPa, temperature rapidly increases, thus decreasing the conductivity.

## ACKNOWLEDGMENTS

We acknowledge funding for this work under the LANL Laboratory Directed Research and Development Project No. LDRD-20050107DR, “Methodologies to measure material response during dynamic loading at microscopic times and length scales.” Los Alamos National Laboratory is operated under U.S. Department of Energy Contract No. DE-AC52-06NA25396. We also acknowledge our colleagues who have provided useful comments in the preparation of this paper, especially Cindy Bolme, Jim Glownia, Brad Holian, and Ramon Ravelo.

[1] M. I. Baskes, S. P. Chen, and F. J. Cherne, Phys. Rev. B 66, 104107 (2002).

[2] O. Schulte and W. B. Holzapfel, Phys. Rev. B 55, 8122 (1997).

[3] D. A. Young, Phase Diagrams of the Elements (University of California Press, Berkeley, CA, 1991).

[4] X. G. Gong, G. L. Chiarotti, M. Parrinello, and E. Tosatti, Phys. Rev. B 43, 14277 (1991).
[5] X. G. Gong, G. L. Chiarotti, M. Parrinello, and E. Tosatti, *Europhys. Lett.* **21**, 469 (1993).

[6] M. Bernasconi, G. L. Chiarotti, and E. Tosatti, *Phys. Rev. B* **52**, 9988 (1995).

[7] S. N. Rapeanu and I. Padureanu, *Physica Scripta* **T57**, 18 (1995).

[8] L. A. Collins, S. R. Bickham, J. D. Kress, S. Mazevet, T. J. Lenosky, N. J. Troullier, and W. Windl, *Phys. Rev. B* **63**, 184110 (2001).

[9] S. Mazevet, J. D. Kress, L. A. Collins, and P. Blottiau, *Phys. Rev. B* **67**, 054201 (2003).

[10] S. Mazevet, J. Kress, and L. A. Collins, *AIP Conf. Proc.* **730**, 139 (2004).

[11] J. N. Fritz and W. J. Carter, “The Hugoniot EOS for liquid Ga and Hg”, Los Alamos National Laboratory Technical Report No. LA-13844-MS (2001).

[12] R. M. Martin, *Electronic Structure* (Cambridge University Press, Cambridge, UK, 2004).

[13] E. Schwegler, M. Sharma, F. Gygi, and G. Galli, *Proc. Natl. Acad. Sci. USA* **105**, 14779 (2008).

[14] S. Mazevet, J. Clerouin, V. Recoules, P. M. Anglade, and G. Zerah, *Phys. Rev. Lett.* **95**, 085002 (2005).

[15] M. P. Desjarlais, J. D. Kress, and L. A. Collins, *Phys. Rev. E* **66**, 025401 (2002).

[16] S. Mazevet, M. Torrent, V. Recoules, and F. Jollet, *High Energy Density Phys.* **6**, 84 (2010).

[17] J. M. Holender, and M. J. Gillan, M. C. Payne, and A. D. Simpson, *Phys. Rev. B* **52**, 967 (1995).

[18] F. J. Cherne, M. I. Baskes, T. C. Germann, R. J. Ravelo, and K. Kadau, in *Shock Compression of Condensed Matter-2003*, edited by M. D. Furnish, Y. M. Gupta, and J. W. Forbes, AIP Conf. Proc. No. 706 (American Institute of Physics, Melville, NY, 2004), p. 281.

[19] K. Kadau, F. J. Cherne, R. J. Ravelo, and T. C. Germann, *Phys. Rev. B* **88**, 144108 (2013).

[20] G. Kresse and D. Joubert, *Phys. Rev. B* **59**, 1758 (1999).

[21] J. P. Perdew and Y. Wang, *Phys. Rev. B* **45**, 13244 (1992).

[22] D. Vanderbilt, *Phys. Rev. B* **41**, 7892 (1990).

[23] G. Kresse and J. Hafner, *J. Phys. Condens. Mat.* **6**, 8245 (1994).

[24] Y. B. Zel'dovich and Y. P. Raizer, *Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena* (Academic Press, New York, NY, 1967).

[25] M. P. Allen and D. J. Tildesley, *Computer Simulation of Liquids* (Oxford University Press, Oxford, UK, 1987).

[26] E. Chisolm and D. Wallace, in *Shock Compression of Condensed Matter-2005*, edited by M. D. Furnish, M. L. Elert, T. P. Russell, and C. T. White, AIP Conf. Proc. No. 845 (American Institute of Physics, Melville, NY, 2006), p. 53.

[27] W. A. Harrison, *Solid State Theory* (McGraw-Hill, New York, 1970).

[28] G. D. Mahan, *Many Particle Physics* (Plenum Press, New York, 1981).

[29] J. Callaway, *Quantum Theory of the Solid State* (Academic Press, New York, NY, 1974).

[30] P. E. Blöchl, *Phys. Rev. B* **50**, 17953 (1994).

[31] J.-B. Maillet, M. Mareschal, L. Soulard, R. Ravelo, P. S. Lomdahl, T. C. Germann, and B. L. Holian, *Phys. Rev. E* **63**, 016121 (2000).

[32] J. H. Rose, J. R. Smith, F. Guinea, and J. Ferrante, *Phys. Rev. B* **29**, 2963 (1984).

[33] A. Ben Abdellah, J. G. Gasser, A. Makradi, B. Grosdidier, and J. Hugel, *Phys. Rev. B* **68**, 184201 (2003).

[34] H. Riedl, T. Persgon, and P-E. Eriksson, *Int. J. Appl. Radiat. Isot.* **30**, 481 (1979).

[35] H. V. Tippelskirch, *Ber. Bunsenges. Phys. Chem.* **80**, 726 (1976).

[36] Edited by D. R. Lide, *CRC Handbook of Chemistry and Physics* (CRC Press, Boca Raton, FL, 2003).