Properties of B4C in the shocked state for pressures up to 1.5 TPa

TL;DR

This study combines density functional theory and molecular dynamics to explore the behavior of boron carbide under extreme pressures up to 1.5 TPa, revealing phase transitions, segregation, and fluid states.

Contribution

It provides the first comprehensive analysis of B4C's properties under ultra-high pressures, including phase transitions and fluid states, using advanced computational methods.

Findings

Discontinuities in the Hugoniot at 50 GPa and 90-100 GPa indicating phase changes.

Good agreement between theoretical predictions and experimental data up to 700 GPa.

Identification of amorphous and atomistic fluid states in B4C at high pressures.

Abstract

Density Functional Theory calculations using the quasi-harmonic approximation have been used to calculate the solid Hugoniot of two polytypes of boron carbide up to 100 GPa. Under the assumption that segregation into the elemental phases occurs around the pressure that the B11Cp(CBC) polytype becomes thermodynamically unstable with respect to boron and carbon, two discontinuities in the Hugoniot, one at 50 GPa and one at 90-100 GPa, are predicted. The former is a result of phase segregation, and the latter a phase transition within boron. First principles molecular dynamics (FPMD) simulations were employed to calculate the liquid Hugoniot of B4C up to 1.5 TPa, and the results are compared to recent experiments carried out at the Omega Laser Facility up to 700 GPa [Phys. Rev. B 94, 184107 (2016)]. A generally good agreement between theory and experiment was observed. Analysis of the FPMD simulations provides evidence for an amorphous, but covalently bound, fluid below 438 GPa, and an atomistic fluid at higher pressures and temperatures.

Tables

Table 1. Table 1: Initial starting conditions for the B 11 C p p {}_{\text{p}} (CBC), B 12 (CCC), and B/C solid state Hugoniots

Structure	${\rho }_0$ (g/cc)	$T_0$ (K)	$P_0$ (GPa)	$E_0$ (kJ/g)
CBC${}_{\text{p}}$	2.529	300	0	-61.6922
CCC	2.471	300	0	-61.0874
B/C (graphite)	2.353	300	0	-60.8083

Equations

$$\gamma =V{\left(\frac{\partial P}{\partial E}\right)}_V$$

$$C_s={\left(\frac{\partial P}{\partial \rho }\right)}_S^{1/2},$$

$$C_s={\left(\frac{\partial P}{\partial \rho }\right)}_T+T\frac{{\left[\frac{{\left(\frac{\partial P}{\partial T}\right)}_{\rho }}{\rho }\right]}^2}{{\left(\frac{\partial E}{\partial T}\right)}_{\rho }}.$$

$$g(r)=4\pi r^2\rho dr.$$

$$\beta (t)=\frac{⟨\vec{B}(0)\cdot \vec{B}(t)⟩}{⟨\vec{B}(0)\cdot \vec{B}(0)⟩},$$

$$E-E_0=\frac{1}{2}(P+P_0)\left(\frac{1}{{\rho }_0}-\frac{1}{\rho }\right),$$

$$T_{\text{Lindemann}}=f^2\frac{\bar{m}{k}_B{\bar{\nu }}^{\frac{2}{3}}{\theta }^2}{9\hslash },$$

Full text

I Introduction

Because of its outstanding hardness, thermodynamic stability, high Hugoniot elastic limit (HEL), low density, and electronic properties boron carbide has many applications, e.g. as a refractory material, in abrasive powders and ballistics, as a neutron radiation absorbent, and in electronics Domnich et al. (2011). However, boron carbide does not behave according to theoretical predictions based on the thermodynamic limit Shirai et al. (2014a); Yao et al. (2015). For example, it displays an abrupt drop in shear strength past $\sim$20 GPa, and it experiences brittle failure under impact Vogler et al. (2004); Dandekar (2001), which complicates the use of this material in real world applications.

