Mechanical Properties and Stability of Body-Centered-Tetragonal C₈ at High Pressures

1 Introduction

Superhard materials have always attracted the attention of researchers because of their high compression strength, chemical inertia and high hardness. Carbon materials, as one of the most important components of superhard materials, have attracted much attention of the researchers. Carbon materials have various allotropes such as graphite, diamond, carbon nanotubes and fullerenes [1], [2], [3]. The reason why carbon can form many allotropes is because the carbon atoms can form sp-, sp²- and sp³-hybridised bonds [4]. The strong covalent sp³-hybridised bonds have always been an important factor in its ability to form superhard materials. Researchers have proposed a number of superhard carbon allotropes [5], [6], [7], [8]. These materials can remain stable under high pressure and exhibit unique properties. Recently, the body-centered-tetragonal C₈ (Bct-C₈) was investigated using the first-principles calculations [9]. Bct-C₈ can be regarded as one of the possible products of compressing carbon nanotube bundles [10], and it has a sp²-hybridised bonding character. In addition, the experiment proves that Bct-C₈ is a semiconductor. It is known that for semiconductors, tailoring electronic properties of semiconductor nanostructures is critical for its applications. Strain has long been used to tune the electronic properties of semiconductors [11], [12]. It was found that nanostructures maintain their integrity under a much higher strain than their bulk counterparts [13], [14], which dramatically expands the strength of the applicable strain to nanostructures [15], [16], [17]. In this paper, we calculated the band structure and density of states (DOS) of Bct-C₈ at 0 GPa. We also calculated the elastic constants at a pressure range from 0 to 100 GPa to ensure the mechanical stabilities of Bct-C₈. The phonon spectrum under 100 GPa is also calculated to ensure the dynamic stabilities at high pressures. Based on the elastic constants, the average acoustic velocity, Debye temperature and elastic anisotropy are also calculated. In addition, the ideal stress-strain relations are also studied.

2 Theoretical Method and Computational Details

Theoretical calculations were carried out using the first-principles calculations based on the density functional theory (DFT) [18], the Perdew–Burke–Ernzerhof (PBE) [19] exchange-correlation functional and projector-augmented wave [20] potentials were employed for self-consistent total energy calculations and geometry optimisation. The calculations were performed using the Vienna Ab-initio Simulation Package (VASP) [21]. The kinetic energy cutoff for the plane wave basis set was chosen to be 900 eV. The k-point samplings with 8 × 8 × 20 in the Brillouin zone were performed using the Monkhorst–Pack scheme. The energy convergence criteria for electronic and ionic iterations were set to be 10⁻⁵ and 10⁻⁴ eV, respectively. With this parameter settings, the calculated energy difference is less than 1 meV in the total energy per atom. Bulk modulus, shear modulus, Young’s modulus and Poisson’s ratio were estimated by using the Voigt–Reuss–Hill approximation [22].

3 Results and Discussion

The optimised geometric structure of Bct-C₈ with the space group I₄/mmm of tetragonal symmetry is shown in Figure 1. The Wyckoff position is 16l (0.288, 0.884, 0), and the lattice constants are a = b = 6.57 Å and c = 2.52 Å. One can see that the calculated lattice parameters are in an excellent agreement with the previous results, indicating that our calculations were reliable. As shown in Figure 1a, there are three kinds of C–C bonds in the crystal structure of Bct-C₈, they are marked as d1, d2 and d3, and the values are 1.52 Å, 1.60 Å and 1.55 Å, respectively. One can see from Figure 1b–d that the structure of Bct-C₈ is composed of a number of carbon nanotubes interlinked longitudinally, and some large holes are evenly distributed in the structure.

Figure 1:

The schematic drawing of the crystal structure of Bct-C₈: (a) the perfective view; (b) the top view; (c) the side view; and (d) the dash lines indicate the unit cell.

Figure 2:

(a) The band structure of body-centered-tetragonal C₈ (Bct-C₈) at 0 GPa, the black curves represent the results calculated with Perdew–Burke–Ernzerhof (PBE) functional, and the red curves represent the results by hybrid Heyd–Scuseria–Ernzerhof (HSE06) functional; (b) the density of states (DOS) of Bct-C₈ at 0 GPa, calculated by PBE functional.

