Sphalerite and wurtzite structured materials are the important members in the family of superhard materials, among which diamond and lonsdaleite are the most prominent representatives. Diamond is widely used in industries because of their excellent mechanical properties [[1], [2], [3], [4], [5], [6], [7], [8], [9]]. With the increasingly stringent service conditions, the demand for materials with excellent strength, toughness and thermal stability becomes more and more urgent. Lonsdaleite may be a latent promising candidate that could meet the demand. Lonsdaleite, as a hexagonal polymorph of diamond whose hardness exceeds that of diamond, is usually called hexagonal diamond because of its crystal structure [10 1¯ 2]. It has been found that the impact diamonds and meteorites containing lonsdaleite possesses mechanical properties superior to diamond [10,11,[13], [14], [15]]. Theoretical predictions indicated that the Vickers hardness of lonsdaleite reaches 152 GPa, which is over 50% harder than that of diamond [16]. This phenomenon is of significant importance because the strong covalent bond solids involved are indispensable to fundamental scientific research and practical applications in many fields [[17], [18], [19]].

At present, lonsdaleite has never been obtained as an independent single crystal [11,13,[20], [21], [22]]. Nemeth et al. [13] pointed out that defects in cubic diamond provide an explanation for the characteristic d-spacings and reflections reported for lonsdaleite. Salzmann et al. [21] provide a quantitative analysis of cubic and hexagonal stacking in diamond samples by X-ray and DIFFaX software package. They believe that the volume fraction of lonsdaleite in the dual-phase or multi-phase solids is less than 60%. Recently, Murri et al. [22] presented MCDIFFaX analyses of X-ray data, combined with experimental and theoretically calculated Raman spectra, and HRTEM analyses of suites of impact diamonds from Popigai to establish the nature and extent of hexagonal stacking within the structures. It is worth mentioning that the mechanical properties of cubic/hexagonal biphase are generally better. Baek et al. [23] studied the impact diamonds rich of lonsdaleite, and suggested several biphasic structures with unique mechanical properties. Tanigaki et al. [14] found that the nano-polycrystalline diamonds rich of hexagonal structure have higher stiffness than single crystal diamond. They confirmed that the stiffness arises from the increase of hexagonal structure by calculating the elastic constants of the models with different volume fractions of hexagonal structure. Chen et al. [24] observed recently a stable two-phase cBN/wBN directly using a transmission electron microscopy (TEM). Although the studies mentioned exhibit the possible structure and unique mechanical properties of diamond/lonsdaleite (or cBN/wBN) biphases, the response and structural reorganization of diamond/lonsdaleite (or cBN/wBN) biphases under severe thermal-mechanical loading conditions have not yet been studied, which should be clarified to provide available information for the design of such kind of high-performance materials.

In this paper, we validate several possible diamond/lonsdaleite biphases and cBN/wBN biphases by investigating their thermodynamic behaviors with first-principles calculations. The stress-strain (σ-ε) and total energy-strain (E-ε) curves of the biphases under various loading conditions are presented, which could provide available information on the structural reorganizations of these biphasic systems during deformation. We also try to explain the unique deformation behavior of the biphasic systems with minimum elastic energy principle.

Our ab initio study was performed using the Vienna ab initio simulation package (VASP) code [25], in which the projector augmented wave (PAW) method [26] is adopted, and the Perdew-Burke-Ernzerhof (PBE) version [27] of the generalized gradient approximation (GGA) was employed to describe the exchange-correlation functional. The electron-ion core interactions of carbon, boron and nitrogen atoms are modeled using pseudo-potential with the electronic configurations 2s²2p², 2s²2p¹ and 2s²2p³, respectively. The plane wave with an energy cutoff of 500 eV is adopted to describe the electronic wave functions and a regular grid with k-point 5 × 9 × 3 is used for the integration in reciprocal space. The convergence tests show that it is suitable for the energy convergence criteria for electronic iterations and Hellmann-Feynman force per atom to be set as 1 × 10^-4 eV and 0.001 eV/Å, respectively. Periodic boundary conditions are applied in the all the three directions of the supercell used. The basic parameters calculated in this work are listed in Table 1, which are in agreement with those from other experimental [12,28] and theoretical [29] work. To consider the thermal stability of the material, we perform standard ab initio MD simulations, and the detailed parameter settings are referred to our previous work [7].