Boron carbide is composed of 12 atom icosahedra cross-linked by three-atom chains. The broad composition range (from 8-20 atom % carbon) Thevenot (1990); Schwetz and Karduck (1991); Gosset and Colin (1991); Kisly et al. (1988), and presence of vacancies makes this a very complex material. In addition, the constituents of a given sample may span a wide composition range and number of defects Werheit (2016). It has been proposed that defects are essential for maximizing the bonding in boron carbide Balakrishnarajan et al. (2007), and they may control the electronic properties Hushur et al. (2016). Because of boron’s propensity to form icosahedral units, at first one may suppose that in the idealized carbon rich B4C stoichiometry all of the boron atoms would be found in the icosahedra, and the carbon atoms in the chains, i.e. B12(CCC). Electronic structure calculations, primarily based upon density functional theory (DFT) Saal et al. (2007); Bylander et al. (1990); Armstrong et al. (1983); Guo et al. (2006); Aydin and Simsek (2009); Vast et al. (2000); PhysRevB.90.064109; Ektarawong et al. (2015, 2016), however, have shown that this is not the most stable atomic arrangement: one of the boron atoms prefers to be located in the center of the three atom chains. There are two inequivalent positions in each icosahedron: those linking the icosahedral cages (commonly referred to as the ‘polar’ sites), and the ones bonded to the three atom chains (the ‘equatorial’ sites). Either a polar or an equatorial boron can be substituted by a carbon atom, and the former configuration, B11C${}_{\text{p}}$(CBC), is energetically preferred.

Under static compression, the icosahedral units are less compressible than the space between them Dera et al. (2014). Raman spectra up to 50 GPa Yan et al. (2009); Guo et al. (2010) and X–ray diffraction measurements up to 126 GPa Fujii et al. (2010) have not observed any evidence indicative of a structural phase transition. Under dynamic compression (until recently Fratanduono et al. (2016) equation of states (EOS) data from shock compression have been available up to pressures of 200 GPa Pavlovskii (1971); Zhang et al. (2006); Gust and Royce (1971); Vogler et al. (2004); Dandekar (2001); Holmquist and Johnson (2006); Ciezak and Dandekar (2009); Grady (1994, 2010); McQueen et al. (1970); Marsh (1980); Grady and Wise (1993)) the behavior is more complex. Anomalies in the shock Hugoniot that are found between 20-60 GPa are suggestive of a phase transition Grady (2015). Pressure–induced amorphization has been observed in shock–loaded samples Chen et al. (2003), as well as in indented Domnich et al. (2002), and mechanically scratched Chen and McCauley (2006) samples. Nanoindentation experiments have found narrow amorphous bands and areas with local disorder as well as evidence for an amorphous phase above 40 GPa Ge et al. (2004); Yan et al. (2006); Reddy et al. (2013). Numerous proposals have been put forward to explain the mechanism of amorphization and the gradual loss of strength that boron carbide experiences under pressure Yan et al. (2006, 2009); Betranhandy et al. (2012); Fanchini et al. (2006); An et al. (2014); An and Goddard III (2015); Reddy et al. (2013); Korotaev et al. (2016a, b); Raucoules et al. (2011); Taylor (2015); Wilkins (1968); Holmquist and Vogler (2009). These include: the bending of the CBC chains, the formation of new B–B bonds that are kinetically trapped, the presence of defects in the CBC chain, as well as the disassembly of the octahedra under shear stress.

The shock Hugoniot of boron carbide has recently been measured up to 700 GPa, and a plateau in the shock velocity was observed that suggests a further phase change above $\sim$130 GPa Fratanduono et al. (2016). However, until now first principles calculations of the shock Hugoniot were available only up to 80 GPa. Using first principles molecular dynamics (FPMD) simulations Taylor showed that the shock properties of the solid phase of boron carbide are polytype dependent Taylor (2015). A discontinuity in the Hugoniot that could be explained by the bending of the three atom chain was observed for B12(CBC), B11C${}_{\text{p}}$(CCC), B11C${}_{\text{p}}$(CCB) at 20-50 GPa, but not for B11C${}_{\text{p}}$(CBC) up to 80 GPa. However, the introduction of a vacancy in the center of the chain in the latter structure resulted in a discontinuity at 16 GPa. It was noted that this is remarkably close to the pressure at which the HEL is observed, $\sim$20 GPa.

