A Modified Theoretical Model of Intrinsic Hardness of Crystalline Solids

Super-hard materials (H_v > 40 GPa) have attracted great attentions because of their practical importance with numerous industrial applications, e.g. cutting and polishing tools, wear resistant coatings, and abrasives^1,2. Over the past several decades, great efforts have been devoted to explore new super-hard materials^1,2,3,4.

Historically, searching for super-hard materials was guided under correlations between hardness and macroscopic properties, such as bulk modulus B⁵ or shear modulus G⁶. However, these empirical correlations are physically questionable, which usually result in misleading. Even though Chen et al.⁷ elegantly modeled the hardness of a material by introducing the Pugh’s ratio⁸, k = G/B, which characterizes the brittleness/ductility of the material, these macroscopic concepts cannot provide insight into the physical origin of hardness. As a consequence, estimating the hardness directly with microscopic parameters may reveal the fundamental factors that control hardness, which is essential for the design and exploration of new super-hard materials.

Up to now, three categories of microscopic models for hardness evaluation have been proposed^{3,4,9,10,11,12}. All of the microscopic models share the same assumption that the hardness equals to the sum of resistance of each bond per unit area to the indenter⁹ with the resistance estimated under different hypotheses^3,4. The microscopic model proposed by Gao at al.⁹ is probably the most popular one with the resistance assumed proportional to the homo-polar energy gap. Later, Gao¹⁰ suggested that the strength of a bond can be more accurately characterized by using average overlap populations per unit volume of the bond, where the overlap populations is evaluated from first-principles calculations. With the advancement on computational technology, the newly developed formula by Gao¹⁰ is getting more and more attractive. However, recently, we found that the hardness values of materials with complex bond types estimated from the formula exhibit dramatically discrepancies with experimental measurements. In the present work, we will reexamine the formula in detail and modify it to capture as many materials as possible.

In microscopic models, it is usually assumed that any complex crystal can be decomposed into a set of pseudo-binary crystals (chemical bonds), and properties of the crystal can be derived from readily accessible parameters associated with chemical bonds, e.g. bond length, valence electron number, ionicity and etc. Strictly speaking, a pseudo-binary crystal represents a pair of neighboring atoms in the crystal, which is different from the traditional concept of a chemical bond, since complex bonding may form between atoms in a crystal. For example, both σ bond and π bond form between carbon atoms in graphite. Nevertheless, we will still call a couple of interacted neighboring atoms as a chemical bond for simplicity in the present work.

The hardness of a multi-component compound can be expressed as the geometrical average of hardness of all pseudo-binary crystals that comprise the compound⁹:

where H^μ and N^μ respectively represents the hardness and number of the μ type bond in the compound.

The hardness of μ type bond equals to the resistance of the bond per unit area to the indenter¹⁰:

where A is a proportional coefficient, N_a is the covalent bond number per unit area, P^μ is the Mulliken overlap population¹³, and v^μ is the bond volume. In the equation, P^μ/v^μ characterizes the strength of the bond, as suggested by Gao¹⁰. For a specific bond, the bond number per unit area N_a is evaluated from its electron number per cubic angstroms as (N_e^μ/2)^2/3, where N_e^μ equals to n^μ/v^μ with n^μ being the number of electrons of the bond. Substitution N_a into equation (2), one obtains:

However, in the original work by Gao¹⁰, H^μ was expressed as:

It reveals that equation (4) is over-simplified, since n^μ of a particular bond is not always 2. To make sure that the model is applicable to crystals comprised of complex bond types, further modifications are still necessary, which will respectively be taken on n^μ, v^μ and the averaging process.

Previously, the number of valence electrons n^μ of a bond follows¹⁴:

where Z_A^μ and N_CA are the valence electron number and coordination number of the A atom constructing μ type bond, respectively. Z_B^μ and N_CB are in analogous to Z_A^μ and N_CA. Equation (5) explicitly assumes that the valence electron of atom A is equally partitioned to bonds surrounding it. This assumption is only reasonable when bonds surrounding the atom are similar in nature. However, in a crystal comprised of complex bond types, the equal partition may result in an unrealistic n^μ that deviates significantly from the true valence electron number of the bond. Take TiB₂ as an example. The electron redistribution map in Fig. 1b displays the characteristics of B-B σ bond, B-B π bond and Ti-B ionic-covalent bond¹⁵. In principle, each B-B σ bond contains 2 valence electrons and B atoms need extra electrons transferred from Ti to form π bond. The rest valence electrons of Ti form ionic-covalent bond between Ti and B. If we assume that valence electrons of Ti are equally partitioned to 3 B-B π bond and 12 Ti-B bond surrounding the Ti atom, then valence electron number for B-B and Ti-B bond is respectively 34/15 = 2 + 4/15 and 4/15. However, if we use equation (5), valence electron numbers for both bonds are 2/3 with N_B = 9 and N_Ti = 12. Obviously, the valence electron number derived from equation (5) is unrealistic, especially when bonds surrounding an atom are significantly different in nature, or electron transfer is involved during the formation of a crystal, or non-bonding electrons (lone pair electrons) exist.

(a) Crystal structure of TiB₂. (b) Electron density difference map, which represents charge redistribution due to formation of chemical bonds, on () plane to illustrate the bonding nature of TiB₂. The map displays the characteristics of B-B σ bond, B-B π bond and Ti-B ionic-covalent bond. For example, the strong B-B σ bond results in a substantial redistribution of electrons into the space between B-B, while the π bond is in accordance with the shoulder to shoulder type of electron redistribution.

Full size image

The second modification is made on the definition of bond volume v^μ. The original definition of bond volume is introduced by Levine¹⁶ in 1973, which assumes that the volume of a bond is proportional to (d^μ)³. d^μ is the bond length. This definition has been broadly adopted in models associated with chemical bonds for over four decades. Here, we suggested that the bond volume (influence region of a bond) is not just correlated to its length, but also proportional to its valence electron number. Accordingly, the cell volume is partitioned to the component chemical bonds with v^μ being:

where Ω represents the cell volume of the crystal.

The third modification is conducted on the averaging process. In equation (1), the hardness is averaged over the number of bonds. Similar to the modification in bond volume, the valence electron number is also emphasized in the averaging process as follows:

Both modification on bond volume and the average process emphasize the crucial role of valence electron populations, which means that both decomposing the crystal into pseudo-binary crystals and estimating properties of the crystal from properties of pseudo-binary crystals should be based on valence electron populations of different bonding states. Comparing to the purely geometrical considerations based on crystal structures, taking the population of valence electrons into consideration is more physical in nature, while it is well-known that properties of a crystal depend strongly on bonding states of valence electrons. For a deep understanding on the modifications, take n^μN^μ as a whole, which accounts for a group of valence electrons occupying the same electron bands. A set of specific electron bands of a crystal are comparable to molecular orbits of a specific chemical bond. In general, a well-defined chemical bond is occupied by 2 valence electrons. Then, we can define M^μ = n^μN^μ/2 as the number of equivalent chemical bonds formed by these valence electrons. At this circumstance, the Mulliken overlap population of an equivalent bond is Q^μ = N^μP^μ/M^μ = 2P^μ/n^μ. Substitution Q^μ and M^μ into equation (4) and the original definition of v^μ results in the hardness of the equivalent chemical bond H^μ being:

Equation (8) is exactly the same as equation (3) with the bond volume defined in equation (6). In addition, take equation (8) and M^μ into equation (1) results in equation (7). It indicates that the modifications are self-consistent and equivalent to decomposing the crystal into a set of well-defined chemical bonds. It is noteworthy that equation (8) reduces to equation (4) that was proposed in the original model¹⁰, when all the bonds in the crystal are well-defined chemical bond, i.e. 2 valence electrons per pseudo-binary crystal. According to equation (8), hard chemical bonds need high covalency of the bond (high 2P^μ/n^μ), short bond length and high valence electron density. These conditions are consistent with other microscopic models^9,10,11,12, which can facilitate the quest for intrinsic super-hard materials.

Crystals from refs 7,9, 10, 11 are selected to check the availability of the modified model. These crystals are classified into three groups, which are corresponding to crystals with zinc blende or wurtzite structure, rock salt structure, and other complex structures. The results are listed in Tables 1, 2 and 3, respectively. The Mulliken overlap population P^μ of a bond was evaluated using first-principles calculations by CASTEP¹⁷. The Vanderbilt-type ultrasoft pseudopotential¹⁸ and exchange-correlation described by generalized gradient approximation¹⁹ were employed. The plane wave cutoff energy was set to be 500 eV. k-points mesh with a separation of 0.03 Å⁻¹ according to Monkhorst-Pack method²⁰ was adopted in the Brillouin zone. For each crystal, the structure was optimized and compared with experimental data to confirm the reliability of the calculation.