Nanocrystalline materials show many unique mechanical properties compared with conventional crystalline materials [[30], [31], [32], [33], [34], [35]]. However, experiments can hardly monitor dynamically the internal structures of materials during deformation, so it is difficult to clarify the deformation mechanism and the evolution of material properties. The first-principles calculation code adopted contains reliable potential functions for more realistic analysis of the material responses, and could capture the deformation behavior and the corresponding internal structures simultaneously, which could provide available information for the design of material structures. In addition, it should be noted that the strength obtained with first principles calculation should be the upper limit of the material strength. With the maturity of synthetic technology, the strength of actual materials would gradually approach the theoretical value. Such as, the strength of a nanoscale single-crystal diamond under large elastic deformation can reach 89-98 GPa [36], close to the theoretical limit. Recently, it was found [37] that the nanoneedle with <001> crystal direction diameter of 60 nm showed high elastic tensile strain (13.4%) and tensile strength (125 GPa), which are close to the theoretical limit predicted by first-principles calculation. All the above results show the significance of the prediction with first principles calculation.

The ideal strength is defined as the stress at the peak of the σ-ε curve obtained by first principles calculation. In this work, the tensile stress-strain curve is calculated by applying incremental strain to a periodic supercell along the prescribed direction, and the stress can be obtained from the DFT stress tensors in VASP calculations.

The relationships between σzx and σzz in the uniaxial and biaxial are given as Eq. (1) [9,38]:

where C is a constant representing the compressive stress applied in the z direction, and C = 0 corresponds to pure shear; ø = 68°is the centerline-to-face angle of the indenter [36] used in Vickers Hardness testing. In this study, we simulate two kinds of biaxial loadings, i.e., σ_zz=σ_zxtan68° and σ_zz=200GPa, respectively.

Figs. S1(a) and (b) in Supporting Information show the changes of the model during the calculations of the uniaxial compressive and pure shear responses, respectively, where the black denotes the initial models, and the pink denotes the relaxed deformed ones. The uniaxial compressive strain is defined as $ε=\frac{l_{0}-l_{ z } }{ l_{0}} $, and pure shear strain is defined as ε=tanα. Under uniaxial compressive loading conditions, a compressive strain is applied along the zz direction, as shown in Fig. S1(a), and all atomic basis vector and all atoms in the model are relaxed until the other five components of the Hallmann-Feynman stress tensor are less than 0.1 GPa. Then another stress component is the stress state under uniaxial compression. Similarly, apply pure shear strain along the xz direction, as shown in Fig. S1(b), to obtain pure shear stress. Fig. S1(c) shows the change of the model subjected to biaxial compressive-shear straining with the ratio determined by Eq. (1), during the deformation the model is relaxed after each increment of loading until the other four components of the Hallmann-Feynman stress tensor are less than 0.1 GPa.

The response of materials simulated in this work can be divided into elastic and plastic stages. The strongly covalent solids usually fail by brittle fracture without obvious plastic characteristic. In elastic stage, σ-ε relation can be expressed as σ_i = C_ij ε_j, where C_ij represents the elastic matrix [39,40]. The independent elastic constants for cubic solids are C₁₁, C₁₂ and C₄₄, with which the elastic matrix can be expressed as Eq. (2):

The bulk moduli B₀, Young’s moduli E, and shear moduli G of diamond and cBN can be calculated with the following relationships (Eqs. (3) and (4)) [41]:

where E_ijk is the Young’s modulus in the [ijk] direction, and S_ij the compliance tensor ([S]=[C]^-1); l_i1, l_j2, and l_k3 are the direction cosines of the direction [ijk]. The Young’s moduli along the three low index directions [100], [110] and [111] can be expressed as Eq. (5):

The expression for shear moduli G is a little more complicated, because it involves the direction and plane of shear. The shear moduli along any directions on the {100} planes are the same, so the shear moduli of the cubic crystal can be expressed as G₁₀₀ = C₄₄; but the shear modulus varies along different directions on the {110} and {111} planes [41]. The shear modulus along the (111)<110> direction is given as Eq. (6):

The elastic constants obtained are summarized in Table I, which coincide with those reported previously [4].

