A Modified Theoretical Model of Intrinsic Hardness of Crystalline Solids

Super-hard materials have been extensively investigated due to their practical importance in numerous industrial applications. To stimulate the design and exploration of new super-hard materials, microscopic models that elucidate the fundamental factors controlling hardness are desirable. The present work modified the theoretical model of intrinsic hardness proposed by Gao. In the modification, we emphasize the critical role of appropriately decomposing a crystal to pseudo-binary crystals, which should be carried out based on the valence electron population of each bond. After modification, the model becomes self-consistent and predicts well the hardness values of many crystals, including crystals composed of complex chemical bonds. The modified model provides fundamental insights into the nature of hardness, which can facilitate the quest for intrinsic super-hard materials.


Theoretical Model and Modifications
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.
Scientific RepoRts | 6:33085 | DOI: 10.1038/srep33085 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 9 : 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 10 : where A is a proportional coefficient, N a is the covalent bond number per unit area, P μ is the Mulliken overlap population 13 , and v μ is the bond volume. In the equation, P μ /v μ characterizes the strength of the bond, as suggested by Gao 10 . 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: 2/3 5 /3 However, in the original work by Gao 10 , H μ was expressed as: 5/3 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 14 : 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 2 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 15 . 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.
The second modification is made on the definition of bond volume v μ . The original definition of bond volume is introduced by Levine 16 in 1973, which assumes that the volume of a bond is proportional to (d μ ) 3 . 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: 3 3 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: n N n N 1/ 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 10 , 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.

Evaluation of the Modified Model
Crystals from refs 7,9-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 17 . The Vanderbilt-type ultrasoft pseudopotential 18 and exchange-correlation described by generalized gradient approximation 19 were employed. The plane wave cutoff energy was set to be 500 eV. k-points mesh with a separation of 0.03 Å −1 according to Monkhorst-Pack method 20 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 2 = 0.984, which is close to the value of 740 suggested by Gao 10 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.
To further verify the capability of the modified model, it is applied to investigate the hardness of TMB 2 s (TM = Ti, Zr, Hf, Re and Os). TiB 2 , ZrB 2 and HfB 2 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 2 has a simple hexagonal structure (space group P6 3 /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 2 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 2 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 2 for example. According to P μ , there are four types of bonds in ReB 2 ( 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 2 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 2 (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 2 s, as shown in Table 4.
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   Figure 4 also reveals that the predicted hardness values for those TMB 2 s are also significantly improved after modification.
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 4 . 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 21 . 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 22 . 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. 23 . 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 23 . 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 24 . 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.

Conclusions
In the present work, three major modifications were introduced to the theoretical hardness model proposed by Gao 10 . 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 μ ) 3 , 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 25 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.