Experimental and Computational Studies on Superhard Material Rhenium Diboride under Ultrahigh Pressures.

An emerging class of superhard materials for extreme environment applications are compounds formed by heavy transition metals with light elements. In this work, ultrahigh pressure experiments on transition metal rhenium diboride (ReB2) were carried out in a diamond anvil cell under isothermal and non-hydrostatic compression. Two independent high-pressure experiments were carried out on ReB2 for the first time up to a pressure of 241 GPa (volume compression V/V0 = 0.731 ± 0.004), with platinum as an internal pressure standard in X-ray diffraction studies. The hexagonal phase of ReB2 was stable under highest pressure, and the anisotropy between the a-axis and c-axis compression increases with pressure to 241 GPa. The measured equation of state (EOS) above the yield stress of ReB2 is well represented by the bulk modulus K0 = 364 GPa and its first pressure derivative K0´ = 3.53. Corresponding density-functional-theory (DFT) simulations of the EOS and elastic constants agreed well with the experimental data. DFT results indicated that ReB2 becomes more ductile with enhanced tendency towards metallic bonding under compression. The DFT results also showed strong crystal anisotropy up to the maximum pressure under study. The pressure-enhanced electron density distribution along the Re and B bond direction renders the material highly incompressible along the c-axis. Our study helps to establish the fundamental basis for anisotropic compression of ReB2 under ultrahigh pressures.


Introduction
Transition metal borides have shown intriguing mechanical and structural properties combining the attractive features of metallic bonding with rigid covalent boron-boron bonding [1][2][3]. In moving across the periodic table from a group IV transition metal boride like TiB 2 to a group VI transition metal boride like ReB 2 , the boron layer transitions from a planar hexagonal net to a more puckered structure. In particular, rhenium diboride (ReB 2 ) has shown desirable mechanical properties with a high average hardness of 30-60 GPa [4][5][6][7] and bulk modulus of 334-360 GPa [4,5], comparable to that of diamonds (442 GPa) [8]. Such materials are useful for their applicability under extreme conditions requiring a combination of high-temperature chemical stability and resistance to plastic deformation. Many superhard materials (hardness above 40 GPa) such as diamonds are prone to oxidation in high-temperature environments and have a propensity for chemical reactivity with transition metals. ReB 2 shows promise as an alternative to diamonds for mechanical uses due to strong covalent bonding between B-B and Re-B atoms and high electron density [4], the compound's stability up to 2000 K, and the ease of machining by electric discharge [9]. In this study, we investigate hexagonal ReB 2  Figure 2 shows the integrated XRD powder data taken at the maximum pressure of 241 GPa with pressure determined using the platinum EOS [12]. The difference curve shown below the powder pattern in Figure 2 resulted from a fit to the hexagonal structure to ReB2. The hexagonal phase of ReB2 was found to be stable to the maximum pressure of 241 GPa. The measured lattice parameters at 241 GPa were a = 2.586 ± 0.004 Å and c = 6.882 ± 0.007 Å. Platinum peaks in Figure 2 are labeled with asterisks (*) and indexed to a face-centered cubic lattice. The platinum lattice parameter at maximum pressure of 241 GPa was measured to be a = 3.490 ± 0.009 Å.  Figure 1b, the X-ray beam was incident along the axis of compression, and scattered X-rays off the sample were captured on a Pilatus 1M detector with X-ray beam size 3.7 µm (vertical) × 3.8 µm (horizontal) FWHM (full width at half maximum) and sample-to-detector distance of 344.63 mm calibrated using the CeO 2 diffraction profile in the Dioptas software. For more information on the optical components of the Beamline 16-BM-D, refer to Park et al. [10]. Structure refinements of lattice parameters were carried out using the GSAS-II software package [11]. The measured pressure-volume data for the sample were fitted to the 3rd order Birch-Murnaghan equation of state (EOS):