In order to verify the reliability of the σ-ε relations, we selected the elastic stages of the tensile and compressive curves of diamond along the <111> direction, as shown in Figs. S2(a) and (b) in Supporting Information. The effective Young’s moduli along an arbitrary direction [ijk] can be calculated with the corresponding strain and the stress, σ_ijk=E_ijkε_ijk. The initial slopes of the σ-ε curves indicate that the Young’s moduli are 1137 GPa and 1198 GPa for tension and compression, respectively, which agree well with the corresponding moduli of 1166 GPa calculated with single-crystal elastic constants. The validation of shear moduli on the planes (001) and (110) are also consistent with those calculated with single-crystal elastic constants, as shown in Fig. S2(c) and (d) in Supporting Information. The shear moduli are 557 GPa and 490 GPa for (001) plane and (111)<110> direction, respectively, which agree well with the corresponding moduli of 555 GPa and 491 GPa calculated with single-crystal elastic constants.

The hexagonal system has six-fold rotation symmetry around [0001]. The independent elastic constants in the hexagonal phase are C₁₁, C₁₂, C₁₃, C₃₃ and C₄₄, and the elastic matrix can be expressed as Eq. (7):

For a material with hexagonal symmetry, the bulk modulus B₀ can be expressed as Eq. (8) [42,43]:

where $\Delta ={{C}_{33}}\left( {{C}_{11}}+{{C}_{12}} \right)-2C_{13}^{2}$.

The Young’s moduli, E, is dependent on crystallographic direction, and E is isotropic in the basal plane, but varies if one moves toward the c axis. The Young’s moduli can be calculated with Eq. (9) [41]:

where l₃ is the cosine of the angle with the basal plane. Both the directions <1210> and <1010> are in the basal plane, where l₃ is equal to 1. The direction <0001> is perpendicular to the base plane, where l₃ is equal to zero. The obtained Young’s moduli E_p (parallel with the basal plane) and E_v (perpendicular to the basal plane) are calculated as Eq. (10):

The calculated E_p and E_v of lonsdaleite and wBN are summarized in Table 1.

Understanding each single-phase structure is the basis for an in-depth analysis of the properties of a biphasic structure. Fig. 1(a) and (b) show the crystallographic stacking sequences in diamond and lonsdaleite, respectively, where the stacking sequence of diamond can be expressed as AαBβCγAαBβCγ… in the [111] direction, which are of typical cubic crystal; the stacking sequence of lonsdaleite and wBN can be expressed as Aαα′A′Aαα′A′… along the [0001] direction, where A' denotes π-rotated A plane [43]. For cubic crystal, the (111)<112> slip system contains the (111)[11 $\bar{2}$]easy shear direction (ESD) and (111)[ $\bar{11}$2] hard shear direction (HSD) [42]. Figs. S3 and S4 in Supporting Information show the stress-strain curves of these four single phases subjected to tensile, compressive and shear strain, respectively, which are in agreement with the existing results [9,16,44,45]. In addition, an important structural feature is the presence of two non-equivalent slip planes, i.e., the (111) plane for diamond and cBN, and the (0001) plane for lonsdaleite and wBN [[46], [47], [48]]. As shown in Fig. 1, there are two nonequivalent candidate (111) planes: a narrowly spaced atomic layer and the widely spaced atomic layer, which are called the glide-set (G) and the shuffle-set (S), respectively.

It has been found in experiments that diamond/lonsdaleite biphasic structure could possess excellent thermal-mechanical properties [13,23]. Fig. 2(e) shows the STEM image from the Canyon Diablo sample, and the existence of hexagonal lonsdaleite is observed, and the schematic diagram of the model is further shown in Fig. 2(f). Thus, we build a set of possible diamond/lonsdaleite biphasic structures, as shown in Fig. 4(a-d), respectively. The lattice mismatch between (111) plane of diamond and (0001) plane of lonsdaleite is 0.03%, which is tiny [23]. In the four diamond/lonsdaleite biphasic models shown in Fig. 2(a-d), the layer with the [111] orientation of diamond and the layer with the [0001] orientation of lonsdaleite are arranged alternately, with the white part denoting diamond while the yellow part denoting lonsdaleite. The four samples are visually named as 4L + D, (3 + 1)L + D, (2 + 1 + 1)L + D and (1 + 1+1 + 1)L + D, respectively. On the other hand, it is observed recently with a TEM [22] that there is a stable dual-phase BN consists of cBN and wBN, as shown in Fig. 2(g) and (h).Therefore, four cBN/wBN biphasic models are also built similarly, as shown in Figs. S6 (c-f) in Supporting Information, named as 4w + c, (3 + 1)w + c, (2 + 1 + 1)w + c and (1 + 1+1 + 1)w + c, respectively. We compare the energies of the diamond/lonsdaleite biphasic models and their single crystal counterparts, and show the results in Fig. S4(a) in Supporting Information, where it can be seen that the energy per atom of each diamond/lonsdaleite biphasic model lies between those of diamond and lonsdaleite, implying the rationality of the existence of the biphases. Similar results can also be found for the cBN/wBN biphases.