In order to determine the coefficient A in equation (4), theoretical values H_T/A determined from crystal structures and bonding properties versus experimental hardness values H_E are plotted in Fig. 2. It is clear that all three sets of data locate at a straight line passing through the origin. By fitting, A is determined to be 693 with R² = 0.984, which is close to the value of 740 suggested by Gao¹⁰ that was derived from the hardness value of diamond. The coincidence of A is not surprising, since the modified model reduces to the original one when n^μ of each bond equals to 2. For simplification, A is taken as 700 in the future. The theoretical hardness values derived from the modified model are also listed in Tables 1, 2 and 3 for comparison, where good agreement is obtained.

Linear fit of the modified model.

Full size image

To further verify the capability of the modified model, it is applied to investigate the hardness of TMB₂s (TM = Ti, Zr, Hf, Re and Os). TiB₂, ZrB₂ and HfB₂ have a simple hexagonal structure (space group P6/mmm), where TM and B atoms are respectively occupy 1a(0, 0, 0) and 2d(1/3, 2/3, 1/2) Wyckoff sites, as shown in Fig. 1a. ReB₂ has a simple hexagonal structure (space group P6₃/mmc), where Re and B atoms respectively occupy 2c(1/3, 2/3, 1/4) and 4f(1/3, 2/3, 0.548) Wyckoff sites, as shown in Fig. 3a. OsB₂ has an orthorhombic structure (space group Pmmn), where Os and B atoms respectively occupy 2a(1/4, 1/4, 0.154) and 4f(0.058, 1/4, 0.632) Wyckoff sites, as shown in Fig. 3b. Different from crystals in Tables 1, 2 and 3, where chemical bonds are well-defined, decomposing these TMB₂ into different kinds of pseudo-binary crystals is not intuitive. Analysis on the decomposition is guided by P^μ. Any pair of atoms with positive P^μ is assumed to be an effective pseudo-binary crystal. Take ReB₂ for example. According to P^μ, there are four types of bonds in ReB₂ (Table 4), B-B bonds, two types of Re-B bonds and Re-Re bonds. The B-B bonds are typical covalent σ bonds with 2 valence electrons per bond. To specify the electron number of other chemical bonds, it is assumed that valence electrons of Re is equally partitioned to all Re-B bonds and Re-Re bonds surrounding it. Therefore, each bond shares 1/2 electron from each Re atom, which means 1/2 electron per Re-B bond and 1 electron per Re-Re bond. In analogous, B-B bonds in OsB₂ are assumed to be typical covalent σ bonds with 2 valence electrons, while Os-B bonds share valence electrons from Os resulting in 1 electron per bond. Decomposition of MB₂ (M = Ti, Zr, Hf) has been introduced above during modification of n^μ, which will not be repeated. With an appropriate allocation of valence electrons to the chemical bonds, the hardness can be predicted by the modified model. The predictions agree well with experimental measurements for these TMB₂s, as shown in Table 4.

Crystal structure of ReB₂ and OsB₂.

Full size image

To verify the improvement of the modifications, hardness values predicted by using the original Gao’s model and the modified model were compared with the experimental measurements, as shown in Fig. 4. In the calculation by using the original model, the proportional coefficient A is also adopted as 700 instead of 740. As stated above, when all the decomposed pseudo-binary crystals occupy 2 valence electrons, the modified model reduces to the original one. Therefore, results for crystals in Tables 1 and 3 obtained from both models are the same, as demonstrated in Fig. 4 that results predicted by different models overlap with each other. However, hardness values for crystals in Table 2 predicted by the modified model are lower than those from the original model. As shown in Fig. 4, without any modification, hardness values of transition metal carbides and nitrides predicted by the original model exhibit a systematic over-estimation. For a rock-salt structure crystal, its unit cell is decomposed into 6 equivalent pseudo-binary crystals with n^μ less than 2. Therefore, hardness values predicted by the modified model will be (n^μ/2)^2/3 times those predicted by the original model. After modification, the predicted values agree well with experiments. Figure 4 also reveals that the predicted hardness values for those TMB₂s are also significantly improved after modification.