Herein, we report the solid state Hugoniot and properties of the B11C${}_{\text{p}}$(CBC) and B12(CCC) polytypes of boron carbide obtained using the quasi-harmonic approximation. In addition, FPMD simulations have been carried out to obtain the EOS of B11C${}_{\text{p}}$(CBC) up to 1.5 TPa and 60,000 K. The results are in good agreement with Omega laser measurements Fratanduono et al. (2016), and they are relevant for ongoing and future inertial confinement fusion experimental efforts at high energy laser facilities such as the National Ignition Facility at Lawrence Livermore National LaboratoryMoore et al. (2016). Due to the technical complexity of such experiments, their interpretation has proven to be challenging Hicks et al. (2009); Knudson et al. (2004, 2015). In the past, theoretical calculations have played an important role in explaining and validating shock-wave measurements Boates et al. (2011); Zhang et al. (2012a, b, 2011); Root et al. (2013); Magyar et al. (2014); Mattsson et al. (2014); Li et al. (2013), and here they provide crucial insight into the behavior of the B4C fluid. An analysis of the FPMD simulations suggests that below 438 GPa (7,683 K) boron carbide behaves like a covalently bonded fluid, but higher pressures and temperatures lead to a liquid-liquid phase transition to an atomistic fluid.

II Computational Details

II.1 Solid State Hugoniot

Benchmark calculations were carried out with the LDA Perdew and Zunger (1981), PBE Perdew et al. (1996) and Am05 Armiento and Mattsson (2005); Mattsson and Armiento (2009) functionals. The thermal expansion was calculated under the quasi-harmonic approximation to obtain the density of the B11C${}_{\text{p}}$(CBC) structure at 300 K. The simulations were conducted with 3- and 4-electron projector augmented wave (PAW) method potentials Blöchl (1994) with 1.7 and 1.5 Bohr core radii for boron and carbon atoms, respectively, and a 600 eV plane-wave cutoff. Comparison of the results with the experimental density of boron carbide revealed that the PBE functional performs best for this system (details are provided in the supplemental information (SI), Fig. S1). Therefore, phonon calculations were carried out for B11C${}_{\text{p}}$(CBC), B12(CCC), boron, and carbon using the Perdew-Burke-Ernzerhof generalized gradient approximation (PBE-GGA) as implemented in VASP Kresse and Hafner (1993, 1994) combined with the phonopy Togo and Tanaka (2015) package under the harmonic approximation in increments of 10 GPa from 0 to 100 GPa.

The B11C${}_{\text{p}}$(CBC) is the most stable structure between 0-50 GPa. Above 50 GPa the segregation of B11C${}_{\text{p}}$(CBC) into elemental carbon and boron becomes enthalpically preferred. In the pressure range examined boron undergoes two structural phase transitions ($\alpha\rightarrow\gamma\rightarrow\alpha-\text{Ga}$) Ogitsu et al. (2013); Sanz et al. (2002); Ogitsu et al. (2009); Oganov et al. (2009); Albert and Hillebrecht (2009); Parakhonskiy et al. (2011); Häussermann et al. (2003) and carbon one (graphite $\rightarrow$ diamond) Correa et al. (2006). The pressures we calculate for these phase transitions, which are provided in the SI (Figs. S10-S11) match well with previous studies. The supercells for the phonon calculations were chosen such that the number of atoms in the simulation cell was always greater than or equal to 96 atoms, typically being $2\times 2\times 2$, or $3\times 3\times 3$ representations of the standard primitive cells. Carbon in the graphite phase at ambient pressure required a $4\times 4\times 2$ representation. Each structure was initially relaxed using a 600 eV plane-wave cutoff, and an automatically generated $\Gamma$-centered Monkhorst-Pack grid where the number of divisions along each reciprocal lattice vector was chosen such that the product of this number with the real lattice constant was 40-45 Å. For the supercell calculations this was reduced to 35 Å. The vibrational internal energy was calculated via phonopy from 0 to 5000 K in increments of 10 K to determine at which temperature the Hugoniot equation was satisfied for each pressure. In addition, the Debye frequency, $\theta_{D}$, was calculated at each pressure for use in the Lindemann equation to estimate the melting curve.

II.2 Liquid Hugoniot and Properties

First-principles molecular dynamics (FPMD) simulations on B11C${}_{\text{p}}$(CBC) have been carried out for pressures up to 1.5 TPa using DFT Kohn and Sham (1965) within the Perdew-Burke-Ernzerhof generalized gradient approximation (PBE-GGA) Perdew et al. (1996) as implemented in VASP. In order to estimate the internal energy of the system at high temperature, the Mermin functional was used Mermin (1965). The simulations were carried out in the canonical (NVT) ensemble using Born-Oppenheimer dynamics with a Nosé-Hoover thermostat. For each pressure ($P$) and temperature ($T$) run, the system was equilibrated within 1-2 ps and simulated using a 0.75 fs ionic time-step. All FPMD simulations were carried out with 120-atom supercells, 3- and 4-electron PAW potentials Blöchl (1994) with 1.7 and 1.5 Bohr core radii for boron and carbon atoms, respectively, and a 700 eV plane-wave cutoff. For the $P$-$T$ points closest to the Hugoniot the FPMD simulations were repeated with harder potentials: 1.1 Bohr core radii and a 900 eV plane-wave cutoff. The resulting corrections to the pressure, energy, and structural parameters are minimal (exact $P$ and $E$ values are given in the SI (Table S3)).

