Introduction

Superhard materials are essential for the automotive, oil and gas, aerospace, and any manufacturing industry that relies on drilling, cutting, and grinding [1–4]. These materials have an experimentally measured Vickers hardness (H v,exp) exceeding 40 GPa [5, 6] and fall into two general classes. The first contains only light main group elements with diamond being the prototypical example due to its outstanding mechanical properties (H v,exp ≈ 96 GPa). Similarly, cubic boron nitride (c-BN) is isoelectronic to diamond, possesses the same crystal structure, and is superhard (H v,exp ≈ 50 GPa) [7]. It also has a distinct advantage over diamond since it can be used to cut ferrous materials without the formation of undesired carbide side-products [3, 8, 9]. Other compounds with a high hardness discovered in this class include B 6O and BC 2N with a H v,exp of 45 and 76 GPa, respectively [10, 11]. In all of these examples, the highly directional, covalent bonding networks are essential for the necessary mechanical response. Nevertheless, their synthesis requires high temperature and high pressure (HTHP) for preparation limiting the wide-scale adoption of most compounds in this class.

Alternatively, combining light elements such as B, C, N, and O with 4d or 5d transition metals comprises the second class of materials with potential for superhardness [9]. The origin of this mechanical behavior stems from the heavy transition metals having an outstanding bulk modulus due to their high density of valence electrons; however, the inherent delocalization of the electrons in the transition metals adversely affects hardness. Incorporating light main group elements with their propensity to form short, directional covalent bonds allows enhanced mechanical response [9]. For example, transition metal borides like ReB 2 [12], OsB 2 [13], and RuB 2 [14] all show ultra-incompressibility and high hardness. Re 2C also yields moderate hardness (H v,exp = 17.5 GPa) [15] while transition metal nitrides like OsN 2, IrN 2 [16], and Re 3N [17] also demonstrate respectable properties. Many of these compounds are approaching the threshold for industrial application while the synthesis most often tends to only require high temperatures making them highly desirable [18, 19].

This second class of materials clearly displays fascinating mechanical response and synthetic accessibility; yet, the availability of the starting elements in many of these compounds is of great concern. The most common transition metals currently used in this class of compounds are W, Re, and Os where the natural (crustal) abundance ranges between 7 ×10 −4 and 1 ppm [20]. This is in stark contrast to the first class of materials that contain main group elements (B, C, N, O, and Si), which are generally considered abundant. Switching to early transition metals like Zr and Y or most 3d transition metals that have crustal abundances between 30 and 5000 ppm is the best option to address sustainability concerns. Additional considerations must also be made for the synthetic routes used in material preparation. For example, the first class of materials requires huge amount of energy for HTHP synthesis that leads to various sustainability and scalable production concerns [21–23]. Instead, a majority of the second class materials, only requires high temperatures for synthesis, reducing energy consumption and eliminating specialized presses or high pressure reactors for preparation.

In the search for new superhard materials, a shift toward the synthesis of cheap, earth-abundant materials that do not require extreme synthetic conditions is required. Yet, identifying the phase-space containing materials that balance the mechanical response with resource and synthesis/processing considerations is not often straightforward. One solution is to capitalize on computational or informatics methods to predict mechanical properties and augment these results with compositional and synthetic information. Combining these distinct foci will allow researchers to predict complex mechanical response like hardness while targeting only the ideal composition space. An informatics-based approach has already been successfully employed in the development of other functional materials like thermoelectrics and batteries [20, 24, 25]. Assembling a database for mechanical materials with an associated interactive website allows a high-information density visualization method to identify correlations between structural and compositional characteristics and the mechanical properties, ultimately directing synthesis. This visualization technique is akin to the Ashby diagrams, which have been widely employed in the materials selection and design optimization process [26]. Performance indices for a given application are developed from materials properties, or functions and combinations thereof, and plotted for different families of material for comparison. Further, introducing resource considerations like elemental abundance, elemental availability, and even reaction conditions, e.g., temperature and pressure, in to this database will ensure the products are not only viable mechanical materials but also that the targeted compositions are synthetically accessible and sustainable.

Results and Discussion

Screening for High-Hardness Materials