Calculated hardness values from models with and without modifications versus the experimental values.

Solid points represent values calculated from the previous Gao’s model, while hollow points are calculated from the modified model. The dash line means H_T = H_E.

Full size image

Before ending, some fundamental aspects on hardness are discussed. It should be noted that experimental measured hardness values usually exhibit significant divergence, since the measurements are very sensitive to many parameters, including loading and unloading speed, applied load, anisotropy of materials, defects in the sample, method of measurement, temperature, etc⁴. As a consequence, a great number of values on hardness are reported for each crystal, which makes selecting the reliable hardness value of a material a great challenge.

Though hardness tests are easy to conduct, interpretation on hardness values are complex. Usually, the experimental hardness value is found decreasing with increasing load, which is referred to as the size effect²¹. When the load reaches a certain level, the measured hardness value will not decrease anymore. This asymptotic value in the hardness-load curve is commonly recommended as the reliable hardness value of a hard and brittle material²². One question arising from the size effect is: will the hardness monotonously increases with the decrease of load and approaches infinite? Despite of the plateau at large loads in the hardness-load curve, another plateau was obtained at small loads during hardness measurements by Wang et al.²³. As illustrated in their work, the asymptotic value associated with small loads is more or less a constant, while the asymptotic hardness value in accordance with large loads depends strongly on microstructures²³. In addition, transition from the constant value to the trend of decreasing with increasing load was found coincident with the onset of cracking around the indentation²⁴. It means that the constant hardness level obtained at small loads is probably the “intrinsic hardness” of a material, which is a measure of the resistance to plastic deformations without initiation of any micro-cracks. In contrast, the asymptotic hardness level at large loads is a complex composite of the resistance to plastic deformation and fracture with the microstructure saturated by micro-cracks. For simplification, the asymptotic hardness level at large loads is called “engineering hardness”.

While microscopic hardness models assume a perfect crystal, the predict hardness value should be close to the “intrinsic hardness” of a material, since only when the material is lightly deformed that the material can still be well characterized as a continuous crystal. Therefore, the predicted hardness values should be compared to hardness values measured at small loads instead of the asymptotic hardness level at large loads. There is no doubt that the “engineering hardness” of a material is a crucial property in its practical uses due to the severe service environment. Even though a material with high “intrinsic hardness” may display low “engineering hardness”, it is essential that a material with high “engineering hardness” should at least contain some components with high “intrinsic hardness”. Therefore, it is desirable to develop microscopic models to explore potential intrinsic hard materials.

In the present work, three major modifications were introduced to the theoretical hardness model proposed by Gao¹⁰. After modification, the model predicts well the intrinsic harness values of many crystals, including those crystals composed of complex chemical bonds. The modifications are:

(1) The valence electron of a chemical bond should be specified based on its bonding nature instead of equally partitioning of valence electrons of an atom to its connecting bonds;

(2) The bond volume v^μ is not only proportional to the cubic power of bond length (d^μ)³, but also proportional to its valence electron number n^μ;

(3) Deriving the hardness of a crystal from the hardness values of chemical bonds should be averaged based on valence electron population.

All these modifications emphasize the crucial role of valence electron populations, which means that properties of a crystal depend strongly on bonding states of valence electrons. Both decomposing the crystal into pseudo-binary crystals and estimating properties of the crystal from properties of pseudo-binary crystals should be based on valence electron populations of different bonding states. The model becomes self-consistent by introducing these modifications, which is equivalent to decomposing a crystal to a set of well-defined chemical bonds with 2 valence electrons. The fundamental idea of these modifications may also be applicable to other models associated with chemical bonds, e.g. models to estimate thermal expansion²⁵ or bulk modulus^26,27. In general, derivations of these models usually start from simple crystals comprised of well-defined chemical bonds, such as crystals in Table 1. Exploring a self-consistent way to define equivalent chemical bonds may directly extend the models suitable for complex crystals.

How to cite this article: Dai, F.-Z. and Zhou, Y. A Modified Theoretical Model of Intrinsic Hardness of Crystalline Solids. Sci. Rep. 6, 33085; doi: 10.1038/srep33085 (2016).