It has been shown that for calculations carried out on fcc carbon using the aforementioned “hard” PAW potentials, and 4 valence electrons the largest deviation in the EOS was $<$0.1% up to 3.5$\times$105 GPa as compared with all-electron calculations using the full-potential linear augmented plane wave (LAPW), and the augmented plane wave plus local orbital (APW+lo) methods Benedict et al. (2014). In addition, the same study found that 4 valence electrons were sufficient for FPMD calculations carried out for $T<$100,000 K. The $P/T$ conditions we consider are lower, further supporting the choice of the computational protocol.

For temperature and density points closest to the Hugoniot, the Grüneisen parameter, sound speed at the shock state, radial distribution functions, and bond auto-correlation functions (BACF) were evaluated from the FPMD trajectories. The Grüneisen parameter, $\gamma$, is a dimensionless parameter describing the effect that changing the volume of a solid has on its vibrational properties. The Grüneisen equation was used to estimate this parameter,

[TABLE]

where $V$ is the volume and $E$ is the internal energy. The calculated values can be compared to those obtained experimentally Fratanduono et al. (2016). The sound speed at the shock state, $C_{s}$, is measured by using the velocity interferometer system for any reflector (VISAR) technique. The sound speed can be calculated using the following equation:

[TABLE]

where $\rho$ is the density and $S$ is the entropy. The sound speed may also be obtained using the approximate equation:

[TABLE]

In studies of the liquid state it is useful to analyze the radial distribution function, $g(r)$, for atom pairs given as

[TABLE]

The bond auto-correlation function (BACF) is defined as

[TABLE]

where the bond vector, $\vec{B}$, is a byte vector with a value of 1 for each bond of a certain type within a cutoff distance, and 0 otherwise. The correlation function describes the average duration of a bond between two atoms. The decay timescale (bond lifetime) is not strongly dependent on the bond cutoff employed: we used 1.5 Å for C-C, 1.6 Å for B-C, and 1.7 Å for B-B bond cutoff distances.

III Results and Discussion

III.1 Solid Hugoniot

To produce the solid state Hugoniot equation of state (EOS) for B11C${}_{\text{p}}$(CBC), B12(CCC), and stoichiometric combinations of segregated boron and carbon (B/C) a series of phonon calculations were conducted under the quasi-harmonic approximation from 0 to 100 GPa. In the solid the total energy, $E$, is given as a sum of the electronic energy, $E_{\text{elec}}$, at 0 K coupled with the temperature dependent vibrational energy, $E_{\text{vib}}$. The $E_{\text{vib}}$ and corresponding temperatures were chosen so that at each pressure the Hugoniot relation,

[TABLE]

was satisfied. The Hugoniot relation shows how the energy, $E_{0}$, density, $\rho_{0}$, and pressure, $P_{0}$, of the initial state are directly related to the final conditions of the shocked state. The initial starting conditions used for each Hugoniot curve are provided in Table 1.

The solid state data was used to construct the five different Hugoniot equations of states shown in Fig. 1. Two of these corresponded to the Hugoniots of the pure B4C polytypes (referred to as CBCp and CCC hereafter). One corresponds to phase segregated boron and carbon using the graphite phase at ambient pressure. Finally, two plots show how the Hugoniot would behave if the polytypes of B4C would undergo a direct phase transition to segregated boron and carbon at 50 GPa, the pressure at which CBCp is predicted to no longer be enthalpically stable using 0 K DFT calculations.

The solid Hugoniots of the two B4C polytypes do not exhibit any discontinuities, which is to be expected since they do not undergo any structural phase transitions in the pressure range considered. The Hugoniots were not determined past $\sim$100 GPa because of the emergence of imaginary frequencies in the phonon calculations by 120 GPa for CBCp. Recently, the solid Hugoniot of several boron carbide polytypes have been obtained via first-principles molecular dynamics simulations up to 80 GPa Taylor (2015). Of the polytypes considered in those calculations only the CBCp structure was also examined herein. Despite the quasi-harmonic method being less robust than molecular dynamics at high temperatures, the plots of shock versus particle velocity obtained using the two methods are in good agreement with each other.