To facilitate the search for novel superhard materials, researchers would ideally like to explicitly calculate hardness using ab initio computation; however, there are two primary obstacles impeding such an approach. First, hardness is influenced by structural components across multiple length scales including, dislocations, impurities, domains, and grains. Second, hardness is not a fundamental material property, but rather a composite property influenced by yield strength, work hardening, ultimate tensile strength, modulus of elasticity, and others. Instead researchers must rely on correlating the crystal structure or electronic structure to the experimentally measured hardness. Methods have already been developed with this approach in mind. For example, it is possible to analyze the sum of individual bond strength by calculating the energy required to promote a valence electron into the conduction band (Fig. 1) [27] or by comparing electronegativities (Fig. 1) [28] with both providing decent correlation to hardness. Alternatively, (intrinsic) hardness can be determined through a measure of valence electron count, interatomic distances, coordination number, and the atomic radius of the component elements (Fig. 1) [28, 29]. These assumptions are reasonable for purely covalent materials while correction factors for ionic bonding or small metallic components can also be accounted for to improve the models and make them more transferable (Fig. 1) [27, 29–31]. Such correction factors are particularly important in metallic-type and polar covalent materials, which make up a large number of the known high-hardness materials. Combining all of these different models in Fig. 1 highlights that the calculated hardness (H calc) using four independent methods compared to H v,exp shows a positive correlation but with a wide scatter, in particular at low hardness. Because these models are all developed based on the assumption of bond localization, they agree well within the scope of covalent-dominant materials, like diamond. Of considerable importance is that these (mostly) semi-empirical calculations used to determine hardness can likely be scaled to examine vast composition space. Nevertheless, many of these relations do not hold for metallic systems due to the delocalization of electrons and difficulty describing the chemical bonding in metals. These models also fail to describe the hardness of highly anisotropic materials, disconnected structures (structures with lower than four coordination numbers), and molecular structures where van der Waals forces exist, degrading their universal predictive ability [29, 32].

Fig. 1
figure 1

Comparing available models for predicting hardness. Comparison between calculated (intrinsic) hardness using multiple semi-empirical models and the experimentally measured Vickers hardness (\(H_{\mathrm {v},\exp }\)) demonstrate reasonable agreement within the scope of studied compounds. H calc is derived from references [27–30, 37] (Color figure online)

Full size image