To further demonstrate the stable existence of the biphases, we compare qualitatively the thermal stability of the diamond/lonsdaleite biphasic models and their single crystal counterparts using ab-initio MD simulations. Fig. 4(a) and (b) show some snapshots of the simulations, which depict the structural changes of diamond and lonsdaleite at T = 1200 K. It should be noted that, before selecting this temperature, we performed simulations for the response of each model in Fig. 3, at various temperatures between 0-1500 K, and found that the structure of single crystal diamond would be basically destroyed as T = 1500 K. At 1200 K, both single crystal diamond and diamond /lonsdaleite structures could basically keep intact, while pure lonsdaleite would be destroyed completely. Therefore, for carbon materials, 1200 K should be an appropriate reference temperature. It can be seen in Fig. 4 that at T = 1200 K the cubic structure of diamond could basically be maintained, while the structure of lonsdaleite has been destroyed completely. More detailed information is shown in Supplementary Movies 1 and 2 obtained with ab-initio MD simulations. To further explore why no pure lonsdaleite could be found in experiment, we simulate the thermodynamic behavior of pure lonsdaleite at room temperature with ab-initio MD, and discuss the results (Section C of Supplementary Information), which could explain why pure lonsdaleite is difficult to exist naturally. In contrast, the thermal stability of the diamond/lonsdaleite biphases is shown in Fig. 4(c-f), where it can be seen that the structural changes are negligible compared to the initial structure, indicating that diamond/lonsdaleite biphases can stably exist, which is consistent with the conclusions from experiments [13,14,23,49]. More detailed description about the ab-initio MD simulations can be found in Supplementary Movies 3-6.

The variations of the thermal stability of cBN/wBN biphases are similar to that of diamond/lonsdaleite biphases. Fig. S6 shows some snapshots from the ab initio MD simulations for the structural changes in cBN/wBN biphases and their single crystal counterparts at T = 2000 K. So far, we have successfully verified the rationality of the existence of the biphasic structure using DFT calculations, which is consistent with experimental observations. In the following subsection the stress and energy responses of the biphasic structure under large shear strains will be investigated.