Calculations were also carried out to determine the Hugoniot of a segregated mixture of B/C, with the B4C stoichiometry. As Fig. 1(A) shows, two abrupt changes in the slope are observed in the pressure versus density plot at $\sim$20 and $\sim$100 GPa. The pressures at which these discontinuities occur correspond to the pressures at which boron is calculated to transition from the $\alpha$ to the $\gamma$ phase, as well as to the $\alpha$-Ga phase, respectively, at 0 K.

Our calculations show that even though the CBCp and CCC structures are dynamically stable to at least 100 GPa, they become thermodynamically unstable with respect to decomposition into elemental boron and carbon above $\sim$50 GPa and $\sim$20 GPa, respectively. This is in reasonable agreement with the results of Fanchini and co-workers, whose DFT calculations with Gaussian basis sets showed that all of the polytypes of boron carbide become metastable past 6-7.5 GPa Fanchini et al. (2006). Fanchini et. al. also estimated the pressures at which the activation barriers to segregation could be overcome and found that whereas B11C(CBC) with the icosahedral carbon either in the polar or equatorial position could withstand stresses of up to 40 GPa, the CCC polytype had a low barrier for segregation and would transition into layers of B12 and graphite near 6 GPa. These results were used to explain the experimental observation of abrupt failures related to the drop in shear strength via shock induced deformation, and the presence of amorphous bands in shock loaded, indented, and scratched experiments Grady (2015); Chen et al. (2003); Chen and McCauley (2006); Domnich et al. (2002); Ge et al. (2004); Yan et al. (2006); Reddy et al. (2013). It should be noted that other computational studies have led to different proposed mechanisms including the fracture of the icosahedra An and Goddard III (2015), and the formation of new bonds upon depressurization Yan et al. (2009) or upon shear deformation An et al. (2014), to explain these observations.

On the other hand, experiments in diamond anvil cells have not observed any phase transitions when boron carbide was compressed up to over 100 GPa Yan et al. (2009); Guo et al. (2010); Fujii et al. (2010). Our phonon calculations indicate that the CBCp and CCC phases are dynamically stable to at least this pressure. This behavior, similar to that of elemental boron Sanz et al. (2002), suggests the presence of a large kinetic barrier to elemental segregation. However, the quickly increasing temperature in shock experiments may allow the polytypes of boron carbide to overcome the kinetic barriers preventing phase segregation at the pressures where they become thermodynamically unstable, i.e. 20-50 GPa. Phase separated boron carbide, or simply boron and carbon, could then undergo phase transitions similar to those found in their pure states.

Therefore, we have also constructed the Hugoniot of the CBCp and CCC polytypes assuming that they both phase segregate into elemental boron and carbon at 50 GPa. As Fig. 1(A) and (C) reveal, this scenario leads to discontinuities in the Hugoniot at the pressure at which segregation is assumed to occur, and the pressure at which boron transitions to the $\alpha$-Ga phase (90-100 GPa). The slight difference between these plots and the one calculated for segregated boron and carbon past 50 GPa stem from the difference in their initial densities (see Table 1).

Recent calculations of the Gibbs free energies predict that at 0 GPa two configurational phase transitions will occur within B4C at around 870 K and 2325 K from an ordered phase with $Cm$ symmetry to disordered phases with $R3m$ and $R\bar{3}m$ symmetry Ektarawong et al. (2014). Monte Carlo simulations were further employed to study these temperature induced phase transitions Yao et al. (2015). Another study found that the B4C phase with monoclinic symmetry becomes thermodynamically equivalent to a rhombohedral B13C2 phase at high temperatures Widom and Huhn (2012). In addition, it was shown that B4C will transform to B2.5C plus $\gamma$-boron in a pressure range of 40-60 GPa Ektarawong et al. (2016). Whereas an ordered form of B2.5C is favored below 500 K, a disordered form is preferred at higher temperatures. At higher pressures $\gamma$-boron and diamond are favored, regardless of the temperature. Only the B2.5C stoichiometry was found to be thermodynamically stable under high pressure with respect to B4C and the elemental phases. These predicted temperature induced transformations, which are not considered herein, are also likely to lead to discontinuities in the solid state Hugoniot. On the other hand, it may be that under shock compression boron carbide initially transitions to an amorphous form.