Instead of relying on semi-empirical based calculations, it is possible to use ab initio calculations to determine a material bulk modulus (B) and shear modulus (G) from high-level computations of atomic forces. These intrinsic material properties are related to hardness [33] and therefore can act as a proxy or screening parameter. For example, as shown in Fig. 2a, plotting the experimental Vickers hardness (H v,exp) as a function of B indicates minor correlation between these two properties. Plotting H v,exp against G (Fig. 2b) results in a nearly linear correlation. Examining Young’s modulus (E), which incorporates B and G, also leads to a linear relationship with the H v,exp, shown in Fig. 2c. From these data, it is evident that shear and Young’s modulus both demonstrate the strongest correlation with experimental hardness; this suggests the feasibility of relating hardness to elastic moduli, which can be easily understood and readily calculated using ab initio methods. A more robust agreement is achieved by plotting H v,exp against a parameterized function involving G and B following the equation f(G/B) = 0.92k 1.137 G 0.708, where k is Pugh’s ratio (G/B) [34, 35]. An excellent match is obtained between the f(G/B) and the experimentally measured hardness, as shown in Fig. 3. Interestingly, Pugh’s ratio was originally developed to distinguish between ductility versus brittleness. The high value of k corresponds to brittle behavior and vice versa. Albeit brittleness is most often associated with high hardness (also evident from the equation); it adversely affects toughness or the ability of a material to absorb energy and is evaluated via impact test, which could be of great application importance [36]. Furthermore, the deformation of extremely ductile solids like gold is governed by plastic flow and dislocation interactions, making it exceedingly difficult to model [33]. As a result, this model does not perform well at low hardness, e.g., < 5 GPa, or for extremely ductile metals. Nevertheless, this model performs well across a wide range of hardness values. The high degree of correlations among the compositional properties like electronegativities or a material’s elastic properties means there are numerous methods for predicting hardness [27–30, 37]. Employing such a broad set of material considerations allows candidate compounds available in large databases such as the Inorganic Crystal Structure Database (ICSD) to be screened for desirable properties thereby focusing synthetic efforts. It is straightforward to include elemental parameters like electronegativity to sort materials in the ICSD into those with potential for high hardness. However, the best agreement is achieved by using the elastic moduli like G and B. Fortunately, first principles calculations can accurately reproduce the stress-strain relationship, which in combination with the Voigt-Reuss-Hill approximation generates these constants [38–40]. Thus, the high-throughput computation and correlation diagrams are valuable to decrease the number of exploratory synthetic experiments conducted, thereby improving discovery rates while reducing overall waste generated. The utility of calculating reliable metrics is realized by developing databases of material properties generated using high-throughput computation within the framework of the Materials Project [41, 42]. In this case, using ab initio calculations of G and B to determine f(G/B) as a proxy for H v,exp, over 1100 compounds can have their potential intrinsic hardness estimated. As shown in Fig. 4, compounds with f(G/B) higher than 20 GPa (dashed line) are considered compounds of interest and worthy of experimental consideration. As summarized in Table 1, some specific compositions indicated by the database are HfBe 5 and MnBe 2, which both demonstrate potential high hardness with an f(G/B) of 25.28 and 33.54 GPa, respectively. More accessible binary compositions that should be considered are Nb 6 C 5 (f(G/B) = 24.46 GPa) and NbIr 3 (f(G/B) = 28.14 GPa). In addition, there are also noteworthy ternary compositions including Hf 2Al 3 C 4 (f(G/B) = 27.37 GPa), ScAl 3 C 3 (f(G/B) = 25.99 GPa), and Sc 3FeC 4 (f(G/B) = 20.0 GPa). These examples demonstrate that it is possible to target compositions using high-throughput first-principle computation thereby improving the rate of novel hard materials discovery.

Table 1 Some materials worth investigation and their corresponding space group predicted based on their high f(G/B)
Full size table
Fig. 2
figure 2

Elastic moduli scale with hardness. Correlation of bulk modulus (B) (a), shear modulus (G) (b), and Young’s modulus (E) (c) with experimental Vickers hardness (\(H_{\mathrm {v},\exp }\)). The dashed line is shown as a guide to the eye (Color figure online)

Full size image
Fig. 3
figure 3

A function of the elastic moduli provides the best agreement. A proposed function to model hardness based on a function of G and B (f(G/B)) is in agreement with experimentally measured Vickers hardness (\(H_{\mathrm {v},\exp }\)). The dashed line is the function f(G/B) [35] (Color figure online)

Full size image
Fig. 4
figure 4

Screening for high-hardness materials. Calculated f(G/B) for the available database. The compounds above the dashed line indicates compounds with a probable f(G/B) > 20 GPa. The inset expands the range above 20 GPa and highlights compositions of interest with an unreported \(H_{\mathrm {v},\exp }\) (Color figure online)

Full size image

Merging Resource and Synthesis Considerations with Mechanical Properties

Screening for materials’ mechanical properties using a combination of materials informatics and first-principle computation is the first step in the development of novel high-hardness materials. It is also critical to incorporate resource considerations including sustainability of compositions as well as the material preparation routes to ensure sustainable material development.