The band structure and the DOS of Bct-C₈ at 0 GPa are calculated to study its fundamental physical and chemical properties. The band structure is shown in Figure 2a, in which the black dashed line represents the Fermi level (E_F). The red dot and green dot [solid dot for PBE functional and the hollow dot for hybrid Heyd–Scuseria–Ernzerhof functional (HSE06) functional] represent the conduction band minimum (CBM) and the valence band maximum (VBM), respectively. The black curve represents the results we calculated using the PBE functional. The CBM locates at Z point and the VBM locates at Γ point. The band gap is 0.44 eV. It is clear that Bct-C₈ is an indirect semiconductor. We know that the band gap calculated by DFT usually should be smaller than the actual values, thus we use the more precise HSE06 functional to correct the band gap of Bct-C₈. As shown in Figure 2a, the red curve represents the results by HSE06 functional. It can be seen that, using both functionals, similar band structures were obtained. The VBM and CBM calculated by HSE06 functional locate at the same positions as those calculated by PBE functional, and the band gap of Bct-C₈ calculated with HSE06 functional increases to 1.37 eV.

The DOS of Bct-C₈ is shown in Figure 2b. We can see that the valence band region can be divided into two parts, in the first part (−22 to −12 eV), the total DOS is mainly contributed by C-s states, and the second part (−12 to 0 eV) is mainly due to the C-p states. In general, for the valence band, the contribution of C-s to the total DOS is getting smaller, and the contribution of C-p is gradually increasing. In addition, the conduction band region is mainly characterised by the C-p states, and the DOS near Fermi level is mainly originated from the C-p orbital electrons.

To study the compressibility of Bct-C₈ under pressure, we calculated the normalised volume V/V₀ ratios (where V₀ is the equilibrium volume at 0 GPa) as the function of pressure, and compared with diamond, c-BN and Imma-C. As Figure 3 shows, from 0 to 100 GPa, the diamond is the most incompressible; Bct-C₈ is more compressible than c-BN and Imma-C, when the pressure increases to 100 GPa, V/V₀ ratio of Bct-C₈ is 0.82, and V/V₀ ratios of c-BN, Imma-C, and diamond are 0.83, 0.84 and 0.85, respectively.

Figure 3:

The V/V₀ of Bct-C₈ compared with diamond, c-BN and Imma-C at pressures.

Elastic properties have always been a very important part of the many properties of materials. From the elastic constant, we can obtain the stability, hardness and anisotropy information regarding the materials [23]. It is known that, for different crystal systems, there are different mechanical stability criteria. For tetragonal symmetry, the criteria of the mechanical stability are given by [24]:

At high pressures, the criteria of mechanical stability are:

where C~11=C11−P,C~33=C33−P,C~44=C44−P,C~66=C66−P,C~12=C12+P,C~13=C13+P, and P is the hydrostatic pressure.

The elastic constants of Bct-C₈ at 0–100 GPa satisfy all the stability criteria above, indicating that the Bct-C₈ is mechanically stable up to 100 GPa. The obtained elastic constants at different pressures are presented in Figure 4a. It is found that the values of C₁₁, C₃₃, C₄₄, C₁₂, and C₁₃ increase with different rate with pressures ranging from 0 to 100 GPa, and C₁₁ and C₃₃ increase faster than other elastic constants. As the pressure increases, C₁₁ and C₃₃ increase from 868 to 1448 GPa and from 961 to 1499 GPa, respectively. C₆₆ keeps decreasing from 132 to 61 GPa. Furthermore, the phonon spectra at 100 GPa are also calculated to ensure the dynamical stability. As shown in Figure 4b, the phonon spectra of Bct-C₈ at 100 GPa have no imaginary frequency in the whole Brillouin zone, indicating that Bct-C₈ is dynamically stable up to at least 100 GPa.

Figure 4:

(a) Elastic constants C_ij (GPa), bulk modulus B (GPa), shear modulus G (GPa) and Young’s modulus E (GPa) as a function of pressure; (b) phonon spectra of Bct-C₈ at 100 GPa.

The bulk modulus describes the volumetric elasticity, and the shear modulus describes a material’s tendency to shear when acted upon by opposing forces. By using the Voigt–Reuss–Hill approximations [22], [23], [24], [25], [26]:

and the B and G can be written as:

The changes in B and G with pressures are also shown in Figure 4a. It is found that the value of B is increasing from 327 to 671 GPa with the pressure increase from 0 to 100 GPa, and as the pressure increases, the value of G increases slightly and then decreases slightly.