Incomplete phase segregation and the occurrence of only a subset of these transitions in some samples, or some regions of a given sample, coupled with the presence of defects and vacancies may explain the discrepancies found in the experimental shock data. The low temperatures frequently employed in diamond anvil cells prevents the material from overcoming the kinetic barriers associated with amorphization and segregation, and we suggest that this is the reason why phase changes have not been observed in these experiments.

The Lindemann criterion was used to predict the melting temperature of the two B4C polytypes using the equation:

[TABLE]

where $\theta$ is the Debye temperature, $\bar{m}$ and $\bar{\nu}$ are the mean atomic mass and volume, $k_{B}$ is the Boltzmann constant, $\hbar$ is the reduced Planck’s constant, and $f$ is a unitless parameter described as the critical ratio of vibrational amplitude to atomic spacing at melting. The Debye frequencies were calculated every 10 GPa, see Fig. 2(A), and converted into Debye temperatures. An experimental melting point for boron carbide at ambient pressure, 2720 K Mukhanov and Solozhenko (2015), was employed to estimate $f$. This yielded values of 0.143 and 0.138 for the CBCp and CCC structures, respectively. It was assumed that $f$ remained constant in the pressure range considered.

Interestingly, both B4C polytypes display a negative slope of melting, Fig. 2(B), with respect to increasing pressure up to at least 50 GPa. This is a relatively rare phenomenon that is found, for example, in compressed lithium Guillaume et al. (2011), and sodium Gregoryanz et al. (2005). The slope of the melting curve in the pure boron phase diagram, on the other hand, is always positive in the same pressure range Oganov et al. (2009). For the CCC polytype $T_{\text{Lindemann}}$ decreases from 2720 K to just under $\sim$2400 K from ambient pressure up to 100 GPa. For the CBCp structure the melting temperature decreases only until 50 GPa where it is predicted to be $\sim$2200 K. At this pressure the slope changes sign and $T_{\text{Lindemann}}$ increases to $\sim$2400 K by 100 GPa. The gradual decrease in melting temperature in the two B4C polytypes is less dramatic than that observed in sodium, which exhibits a melting point decrease from $\sim$1000 K at 31 GPa to $\sim$300 K at 118 GPa Gregoryanz et al. (2005). Our results suggest that the differences in calculated melting curves for various structural candidates could potentially be used to help characterize a given sample of boron carbide. However, more accurate techniques to estimate melting, such as two-phase simulations, would be required to confirm these preliminary results. For example, incongruent melting could occur in B4C, as in alloys.

III.2 Liquid Properties and Structure

The specific EOS points that were used to interpolate the Hugoniot of liquid B4C with the initial starting condition of $\rho_{0}=$ 2.529 g/cc, T${}_{0}=$ 300 K, and P0 at ambient pressure are provided in the SI (Table S4). For each of the selected densities, the data can be readily used to accurately determine the pressures and the temperatures that satisfy the Hugoniot relation, Eq. 6. For this purpose, we used a linear interpolation of Eq. 6 on a grid of temperatures. In Fig. 3 the results are directly compared with experimental measurements obtained by laser compression at the Omega laser facility Fratanduono et al. (2016). In an experimental shock experiment the VISAR technique is capable of measuring the shock and particle velocities, as well as the sound speed. Because it is easier to accurately determine pressures as opposed to densities in experiment, much of the data below is discussed in terms of the pressure. However, we note that Fig. 3(A) can be employed to determine the density that a given pressure corresponds to.

In Fig. 4 the temperature, density, and pressure obtained using the FPMD calculations along the Hugoniot have been converted to sound speeds, to directly compare with the available experimental data. At lower shock and particle velocities the agreement between theory and experiment is excellent. However, the experimental Hugoniot is softer at higher pressures. This disagreement could be explained by errors associated with the quartz reference used in the impedance matching analysis, systematic errors in the experimental measurements, or uncertainties in the metrology of the B4C sample. Studies to support the differences in properties based on the experimental procedure are still ongoing and inconclusive at this point. Despite the slight deviation in the shock and particle velocities the agreement between the FPMD and Omega laser results is good, so that the FPMD calculations provide a starting point for comparison with higher energy laser experiments.