It was found recently in experiments that the mechanical properties of diamond/lonsdaleite biphases are superior to diamond [12,14,23]. Similar conclusion could also be obtained for cBN/wBN biphases [50]. To explore the strengthening mechanism in such biphasic structures, we investigate the responses of these biphasic models subjected to pure shear strain. Fig. 4(a) shows the comparison between the shear stress-strain (σ-ε) curves of diamond/lonsdaleite biphases and those of diamond under pure shear along ESD and along HSD. It can be seen that the σ-ε curves of biphasic models show a distinct zigzag except the (3 + 1)D + L model. We first select 4L + D to analyze this distinct zigzag phenomenon. The σ-ε curve of 4L + D can be divided into four stages, as shown in Fig. 4(a). Stage I: σ increases from ε₀ = 0 to ε₁ = 0.24 when the first stress peak of 105.7 GPa appears, which is higher than the peak stress of 90 GPa of diamond subjected to pure shear in ESD. The change in the atomic structure in Stage I can be seen in Fig. 4(c) and (d), where the structure remains biphasic. Stage II: the σ-ε curve of 4L + D drops sharply from the peak at ε₁ = 0.24 to the valley at ε₂ = 0.25. Fig. 4(d) and (e) show the atomic structures at these two strain states, where the partial slips occur in the lonsdaleite region, leading to the phase transition from hexagonal to cubic. Phase transition takes place by the break and recombination of the covalent bonds on the glide-set plane, i.e., as ε increases from 0.24 to 0.25, Bonds C₁-C₂ and C₁-C₃ break, followed by the formation of Bonds C₁-C₃ and C₁-C₄. To clarify this process, we calculate the generalized stacking fault energy (GSFE) curves for (111)<112> slip systems on the glide-set and on the shuffle-set planes, as shown in Fig. 5(a), where the energy barrier (γ_Ug) on the glide-set plane is 5.58 J/m², which is smaller than that on the shuffle-set plane (8.58 J/m²), indicating that the slip should be dominated by the partial slip on the glide-set plane, which coincides with the conclusion that the bond break and recombination on the glide-set plane. Fig. 5(b-d) show some key atomic configurations, and it is obvious that the cubic configuration has lower total energy. Stage III spans from ε₃ = 0.25 to ε₄ = 0.40, when the cubic structure [Fig. 4(e)] continues to be sheared and the stress rises quickly to the second peak of 130 GPa, which is close to the stress peak (134.3 GPa) of diamond subjected to shearing in HSD and about 1.23 times the stress of the first peak. The rapid increase of stress in this stage can be attributed to the slip along HSD of diamond. Stage IV spans from ε₄ = 0.40 to ε₅ = 0.41 when the stress falls quickly due to the graphitization of the structure, as shown in Fig. S8(f) and (g). The Supplementary Movie 7 describes in details the structural changes during the shear deformation, where the transition from biphasic phase to cubic phase can clearly be observed.

It is worth to discuss in more details the strengthening mechanism related to the phase transition in the biphasic structure. During the loading process, the material always deforms along the slip system of the lowest resistance. For diamond, the slip should occur along the ESD instead of HSD, corresponding to the maximum shear stress of 90 GPa instead of 130 GPa. In the biphasic structure, initial slip should always occur in the ESD in the lonsdaleite region, inducing the phase transition; thereafter, further slip would occur in the HSD, leading to the marked increase in the maximum shear stress. It is noteworthy that this strengthening mechanism could not only improve the material strength but also increase its ductility.

As shown in Fig. 4(a), the extreme stresses of (1 + 1+1 + 1)L + D and (2 + 1 + 1)L + D are close to the peak stress (134.3 GPa) of diamond subjected to shear straining along HSD. They also follow the strengthening mechanism by phase transition. The strengthening mechanism of cBN/wBN biphases is similar, where the stresses at the first peaks are close to the maximum stress of cBN sheared along the ESD, while the maximum stresses are close to the peak stress of the cBN sheared along the HSD, as shown in Fig. 4(b).

It should be noted that (3 + 1)L + D does not experience phase transition, therefore, the σ-ε curve does not show a zigzag shape, and the maximum stress of 98.8 GPa is lower than that of other biphasic structures, as shown in Fig. 4(a). This phenomenon can be ascribed to that pure shear stress could not confine the deformation of the model in z-direction, resulting in the failure of the structure in z-direction at large strain. To further verify this speculation, we calculate the responses of the (3 + 1)L + D model subjected to pure shear and to combined shear-compressive deformations (Eq. (1)), respectively, and the results obtained are shown in Section D of Supplementary Information. The biaxial stress state restricts z-direction deformation of the (3 + 1)L + D, so the release of stress is achieved by partial slip on the glide-set plane in the lonsdaleite region. Therefore, under the state of biaxial stress, the response of the (3 + 1)L + D model also follows the strengthening mechanism by phase transition.

We can also interpret the rationality of the strengthening mechanism related to phase transition from the viewpoint of energy. As shown in Fig. 6, we extract the energy of the structure in each loading step, and we can find that the energy-strain (E-ε) curves have the same trend as the σ-ε curves, indicating that the structure has lower energy after phase transition so that it could exist more stably and have the ability to withstand further deformation.