In material science, ductility is a solid material’s ability to deform under tensile stress. If a material is brittle, when subjected to stress, it will break without significant deformation (strain). Additionally, these material properties are dependent on pressure. The ratio of bulk modulus to shear modulus B/G proposed by Pugh [27] shows a ductile or brittle feature of a material. The high B/G ratio indicates that the material is ductile, and the low value indicates that the material is brittle. Concretely, If B/G > 1.75, the material behaves in a ductile way [28]; otherwise, the material behaves in a brittle way. The B/G ratio as a function of pressure is calculated, and the results are illustrated in Figure 5. From Figure 5, it can be seen that the B/G ratio of Bct-C₈ at 0 GPa is 1.08, which is less than 1.75, thus Bct-C₈ at 0 GPa is brittle. With the pressure increasing from 0 to 100 GPa, the value of B/G ratio increases monotonously. When the pressure increases to 70 GPa, the B/G ratio starts to be greater than 1.75, and the B/G ratio is 2.23 when the pressure reaches 100 GPa. This indicates that when the pressure increases from 0 to 100 GPa, Bct-C₈ changes from being brittle to ductile.

Figure 5:

Calculated B/G ratio and Poisson’s ratio ν of Bct-C₈ at pressures.

To obtain more information about the elastic properties, Young’s modulus E and Poisson’s ratio ν were calculated. As a measure of the stiffness of a solid material, Young’s modulus E is defined as the ratio between stress and strain. When a material receives a tension or compression in unidirectional, the absolute value of the ratio of transverse contraction strain to longitudinal extension strain is called Poisson’s ratio ν. They are given by [22], [29]:

As shown in Figure 4a, at 0 GPa, the Young’s modulus E of Bct-C₈ is 696 GPa, and the Young’s modulus increases with the pressure. When the pressure increases to 70 GPa, E reaches a maximum value of 811 GPa, and then, as pressure increases, the E decreases. As known, the larger the value of E, the stiffer is the material. We can see that at 70 GPa, Bct-C₈ possesses the largest stiffness. From Poisson’s ratio, we can obtain a lot of information about the bonding forces characteristics of a material. The ν reflects the stability of a material against shear deformation, and the material with bigger ν shows better plasticity. The typical value of ν is 0.1 for covalent materials and 0.33 for metallic materials [30]. From Figure 5, we can see that the ν of Bct-C₈ increases from 0.15 to 0.31 with the pressure increase, indicating that with the pressure increase, the directionality degree of covalent bonding becomes weaker and the plasticity becomes better.

There are different physical and chemical properties of the crystal in different directions, which are represented as the crystal anisotropy. Most materials exhibit anisotropic behaviour. Elastic anisotropy research is an important application in material science. In this paper, we mainly discussed the anisotropy of elastic modulus of materials. The 3D figures of the Young’s modulus for all possible directions of Bct-C₈ are calculated and shown in Figure 6. For an isotropic material, the 3D figures are represented as a sphere, and the deviation degree from the sphere represents the anisotropy [31]. At 0 GPa, the maximum of Young’s modulus E_max is 960 GPa, and the minimum of Young’s modulus E_min is 414 GPa. At 100 GPa, the E_max and E_min are 1456 and 232 GPa, respectively. The average value of all directions at 0 and 100 GPa are 745 and 889 GPa, and the ratio E_max/E_min are 2.32 and 6.28, respectively. These results indicate that the anisotropy increases when the pressure increases.

Figure 6:

(a) The direction dependence of Young’s modulus (GPa) of Bct-C₈ at 0 GPa; (b) 100 GPa.

The acoustic velocity is of great significance to the study of the chemical bonding characteristics, and the symmetry and propagation direction of the crystal determine the acoustic velocity. Using the elastic constants, the phase velocities of pure transverse and longitudinal modes of the Bct-C₈ can be calculated by following Brugger’s procedure [32]. In the principal directions, the acoustic velocities in a tetragonal symmetry are given by the following expressions: for [100],

for [010],

for [110],