The Grüneisen parameter, $\Gamma$, is one of the few experimental quantities that can be extracted during a laser shock experiment. Using the results of our FPMD calculations we have obtained this parameter on the principle Hugoniot via Eq. 1, as shown in Fig. 5. The maximum in $\Gamma$ occurs at $\sim$500 GPa (5.14 g/cc, 10,242 K). For every density considered $\Gamma$ remains below 1. It is almost twice as large as the value extracted from experiment, and falls just outside of the experimental error bars. $\Gamma$ is related to the slope of the Hugoniot and the sound speed at the shock state, $C_{s}$. Since the experimental and theoretical $C_{s}$ values agree well, but the slopes of the Hugoniot differ, the theoretically determined $\Gamma$ values will also be different.

Fig. 6 A shows the radial distribution function, $g(r)$, of the liquid at 333 GPa (4.73 g/cc, 3,898 K). The $g(r)$ calculated for a wide range of densities can be found in the SI (Figs. S12-S28). A crystal can possess short (atom-atom), mid (molecular units), and long range (crystalline) order. As a crystal transitions into the liquid phase the long range order disappears. At pressures lower than 438 GPa (5.00 g/cc, 7,683 K) the radial distribution functions show clear mid-range order with crisp first and second coordination shells for all of the atom-atom pairs as a result of the strong covalent interactions. According to Hume-Rothery theory a change in the atomic composition can lead to the electronic density of states exhibiting pseudo-gap behavior where the Fermi level falls on a minimum in the density of states. Similarly, the presence of a pseudo-gap caused by a related Peierls-like mechanism in the short and mid-range atomic ordering has been observed for compressed liquid sodium Raty et al. (2007). As the density of states figures in the SI (Figs. S49-S65) illustrate, the pseudo-gap gradually closes at higher densities as the mid-range order in boron carbide disappears. As the pressure and temperature increases the second coordination shells in the pair correlation functions become less prominent, the liquid becomes more amorphous, and mid-range order becomes less prevalent. This is a direct consequence of the disappearance of the molecular motifs, i.e. icosahedra and three atom chains, within the liquid.

In addition, despite having quite distinct radial distribution functions at lower densities, as the density of the liquid increases the $g(r)$ of the individual atom-atom pairs start to become more uniform. At 1200 GPa (6.16 g/cc, 41,169 K) there is a quick increase in all of the radial distribution functions between 1.0 and 1.5 Å, i.e. near the van der Waals radii of the boron and carbon atoms, followed by a rapid leveling off to a constant value, unlike what is observed at lower densities. This suggests that the boron and carbon atoms start to behave very similarly at this density, and that the atoms are moving around quickly enough to prevent any real coordination shells from forming for any significant length of time.

The first coordination number for each type of atom-atom pair, which was obtained by integrating the first $g(r)$ peak, is plotted as a function of density in Fig. 6 B. As would be expected, as the pressure increases the coordination numbers of the first coordination shells also generally increase, though not significantly. For example, in going from ambient conditions to 1200 GPa (6.16 g/cc, 41,169 K) the average coordination numbers change from 7.00 to 8.85, 0.66 to 1.44 and 1.42 to 2.13 for B-B, C-C and B-C pairs, respectively.

The BACF, defined in Eq. 5, are provided in Fig. 6(C) at 333 GPa (4.73 g/cc, 3,898 K) for the B-B, B-C and C-C atom pairs. The fluid is highly reactive, and any species present in the liquid at specific MD time steps would have very short lifetimes. Indeed, we see that by 1 picosecond nearly all of the atom-atom correlations have disappeared. Because C-C bonds are stronger than B-C and B-B bonds, the C-C correlation times are slightly longer than those of the other two atom pair types. As the remaining BACFs, which are provided in the SI (Figs. S29-S45), show, the lack of correlation between atom pairs becomes even more prominent at higher pressures and temperatures where the BACFs reduce to zero nearly instantly and the correlations for all of the atom pairs behave similarly. The loss of correlation with increasing density can be explained based on the comparison of the atomic interaction energy versus the kinetic energy of the system. If the overall interaction energy for the liquid system remains roughly constant, but the kinetic energy of the liquid increases as the density increases, then the system’s behavior becomes consistent with the hard sphere model Hansen and McDonald (1986); Kirkwood et al. (1950). In this model, as the kinetic energy, or temperature, increases at a constant density the attractive part of the interatomic potentials becomes irrelevant leaving only the short range repulsive portion to describe atomic ordering, resulting in a loss of atom-atom correlations at higher temperatures. This hard sphere behavior is observed in both the $g(r)$ and BACF plots of boron carbide along the Hugoniot. Moreover, higher temperatures result in an increase of the entropy, which will eventually become significantly larger than the energy associated with a covalent bond. When this occurs, the system assumes ideal gas behavior, which is consistent with the hard sphere model.