Biaxial compressive-shear loading was suggested recently for simulating the applied stress state during indentation, which simplifies the three-dimensional stress state of the indentation into a biaxial plane stress state containing shear stress (σ_zx) and normal compressive stress (σ_zz) components [38,51]. We first apply a biaxial state of stress containing varying shear stress σ_zx to the diamond/lonsdaleite biphasic models under fixed σ_zz = 200 GPa, and show the shear stress-strain curves (σ-ε curves) in Fig. 7(a). It can be seen in Fig. 7(a) that the σ-ε curves of all the biphasic models exhibit distinct zigzag and the maximum shear stresses are all close to the peak stress (200.7 GPa) of diamond subjected to biaxial stress with σ_zz = 200 GPa along the HSD. It should be noted that the value of σ_zx is increased because of the application of σ_zz that could act as the constraint on the deformation in z direction. To further understand the impact of the constraint, we use the (3 + 1)L + D model as an example, which does not experience phase transition under pure shear. The σ-ε curve of (3 + 1)L + D is shown in Fig. 7(a), where two significant drops appear between ε₁ = 0.22 and ε₂ = 0.23 and between ε₃ = 0.30 and ε₄ = 0.31. The key structural snapshots of (3 + 1)L + D at ε₀ = 0, ε₁ = 0.22, ε₂ = 0.23, ε₃ = 0.30 and ε₄ = 0.31 are shown in Fig. 7(c-g), respectively, which illustrates that the sudden drop in stress is due to the partial slips that take place in the lonsdaleite region. The partial slip eventually leads to phase transition. Thereafter, further slip would occur in the HSD of diamond, as shown in Fig. 7(g). The maximum stress of cBN/wBN biphases is also close to that of cBN along the HSD, as shown in Fig. 7(b). The root cause of the stress enhancement can also be attributed to phase transition.

We also calculate the responses of these biphasic models under Vicker shear strain loading, i.e., σ_zz=σ_zxtan68° is applied to the biphasic models subjected to shear strain, ε_zx. The shear stress-strain curves and the corresponding key structural snapshots are shown in Fig. 8, where for all the biphasic structures, the enhancements of both strength and ductility can be attributed to phase transition. We extract the key structural snapshots of 4L + D subjected to Vicker shear strain to gain an insight into the enhancement mechanisms, as show in Fig. 8(c-f). The transition from biphasic structure to cubic one is caused by the partial slip on the glide-set planes in the lonsdaleite region, which suggests that the phase transition induced strengthening is also applicable to the biphasic structure subjected to biaxial loading.

Compared with pure shear straining, the application of σ_zz could restrict the deformation in z-direction, resulting in the increase of σ_zx. The σ_zx of several key systems under uniaxial and biaxial stress states are summarized and listed in Table 2. The existence of σ_zz does contribute to the enhancement of σ_zx. In addition, biaxial stress state is successfully explored for simulated indentation stress distribution, although it could not fully represent the indentation stress state, it should be closer to the case of indentation than pure shear.

The mechanical and thermal properties of diamond/lonsdaleite biphasic structure were studied using first principles calculations. The atomic reconfiguration by partial slip increases both the strength and strain range, which challenges the traditional concept that the strength and ductility of a material can hardly be enhanced simultaneously. The atomic stacking sequence of the biphasic structure differs from that of cubic diamond, resulting in the appearance of weak atomic layers with low stacking fault energy. Partial slip would occur on the weak atomic layers of the lonsdaleite region when the stress reaches a critical value, which would transform the model structure from biphasic to cubic. The further shear strain would induce the slip along the hard shear direction of the cubic structure, accounting for the extremely high strength of biphasic structure. Meanwhile, ab initio molecular dynamics simulations reveal a stabilization mechanism that a biphasic structure would contribute to the stability of the material, which confirms the conclusion from experiments. These results could not only validate the rationality of biphasic structure repeatedly observed in experiments, but also provide the strengthening mechanism related to phase transition, leading to an explanation for the biphasic structure with unique mechanical properties.

This work was financially supported by the National Natural Science Foundation of China (Nos. 11932004 and 11802045), the National Postdoctoral Program for Innovative Talents (No. BX20190039), the Postdoctoral Program for Innovative Talents of Chongqing (No. CQBX201804), and the Natural Science Foundation of Chongqing (No. cstc2019jcyj-bshX0029). First principles calculations were carried out at Supercomputing Center of Lv Liang Cloud Computing Center in China.

Supplementary material related to this article can be found, inthe online version, at doi:https://doi.org/10.1016/j.jmst.2020.03.005.