where v_l is the longitudinal acoustic velocity, v_t1 and v_t2 are the first and the second transverse mode, respectively. The density of the structure is ρ. The acoustic velocities are calculated based on the elastic constants, so the anisotropic properties of the materials will also be reflected in the acoustic velocities. Calculated the average longitudinal acoustic velocity is vlm=(B+4G/3)/ρ , and the average transverse acoustic velocity is vtm=G/ρ . The average acoustic velocity is vm=[(2/vtm3+1/vlm3)/3]−1/3. Based on the calculation results of the average acoustic velocity, we can obtain the Debye temperature:

where N_A is Avogadro’s number, h and k_B are the Planck and Boltzmann constants, respectively, n is the total number of atoms in the formula unit, M is the mean molecular weight, and ρ is the mass density. The calculated acoustic velocities and Debye temperatures of Bct-C₈ are shown in Table 1. The density at 100 GPa is larger than that at 0 GPa, meanwhile, the longitudinal acoustic velocities at 100 GPa along different directions are larger than that at 0 GPa. As the pressure increases, the average acoustic velocity and Debye temperature increase, the average acoustic velocity increases by 48 % from 0 to 100 GPa, and the Debye temperature increases by 59 % from 0 to 100 GPa. Due to the influence of Debye temperature on the strength of covalent bonds in solids, we can conclude that the covalent bong strength of Bct-C₈ will be stronger as the pressure increases.

Fracture toughness K_IC measures the resistance of a material against crack propagation, and is one of the most important mechanical properties of materials [33]. We use the following empirical formula for calculating fracture toughness of Bct-C₈ [34]:

where V₀ is the volume per atom (in m³), G and B are shear and bulk moduli (in MPa), and the unit of K_IC is in MPa ⋅ m^1/2. At 0 GPa, the calculated K_IC of Bct-C₈ is 4.29 MPa ⋅ m^1/2, which is 47 % smaller than diamond (6.33 MPa ⋅ m^1/2), and 39 % smaller than c-BN (5.97 MPa ⋅ m^1/2) [34].

To better characterise the intrinsic hardness of Bct-C₈, the ideal stress-strain curves for large stains were calculated, as seen in Figure 7. For tetragonal phase, the main high symmetry directions, including [001], [100], [110] and [111], are usually considered as the tension directions. The calculated results of ideal tensile strength are shown in Figure 7a. It can be seen that the weakest ideal tensile strength is in the [100] direction with the value of 72 GPa, and the strongest stress response is in the [001] direction with a 93 GPa peak tensile stress. In addition, the responses in the [110] and [111] directions are 73 and 76 GPa, respectively. This means that Bct-C₈ would first cleave in the (100) plane under tensile loadings. To investigate ideal shear strength of Bct-C₈, various shear-sliding planes were systematically studied along different directions under shear deformation, as shown in Figure 7b. It can be seen that the range of shear strength is from 70 to 99 GPa, where the weakest shear strength is along the (100)[010] shear direction with 70 GPa, and the strongest shear strength is in the (110)[001] direction with 99 GPa. That indicates that Bct-C₈ will first cleave in the (110)[001] direction under shear strength. The shear strength in the (001)[100], (100)[010], (010)[101], (110)[001], (110)[1–10], (110)[1–11] and (111)[–1–12] directions are 98, 70, 71, 99, 73, 73 and 84 GPa, respectively. It is clearly seen that the variance between the weakest shear strength of Bct-C₈ and its weakest tensile strength is slight. The value of weakest shear strength is much larger than the superhard criteria (H_V > 40 GPa), which means Bct-C₈ is a potential superhard material.

Figure 7:

(a) Tensile strengths; (b) shear strengths of Bct-C₈.

4 Conclusions

In summary, systematic DFT calculations were performed on Bct-C₈, including the structural stability, elastic anisotropy, electronic properties, tensile and shear strengths. The calculations of elastic constants and phonon spectra show that Bct-C₈ is mechanically and dynamically stable at least up to 100 GPa. The calculated anisotropy of elastic modulus shows that with the pressure increase, the elastic anisotropy increases. The average acoustic velocity and Debye temperature are all calculated based on the elastic constants; the results show that they all increase with the pressure increase. When the pressure increases from 0 to 100 GPa, Bct-C₈ changes from being brittle to being ductile. As the ideal stress-strain curves show, the Bct-C₈ would first cleave in the (100) plane with a 72 GPa peak tensile stress. In addition, the weakest shear strength is along the (100)[010] shear direction with 70 GPa. So we can conclude that the Bct-C₈ is an intrinsic superhard material.