Results
Here, V is the measured volume at high pressure and V 0 is the ambient pressure volume with x = V 0 /V; K 0 and K 0 ' are the bulk modulus and its first pressure derivative, respectively. To determine the initial volume V 0 of the ReB 2 sample, ambient pressure XRD measurements were separately recorded of the starting material, and the lattice parameters were determined to be a 0 = 2.901 Å and c 0 = 7.482 Å. The platinum EOS used was calibrated up to 550 GPa from Yokoo et al. [12] using the 3rd order Birch-Murnaghan EOS and employed as a pressure marker using K 0 = 276.4 GPa and K 0 ' = 5.12 with the platinum lattice parameter a = 3.924 Å at ambient pressure.
First-principles calculations are based on density functional theory (DFT) [13], which dictates that the ground state energy (or potential) of interacting electrons is a functional of charge density. The DFT potential is constructed as the sum of external potential due to atomic nuclei, which are seen as fixed by electrons within the Born-Oppenheimer approximation [14], and an effective potential due to electron interactions. The resulting electronic ground state is obtained by solving self-consistently one-electron Schrödinger-like equations known as Kohn-Sham equations [15]. Here, we used the DFT software VASP (Vienna Ab initio Simulation Package, version 5.4.4) [16,17], in which a plane-wave basis set and pseudopotential method are adopted. In our calculations, we employed the projector augmented wave (PAW) [18,19] method and the Perdew-Burke-Ernzerhof generalized gradient approximation (PBE-GGA) [20] functional. Charge carriers in the Re:5d 6 6s 1 and B:2s 2 2p 1 configurations were treated as valence electrons, and the valence wave functions were expanded in a plane wave basis up to a kinetic energy of 420 eV. The Monkhorst-Pack k-point sampling of the Brillouin zone [21] was chosen by a Γ-centered k-point mesh with a fine resolution = 0.01 × 2π/Å (33 × 33 × 13). The convergence criteria for self-consistent field and structure relaxation were set to 10 −6 eV/unit cell and 10 −3 eV/Å, respectively. For each given external pressure point, we first performed a structure optimization calculation in the hexagonal phase with fully relaxed lattice parameters and atomic positions. The theoretical lattice parameters at ambient conditions are a 0 = 2.913 Å and c 0 = 7.504 Å, which are within a 0.5% error margin compared to the corresponding experimental values. After the structure relaxation, we then performed calculations with lattice distortion to obtain the crystal's elastic tensor, which provided information on mechanical properties such as bulk and shear moduli, as well as crystal anisotropy. The bulk modulus computed by DFT with the Voigt-Reuss-Hill approximation [22] is K 0 = 357 GPa at ambient conditions, which agrees within a 2% error margin with the value K 0 = 364 GPa obtained by fitting the experimental P-V curve to the 3rd order Birch-Murnaghan equation. The theoretical structural visualization and charge distribution were plotted by the VESTA software (version 3.4.8) [23]. Figure 2 shows the integrated XRD powder data taken at the maximum pressure of 241 GPa with pressure determined using the platinum EOS [12]. The difference curve shown below the powder pattern in Figure 2 resulted from a fit to the hexagonal structure to ReB 2 . The hexagonal phase of ReB 2 was found to be stable to the maximum pressure of 241 GPa. The measured lattice parameters at 241 GPa were a = 2.586 ± 0.004 Å and c = 6.882 ± 0.007 Å. Platinum peaks in Figure 2 are labeled with asterisks (*) and indexed to a face-centered cubic lattice. The platinum lattice parameter at maximum pressure of 241 GPa was measured to be a = 3.490 ± 0.009 Å. Powder diffraction of ReB2 indexed to a hexagonal phase at a pressure of 241 GPa. Shown below the data curve is the difference curve as a result of Rietveld refinement. The platinum peaks (labeled with asterisk *) were indexed to a face-centered cubic phase and its measured volume was used in the calculation of pressure. Figure 3a shows the volume compression for ReB2 in two separate compression experiments with maximum pressures to 105 GPa in Experiment I (Expt. I) and 241 GPa in Experiment II (Expt. II). Figure 3b shows the volume compression for low-pressure data that exhibit uniaxial compression transitioning into non-hydrostatic compression at around 35 GPa. Both Expt. I and II showed transition zones from uniaxial compression to non-hydrostatic. For data points below the transition zone, the samples showed a fairly linear volumetric compression that is similarly seen in elastic samples before yielding to plastic deformation above the transition zone. It can be inferred from Figure 3 that the sample yielding at 30 GPa is a measure of the uniaxial compression strength, or the material's resistance to change before yielding. This is in agreement with the sample's average hardness, or its resistance to deformation being between 30 and 60 GPa [4][5][6][7]. The bulk modulus and its pressure derivative are taken from Equation (1) by fitting to the non-hydrostatic curve above 35 GPa, and they were determined to be K0 = 364 GPa and K0´ = 3.53, respectively.  Powder diffraction of ReB 2 indexed to a hexagonal phase at a pressure of 241 GPa. Shown below the data curve is the difference curve as a result of Rietveld refinement. The platinum peaks (labeled with asterisk *) were indexed to a face-centered cubic phase and its measured volume was used in the calculation of pressure. Figure 3a shows the volume compression for ReB 2 in two separate compression experiments with maximum pressures to 105 GPa in Experiment I (Expt. I) and 241 GPa in Experiment II (Expt. II). Figure 3b shows the volume compression for low-pressure data that exhibit uniaxial compression transitioning into non-hydrostatic compression at around 35 GPa. Both Expt. I and II showed transition zones from uniaxial compression to non-hydrostatic. For data points below the transition zone, the samples showed a fairly linear volumetric compression that is similarly seen in elastic samples before yielding to plastic deformation above the transition zone. It can be inferred from Figure 3 that the sample yielding at 30 GPa is a measure of the uniaxial compression strength, or the material's resistance to change before yielding. This is in agreement with the sample's average hardness, or its resistance to deformation being between 30 and 60 GPa [4][5][6][7]. The bulk modulus and its pressure derivative are taken from Equation (1) by fitting to the non-hydrostatic curve above 35 GPa, and they were determined to be K 0 = 364 GPa and K 0´= 3.53, respectively.