In terms of chemical composition, scarcity (inverse crustal abundance, ζ) and the supply and demand risk (Herfindahl-Hirschman index), based on geological and geopolitical data, provide a quantitative measure of resource economic and sustainability factors. Using the United States Geological Survey (USGS) commodity statistics, the HHI parameter can be calculated as sum squared of market fraction (x i ) for a given country, based on their production (HHIP) or geological reserves (HHIR) of each element [20]. Here, for each composition, the weighted average HHIP values were calculated using weight fraction of each element in the chemical formula. The U.S. Department of Justice and the Federal Trade Commission define markets as unconcentrated, highly concentrated, or moderately concentrated for a given commodity when HHI values are below 1500, over 2500, and between these limits, respectively [43, 44]. Borrowing crustal abundance from the CRC Handbook elemental scarcity, [45] which is crustal abundance in inverse parts per million (ppm −1), a weighted average for each compound’s scarcity was calculated based on the weight fraction of elements in the chemical formula [20]. Known materials with a range of hardnesses can first be analyzed using these data by developing high-information density plots. In this case, Fig. 5 illustrates elemental scarcity and HHIP, with the marker size proportional to H v,exp and the marker color corresponding to different classes of superhard materials. The first class is most often HTHP synthesis while the second class corresponds to mainly high temperature. It is immediately evident that most known high-hardness materials possess scarcity of over 10000 ppm −1, and they are classified as materials with at least a moderate supply and demand risk. For comparison’s sake, this means that most of these compounds are more scarce than nickel (≈11000 ppm −1). It is possible to use Fig. 5 and simply look for compounds with large markers located in the lower-left corner of the plot. This would suggest the material has a high hardness and is both abundant and readily available. Such an analysis indicates that main group materials such as diamond, c-BC 5, c-BN, and BC 2N are the best choices with TiC, SiC, and TaB 2 as alternative abundant materials even though they are significantly less hard. This analysis also shows that compositions like WB 4 do not have a high scarcity; however, this composition contains elements that make it highly concentrated and possibly at risk. Many of the transition metal borides and nitrides show even worse indicators than WB 4 with a ζ that varies between 10000 and 100000 ppm −1 due to the presence of transition metals such as Re, Ru, and Os despite achieving hardnesses as large as 30 GPa. This sorting diagram can be expanded to also select optimal compositions from the computational database. Figure 6 plots ζ against the HHIP with the marker size for the compounds with a known H v,exp (gray circles) and the marker size of the predicted materials scaling with f(G/B) (red circles). As highlighted by the vertical dashed line compositions to the left in the plot, all contain elements that do not possess a supply and demand risk. Interestingly, nearly all of these compositions are sustainable with a scarcity index > 1 ×10 6 ppm −1. For example, the binary composition VAl 3 (f(G/B) = 20 GPa) has a scarcity of 2.5 ×10 3ṗpm −1 and an HHIP of 2.2 ×10 3 making it more abundant than compounds like c-BN and less of an economic risk than SiC. The diagram also highlights ternary phases like Ti 2AlC, which has a f(G/B) = 22 GPa and has a scarcity and HHIP nearly the same as TiC. Employing this visual sorting diagram clearly makes selecting novel, sustainable, high-hardness compounds self-evident.

Fig. 5
figure 5

Analyzing sustainability parameters. Elemental scarcity (ζ) is plotted against HHIP with marker size scaling with reported experimental hardness (\(H_{\mathrm {v},\exp }\)). The color of the symbol indicates material class with the first class being HTHP synthesis (yellow) and the second class mostly requiring only high temperature (blue) (Color figure online)

Full size image
Fig. 6
figure 6

Screening for sustainable compositions. Elemental scarcity (ζ) is plotted against HHIP. The red marker size scales with the predicted hardness based on f(G/B) and the gray circles are the experimental Vickers hardness (\(H_{\mathrm {v},\exp }\)). Left of the vertical dashed line are the ideal compositions (Color figure online)

Full size image