Phonon calculations showed that at ambient pressure the highest vibrational frequency, which corresponded to the CBC stretching mode, was $\sim$1600 cm*-1* in the CBCp polytype. A single vibration with this frequency would occur in 2.1 femptoseconds. Using this time scale, and the BACFs, the number of vibrations that would occur for each atom pair type before it breaks can be calculated. This was done along the Hugoniot for BACF values of 66%, 50%, 33%, and 20%. Fig. 6D shows the number of vibrations that would occur using the time associated with breaking 50% of the bonds between each atom pair type as a function of temperature. At first the number of vibrations decreases exponentially from 333 GPa (3898 K) to 438 GPa (7683 K). At higher temperatures and densities the number of vibrations for every atom-atom pair type is nearly constant, decaying very gradually. As shown in the SI (Figs. S46-S48), this trend is the same for all of the BACF values employed.

This behavior suggests that liquid B4C undergoes a liquid-liquid transition near 438 GPa (5.00 g/cc, 7,683 K). Initially, below this density, the liquid behaves like a covalently bonded fluid where ordering between atoms remains intact over long periods of time, and atom pairs can undergo a significant number of vibrations. Past 438 GPa (5.00 g/cc, 7,683 K) the liquid becomes atomistic as the time scale of the interaction between atom pairs decreases. This is expected in fluids at very high temperatures where the diffusion of the atoms through the liquid is high. Interestingly, this density correlates with the peak in the Grüneisen parameter observed in Fig. 5, which, along with the $g(r)$ plots, provides further evidence for the liquid-liquid phase transition.

Diffusivity constants, provided in the SI (Figs. S66-S68), confirmed that the boron carbide is in fact in the liquid state. It was found that as the densities and temperatures increase, the rate of diffusion of the atoms also increases preventing any long term correlation.

IV Conclusions

The solid state Hugoniot equation of states of two polytypes of boron carbide, B11C${}_{\text{p}}$(CBC) and B12(CCC), were computed to 100 GPa using the quasi-harmonic approximation. In addition, Hugoniot plots assuming segregation into stoichiometric ratios of boron and carbon at the pressure that B11C${}_{\text{p}}$(CBC) is found to become thermodynamically unstable with respect to the elemental phases, 50 GPa, were constructed. The latter displayed discontinuities at 50 GPa and 90-100 GPa, the pressure at which boron transitions to the $\alpha$-Ga phase. The results for B11C${}_{\text{p}}$(CBC) were in good agreement with those computed via first-principles molecular dynamics (FPMD) calculations up to 80 GPa Taylor (2015).

Further, we employed first principles molecular dynamics (FPMD) simulations to compute the Hugoniot of B11C${}_{\text{p}}$(CBC) up to 1.5 TPa (53,906 K). The data was compared with Omega laser measurements obtained up to 700 GPa Fratanduono et al. (2016). At lower shock and particle velocities the agreement between experiment and theory is excellent. At higher pressures the experimental Hugoniot is somewhat softer than the theoretically calculated one; the reason for this is still unclear. Excellent agreement between theory and experiment is also attained for the sound speed at the shock state, as a function of pressure.

The radial distribution functions and bond auto-correlation functions obtained from our FPMD simulations suggest the presence of clear mid-range order and strong covalent B-B, C-C and C-B interactions below 438 GPa (7,683 K). At higher pressures and temperatures the mid-range order becomes less prevalent and the liquid becomes more amorphous. Around 1200 GPa (41,169 K) the boron and carbon atoms move around quickly enough to prevent the formation of distinct coordination shells. The loss of correlation between atom pairs at higher pressures is suggestive of a liquid-liquid transition into a hard-sphere atomistic liquid.

V Acknowledgments

A.S. acknowledges financial support from the Department of Energy National Nuclear Security Administration under Award Number DE-NA0002006. The Center for Computational Research (CCR) at the University of Buffalo is acknowledged for computational support.