Results
Materials 2020, 13, x FOR PEER REVIEW 4 of 11

Figure 2.
Powder diffraction of ReB2 indexed to a hexagonal phase at a pressure of 241 GPa. Shown below the data curve is the difference curve as a result of Rietveld refinement. The platinum peaks (labeled with asterisk *) were indexed to a face-centered cubic phase and its measured volume was used in the calculation of pressure. Figure 3a shows the volume compression for ReB2 in two separate compression experiments with maximum pressures to 105 GPa in Experiment I (Expt. I) and 241 GPa in Experiment II (Expt. II). Figure 3b shows the volume compression for low-pressure data that exhibit uniaxial compression transitioning into non-hydrostatic compression at around 35 GPa. Both Expt. I and II showed transition zones from uniaxial compression to non-hydrostatic. For data points below the transition zone, the samples showed a fairly linear volumetric compression that is similarly seen in elastic samples before yielding to plastic deformation above the transition zone. It can be inferred from Figure 3 that the sample yielding at 30 GPa is a measure of the uniaxial compression strength, or the material's resistance to change before yielding. This is in agreement with the sample's average hardness, or its resistance to deformation being between 30 and 60 GPa [4][5][6][7]. The bulk modulus and its pressure derivative are taken from Equation (1) by fitting to the non-hydrostatic curve above 35 GPa, and they were determined to be K0 = 364 GPa and K0´ = 3.53, respectively.  The axial compression of the a and c lattice parameters for ReB 2 are shown in Figure 4, with the measured lattice parameters at 241 GPa being a = 2.586 ± 0.004 Å and c = 6.882 ± 0.007 Å. For both experiments, the c-axis showed a strong incompressibility as c/c 0 = 0.920 ± 0.001 at 241 GPa, not even 10% compression. In comparison, there was strong anisotropy between the a-axis and c-axis that persisted throughout the entirety of both experiments, and the anisotropy increased with pressure to maximum compression of a/a 0 = 0.891 ± 0.001. The maximum volume compression at 241 GPa was measured to be V/V 0 = 0.731 ± 0.004.
Materials 2020, 13, x FOR PEER REVIEW 5 of 11 The axial compression of the a and c lattice parameters for ReB2 are shown in Figure 4, with the measured lattice parameters at 241 GPa being a = 2.586 ± 0.004 Å and c = 6.882 ± 0.007 Å. For both experiments, the c-axis showed a strong incompressibility as c/c0 = 0.920 ± 0.001 at 241 GPa, not even 10% compression. In comparison, there was strong anisotropy between the a-axis and c-axis that persisted throughout the entirety of both experiments, and the anisotropy increased with pressure to maximum compression of a/a0 = 0.891 ± 0.001. The maximum volume compression at 241 GPa was measured to be V/V0 = 0.731 ± 0.004. To simulate the high-pressure experiments, we performed structure relaxation DFT calculations with the GGA functional up to 250 GPa. The DFT-GGA axial and volume compressions simulated under hydrostatic pressure are shown in Figure 5a,b, respectively. The a/a0 and c/c0 curves with pressure are both concave up, suggesting that the upturn or concave down behavior observed experimentally at pressures between 5 and 35 GPa (Figure 3b) is related to a non-hydrostatic condition. Figure 5a also shows that the lattice parameter a is more compressible than c. In particular, the DFT-GGA value c/c0 near 240 GPa is 0.919, which is in excellent agreement with the experiment. The DFT-GGA value a/a0 near 240 GPa is 0.877, which underestimates the experimental value of 0.891. In addition to non-hydrostatic condition, the theory-experiment deviation at high pressure is most likely due to the employed GGA functional. In particular, compared with previous local density approximation (LDA) studies at 100 GPa [24], while in both LDA and GGA values c/c0 = 0.954, the LDA ratio a/a0 = 0.932 is larger than the GGA value of 0.928. Regardless of the functional being employed, it is clear that an anisotropic compression behavior persisted up to the maximum pressure under study: a/c began with 0.388 at ambient conditions and decreased monotonically to 0.370 (0.375) in theory (experiment) near 240 GPa. To simulate the high-pressure experiments, we performed structure relaxation DFT calculations with the GGA functional up to 250 GPa. The DFT-GGA axial and volume compressions simulated under hydrostatic pressure are shown in Figure 5a,b, respectively. The a/a 0 and c/c 0 curves with pressure are both concave up, suggesting that the upturn or concave down behavior observed experimentally at pressures between 5 and 35 GPa (Figure 3b) is related to a non-hydrostatic condition. Figure 5a also shows that the lattice parameter a is more compressible than c. In particular, the DFT-GGA value c/c 0 near 240 GPa is 0.919, which is in excellent agreement with the experiment. The DFT-GGA value a/a 0 near 240 GPa is 0.877, which underestimates the experimental value of 0.891. In addition to non-hydrostatic condition, the theory-experiment deviation at high pressure is most likely due to the employed GGA functional. In particular, compared with previous local density approximation (LDA) studies at 100 GPa [24], while in both LDA and GGA values c/c 0 = 0.954, the LDA ratio a/a 0 = 0.932 is larger than the GGA value of 0.928. Regardless of the functional being employed, it is clear that an anisotropic compression behavior persisted up to the maximum pressure under study: a/c began with 0.388 at ambient conditions and decreased monotonically to 0.370 (0.375) in theory (experiment) near 240 GPa.