Nevertheless, these plots have several shortcomings. First, it implies that scarcity and HHIP are equally important when selecting superhard materials. This may not be the case. Although high market concentration can lead to commodity speculation, supply risk, and increased cost, the primary resource driver of cost is scarcity, and therefore, it should receive greater weight in material selection. Second, the analysis from Fig. 5 alone ignores difficulty in synthetic preparation. In fact, some of the hardest materials require more extreme HTHP routes to be synthesized. Clearly, in screening for superhard materials, there is need to quantitatively account for difficulty in material preparation alongside hardness indicators and resource considerations such as scarcity or HHIP. One way to achieve this is by including synthesis variables such as temperature and pressure in the high-information density plots. One such example is shown in Fig. 7, where scarcity is plotted against calculated f(G/B) with marker size now scaling with synthetic pressure and color indicating synthetic temperature. Replacing H v,exp with f(G/B) provides the opportunity to predict sustainable, hard materials by directly calculating the elastic modulus of crystal structure, which allows property prediction a priori. Materials on the right hand side of the resulting plot have the largest f(G/B) and therefore higher hardness. It is also clear that the application of both high pressures and high temperatures is required to prepare these materials. For example, c-BC 5 is synthesized at the extremely severe conditions of 2200 K and 24 GPa. Even diamond, although better in comparison to c-BC 5, γ-B 28, and BC 2N, still requires specialized sintering hot presses to achieve 1473 K and 5 GPa during synthesis. This in turn leads to the high cost of synthetic diamond cutting tools despite carbon itself being abundant with unconcentrated production. The alternative is to focus on high-hardness materials that can be synthesized via pressureless sintering methods, induction heating, or arc-melting. These compounds including OsB and RuB 2 need the highest temperatures while RuO 2 and InSb requires the lowest temperature. Materials that offer a compromise between good hardness, scarcity, and synthetic conditions are TiB 2, ZrB 2, and BP. Of these materials, we find that BP and TiB 2 benefit from only moderate supply risk.

Fig. 7
figure 7

Incorporating material preparation parameter space. Scarcity plotted against f(G/B) with marker size scaling with synthetic pressure and marker color corresponding to synthetic temperature. A solid color fill represents no supply risk (unconc.) via HHIP, a gray fill is a moderate risk, and no fill is high risk (conc.) (Color figure online)

Full size image

Conclusions

Merging computational proxies for a hardness and basic materials informatics through the construction sorting diagrams allows researchers to scan vast inorganic crystal structures databases using first-principle calculations. This method will expedite the search for superhard materials by targeting ideal composition and crystal structures. Further, incorporating economical and environmental consideration into these databases is essential to narrow the research focus to abundant and environmentally benign elements. Combining all of the resulting data on high-information density plots generate an ideal tool that clearly illustrates the balance between material performance and incorporating only elements that are not geologically or geopolitically at risk. These databases can also be expanded to include synthetic accessibility; however, this approach is limited to compiling experimental reports of reaction conditions rather than based only on ab initio calculations. Together, this research provides a unique method for directing future synthetic efforts by identify trends in material class and composition as well as observe potential new “design rules” for new superhard materials.

Methods

The dataset used in this work including the mechanical properties and reaction conditions, e.g., temperature and pressure, were manually extracted from published literature [33, 41, 42]. For graphical data representations in figures, for example, data mining was realized using free software such as PlotDigitizer [46] and DataThief [47]. HHI values based on production and reserves for each element were calculated using USGS data from commodity statistics factsheets from 2011 where available and from 2010 or 2009 where unavailable [20, 48, 49]. The crustal abundance for each element was adapted from the CRC Handbook [45]. Batch calculations of scarcity and availability (HHI) were made using a web tool developed during this work hosted at http://www.scarcitycalculator.com.

The elastic moduli can be accurately reproduced from first-principle calculations using the stress-strain relationship in combination with the Voigt-Reuss-Hill approximation [38–40]. The elastic constants (C i j are determined using the Vienna ab initio simulation package (VASP) [50, 51] and the arithmetic relations outlined by Wu et al. [52]. All VASP calculations were performed using the projector augmented wave (PAW) [53] pseudopotentials with exchange and correlation described by the Perdew-Burke-Ernzerhof generalized gradient approximation (PBE-GGA) [54]. The structures were initially fully relaxed (c/a ratio, unit cell volume, and atomic positions) to ensure the crystal structure is at an electronic ground state. The elastic tensors were then determined based on six finite distortions of the crystal using displacements of ±0.015 Å. The energy cutoff of the planewave basis set was fixed to 700 eV, and a minimum of 100 k-points were used to ensure electronic convergence. Details of the high-throughput computations adapted from the materials project can be found elsewhere [41, 42].

Data processing and visualization is achieved using a method originally described by Gaultois and coworkers for thermoelectric materials [20]. In that work, high-information density plots were created from data sets using the commercial Highcharts code by encoding vectorized information as abscissa, ordinate, marker size, and marker color. Seshadri and Sparks then described the creation of a generic visualization tool which was used for the present work [55].