Discussion
The addition of interstitial covalently bonded boron atoms to high-electron-density transition metals such as Re and Os has given a family of transition metal diborides with desirable mechanical properties. Re and Os, being one column away from each other on the periodic table, share similar properties, although Re has a slightly smaller electron density and Os has a higher hardness and incompressibility [25,26]. OsB2 was shown to have a comparable bulk modulus (342-365 GPa) to ReB2, but ReB2 is considered somewhat superior due to shorter metallic bonds [27,28]. There is also a noticeable difference between the brittleness and ductility of the two materials. Pugh [29] introduced the ratio between the shear modulus and the bulk modulus (G/K) to distinguish a material's ductile or brittle behavior. A low (high) G/K value is correlated with ductility (brittleness). Based on the elastic and plastic properties of pure polycrystalline simple metals, an empirical value of the brittle-to-ductile transition is 0.571. Figure 6a shows the bulk and shear moduli computed by DFT using the GGA functional. While both G and K are enhanced by pressure, the rate of increase for K is larger, indicating that the G/K decreases with pressure, as seen in Figure 6b. In particular, G/K changes from 0.762 at 0 GPa to 0.627 at 240 GPa. In comparison, the G/K value is substantially smaller in OsB2 [30]. In our calculation, the G/K values for OsB2 at 0 GPa and 240 GPa are respectively 0.528 and 0.478, both of which are below the critical brittle-to-ductile transition value 0.571, showing that OsB2 is more ductile than ReB2. Another relevant quantity is the Poisson's ratio, which can be obtained by (3K − 2G)/[2(3K + G)]. Based on Frantsevich's rule [31], a material is brittle if its Poisson's ratio is less than 1/3; otherwise, the material is ductile. Figure 6b also shows that the Poisson ratio for ReB2 increases with pressure, indicating an enhanced ductility.

Discussion
The addition of interstitial covalently bonded boron atoms to high-electron-density transition metals such as Re and Os has given a family of transition metal diborides with desirable mechanical properties. Re and Os, being one column away from each other on the periodic table, share similar properties, although Re has a slightly smaller electron density and Os has a higher hardness and incompressibility [25,26]. OsB 2 was shown to have a comparable bulk modulus (342-365 GPa) to ReB 2 , but ReB 2 is considered somewhat superior due to shorter metallic bonds [27,28]. There is also a noticeable difference between the brittleness and ductility of the two materials. Pugh [29] introduced the ratio between the shear modulus and the bulk modulus (G/K) to distinguish a material's ductile or brittle behavior. A low (high) G/K value is correlated with ductility (brittleness). Based on the elastic and plastic properties of pure polycrystalline simple metals, an empirical value of the brittle-to-ductile transition is 0.571. Figure 6a shows the bulk and shear moduli computed by DFT using the GGA functional. While both G and K are enhanced by pressure, the rate of increase for K is larger, indicating that the G/K decreases with pressure, as seen in Figure 6b. In particular, G/K changes from 0.762 at 0 GPa to 0.627 at 240 GPa. In comparison, the G/K value is substantially smaller in OsB 2 [30]. In our calculation, the G/K values for OsB 2 at 0 GPa and 240 GPa are respectively 0.528 and 0.478, both of which are below the critical brittle-to-ductile transition value 0.571, showing that OsB 2 is more ductile than ReB 2 . Another relevant quantity is the Poisson's ratio, which can be obtained by (3K − 2G)/[2(3K + G)]. Based on Frantsevich's rule [31], a material is brittle if its Poisson's ratio is less than 1/3; otherwise, the material is ductile. Figure 6b also shows that the Poisson ratio for ReB 2 increases with pressure, indicating an enhanced ductility. We next address the strong lattice anisotropy observed in ReB2. Interestingly, the hexagonal ReB2 and orthorhombic OsB2 both showed anisotropic behavior in lattice parameters, with the c-axis being the most incompressible [5,27]. The observed extreme anisotropy of ReB2 shown in Figure 4 is likely attributed to the high electron density of Re and the high density of states (DOS) at the Fermi level (EF) [4], which result in increased Coulomb repulsion with pressure. Figure 7a shows our computed DOS for ReB2 at 0 GPa (top panel) and 240 GPa (bottom panel). The DOS plots indicate that at EF (denoted by the vertical dashed red line), the spectra have a dominant contribution from the Re atom. At high pressure, the spectral contributions at EF from Re and B atoms both increase, as seen in Figure 7b. The enhanced DOS also suggests an increased metallic bonding, or a reduced covalent bonding, which is consistent with the iso-surface charge density plots in Figure 8: at ambient conditions, ReB2 possesses strong covalent bonds between B-B atoms. When external pressure increases, the hybridizations between Re-B and B-B atoms both increase, leading to an enhanced DOS at EF and a reduced directional bonding. The increased Re-B bonding states near EF also can lead to enhanced bulk and shear moduli at high pressure [32].  We next address the strong lattice anisotropy observed in ReB 2 . Interestingly, the hexagonal ReB 2 and orthorhombic OsB 2 both showed anisotropic behavior in lattice parameters, with the c-axis being the most incompressible [5,27]. The observed extreme anisotropy of ReB 2 shown in Figure 4 is likely attributed to the high electron density of Re and the high density of states (DOS) at the Fermi level (E F ) [4], which result in increased Coulomb repulsion with pressure. Figure 7a shows our computed DOS for ReB 2 at 0 GPa (top panel) and 240 GPa (bottom panel). The DOS plots indicate that at E F (denoted by the vertical dashed red line), the spectra have a dominant contribution from the Re atom. At high pressure, the spectral contributions at E F from Re and B atoms both increase, as seen in Figure 7b. The enhanced DOS also suggests an increased metallic bonding, or a reduced covalent bonding, which is consistent with the iso-surface charge density plots in Figure 8: at ambient conditions, ReB 2 possesses strong covalent bonds between B-B atoms. When external pressure increases, the hybridizations between Re-B and B-B atoms both increase, leading to an enhanced DOS at E F and a reduced directional bonding. The increased Re-B bonding states near E F also can lead to enhanced bulk and shear moduli at high pressure [32]. We next address the strong lattice anisotropy observed in ReB2. Interestingly, the hexagonal ReB2 and orthorhombic OsB2 both showed anisotropic behavior in lattice parameters, with the c-axis being the most incompressible [5,27]. The observed extreme anisotropy of ReB2 shown in Figure 4 is likely attributed to the high electron density of Re and the high density of states (DOS) at the Fermi level (EF) [4], which result in increased Coulomb repulsion with pressure. Figure 7a shows our computed DOS for ReB2 at 0 GPa (top panel) and 240 GPa (bottom panel). The DOS plots indicate that at EF (denoted by the vertical dashed red line), the spectra have a dominant contribution from the Re atom. At high pressure, the spectral contributions at EF from Re and B atoms both increase, as seen in Figure 7b. The enhanced DOS also suggests an increased metallic bonding, or a reduced covalent bonding, which is consistent with the iso-surface charge density plots in Figure 8: at ambient conditions, ReB2 possesses strong covalent bonds between B-B atoms. When external pressure increases, the hybridizations between Re-B and B-B atoms both increase, leading to an enhanced DOS at EF and a reduced directional bonding. The increased Re-B bonding states near EF also can lead to enhanced bulk and shear moduli at high pressure [32].  and (2/3, 1/3, 1/2 ± z), where z is 0.0476 (0.0452) for P = 0 (240) GPa. The Re and B atoms are aligned along the c-axis. Under compression, as seen in Figure 8, the electron density is centered along the Re-B bonds, which are parallel to the c-axis. The strong electron Coulomb repulsion between charge density distributed along the Re-B bond direction makes the material highly incompressible along the c-axis. A strong anisotropy between the crystal a and c lattice parameters also suggests that the highest hardness in ReB2 single crystals is along the c-axis [5,33]. Iso-surface levels were set to be 0.1 a0 -3 , with a0 the Bohr radius. The calculations were based on the VASP software and its CHGCAR file, which contains the lattice vectors, atomic coordinates, the total charge density multiplied by the volume on the fine fast Fourier transform-grid, and the projector augmented wave one-center occupancies. The theoretical structural visualization and charge distribution were plotted by the VESTA software.
In addition to anisotropy in the lattice parameters, it is important to consider crystal elastic anisotropy, which is related to the occurrence of micro-cracks in materials [34,35]. Figure 9a shows the five independent elastic constants computed by DFT using the GGA functional for hexagonal ReB2 as a function of pressure. It is seen that c11 and c33 are largely enhanced upon compression, compared to the other elastic constants. Also, c33 is larger than c11, indicating that the c-axis is the least compressible. There are other ways to represent the level of elastic anisotropy of a material. In a hexagonal crystal, the following three parameters can be used [24,36]: ∆p = c33/ c11, ∆s1 = (c11 + c33 − 2c13)/4c44, and ∆s2 = 2c44/(c11 − c12). These three parameters would be equal to unity for isotropic compressibility. Figure 9b shows that the computed ∆p, ∆s1, and ∆s2 are all larger than 1 up to the maximum pressure under study. These results are consistent with previous lower-pressure theoretical studies [24,37], indicating a strong elastic anisotropy of ReB2, where its c-axis compressibility is smaller than that along the a-axis. the VASP software and its CHGCAR file, which contains the lattice vectors, atomic coordinates, the total charge density multiplied by the volume on the fine fast Fourier transform-grid, and the projector augmented wave one-center occupancies. The theoretical structural visualization and charge distribution were plotted by the VESTA software.
In the hexagonal phase of ReB 2 (with space group P6 3 /mmc), the two Re atoms are located at Wyckoff positions (1/3, 2/3, 1/4) and (2/3, 1/3, 3/4), and the four B atoms are located at (1/3, 2/3, ±z) and (2/3, 1/3, 1/2 ± z), where z is 0.0476 (0.0452) for P = 0 (240) GPa. The Re and B atoms are aligned along the c-axis. Under compression, as seen in Figure 8, the electron density is centered along the Re-B bonds, which are parallel to the c-axis. The strong electron Coulomb repulsion between charge density distributed along the Re-B bond direction makes the material highly incompressible along the c-axis. A strong anisotropy between the crystal a and c lattice parameters also suggests that the highest hardness in ReB 2 single crystals is along the c-axis [5,33].
In addition to anisotropy in the lattice parameters, it is important to consider crystal elastic anisotropy, which is related to the occurrence of micro-cracks in materials [34,35]. Figure 9a shows the five independent elastic constants computed by DFT using the GGA functional for hexagonal ReB 2 as a function of pressure. It is seen that c 11 and c 33 are largely enhanced upon compression, compared to the other elastic constants. Also, c 33 is larger than c 11 , indicating that the c-axis is the least compressible. There are other ways to represent the level of elastic anisotropy of a material. In a hexagonal crystal, the following three parameters can be used [24,36]: ∆ p = c 33 / c 11 , ∆ s1 = (c 11 + c 33 − 2c 13 )/4c 44 , and ∆ s2 = 2c 44 /(c 11 − c 12 ). These three parameters would be equal to unity for isotropic compressibility. Figure 9b shows that the computed ∆ p , ∆ s1 , and ∆ s2 are all larger than 1 up to the maximum pressure under study. These results are consistent with previous lower-pressure theoretical studies [24,37], indicating a strong elastic anisotropy of ReB 2 , where its c-axis compressibility is smaller than that along the a-axis.

Conclusions
Ultrahigh non-hydrostatic compression studies were carried out on a superhard material ReB 2 for the first time to a pressure of 241 GPa. The equation of state determined from the non-hydrostatic pressure-volume curve above 35 GPa yielded a bulk modulus and pressure derivative of K 0 = 364 GPa and K 0´= 3.53, respectively. Substantial anisotropy of the lattice parameters was indicated to increase with pressure up to the maximum pressure, with a/a 0 = 0.891 and c/c 0 = 0.919 at 241 GPa, showing~3% difference in axial compression. The results from density functional theory simulations for anisotropic compression, equation of state, and elastic constants were in good agreement with the experimental data. The superhard and ultra-incompressible features of ReB 2 render it a promising material for wide ranges of applications in extreme environments.