Effects of Vacancies on the Structural, Elastic, Electronic and Thermodynamic Properties of C11_b-VSi₂ by First-Principles Calculations

Abstract

The effects of V and Si vacancies on structural stability, elastic properties, brittleness-toughness transition, Debye temperature and electronic properties of tetragonal C11_b-VSi₂ are investigated using the first-principles calculations. The vacancy formation energy and phonon dispersions confirm that perfect C11_b-VSi₂ and C11_b-VSi₂ with different atomic vacancies are thermodynamically and dynamically stable. The C11_b-VSi₂ with V-atom vacancies is more stable than that with Si-atom vacancies. The introduction of different atomic vacancies enhances the elastic modulus and its anisotropy of C11_b-VSi₂. The electron density difference and densities of state of perfect VSi₂ and VSi₂ with different vacancies are calculated, and the chemical bonding properties of perfect VSi₂ and VSi₂ with vacancies are discussed and analyzed. Additionally, the results show that the chemical bond strength of VSi₂ is enhanced by the introduction of vacancies. Finally, Debye temperatures of perfect VSi₂ and VSi₂ with vacancies are also calculated.

1. Introduction

Transition metal (TM) disilicides have excellent properties such as high strength, high melting point, high temperature creep resistance and good oxidation resistance in high-temperature environments, and have a wide range of promising applications in high temperature devices and semiconductors [1,2,3,4]. Moreover, transition metal disilicide is a promising high-temperature structural material because of its excellent overall performance at high temperatures. In addition, with the rapid development of the aerospace industry, the demand for structural materials with low density and good high-temperature properties is becoming stronger. Although the strength of transition metal disilicide can meet the requirements of high-temperature materials, its inherent brittleness seriously affects its application. Currently, Xiang et al. [5] investigated the high thermal emission of transition metal disilicides by first-principles calculations. The results showed that the addition of SiO₂ significantly improved their high-temperature properties. Zhang et al. [6] systematically investigated the phase stability, mechanical properties, and electronic structure of transition metal disilicides by first-principles. The results show that for tetragonal TiSi₂, NiSi₂, CuSi₂ and ZnSi₂, they are either thermodynamically unstable or mechanically unstable.

In addition, TMSi₂ (TM = V, Nb and Ta) disilicide is an important member of the transition metal disilicide family. Additionally, many of its properties have been extensively studied, such as magnetization [7], electronic structure and mechanical properties. In previous studies, we could determine that TMSi₂ (TM = V, Nb and Ta) had thermodynamic stability and optical anisotropy [8]. Moreover, VSi₂ is of interest because of its more excellent optical anisotropy [9]. In previous studies, the structural stability, mechanical and thermodynamic properties of VSi₂ have been investigated by experimental and theoretical methods, and it has been shown that VSi₂ has excellent comprehensive properties [10,11,12,13]. It is well known that VSi₂ consists of two crystal structures, namely, the tetragonal C11_b structure and the hexagonal C40 structure. However, for most of the TMSi₂ disilicides, the C40 structure is more stable than the C11_b structure, and the formation energy differences in the C11_b/C40 competition are less than 0.1 eV/atom [14], indicating that the phase stability of C11_b is very close to that of C40. The C40 structure compound has desirable physical properties such as high melting point and low density, which are more attractive in the field of microelectronic devices. The C11_b structure compound has good overall mechanical properties at higher temperatures, but it has brittle fractures at room temperature [15,16]. Therefore, it is necessary to investigate the overall properties of C11_b-VSi₂ in order to reduce its brittle behavior and obtain high strength. Different crystal structures and various defects (including vacancies, dislocations, etc.) combine to determine the elastic properties of the material [17,18]. Therefore, vacancies play a crucial role in the mechanical properties of high-temperature materials. At present, the physical properties of VSi₂ containing vacancy defects are less studied. Therefore, it is essential to study the effect of vacancies on the overall properties of high-temperature alloys [19].

In this article, we study the effect of different atomic vacancies on the mechanical performance of the tetragonal C11_b-VSi₂ by first-principles calculations based on density functional theory (DFT) [20]. This theory has been widely used to predict the electronic, mechanical, and thermal properties of solids [21,22,23,24,25,26,27,28]. On the basis of the VSi₂ vacancy model, we have studied five different vacancies. To check the stability of the vacancies, the vacancy formation energy and phonon dispersion were calculated for the different vacancies. The effects of vacancies on the elastic properties, tough and brittle behavior and electronic properties of the material were investigated in detail, and it was concluded that vacancies can improve the brittle behavior of C11_b-VSi₂.

2. Calculation Method

In this paper, based on density functional theory (DFT), first-principles calculations are performed on the structure, elastic properties, electronic properties and chemical bonds of perfect C11_b-VSi₂ and C11_b-VSi₂ with different vacancies through the CASTEP [29] code. The interactions between the ionic nuclei and valence electrons are described in terms of a ultrasoft pseudopotentials. The electronic configurations of the pseudo-atoms are V 3s²3p⁶3d³4s² and Si 3s²3p², respectively. To verify that the total energy in the ground state is convergent, a plane wave basis is set to an electronic wave function with a cut-off energy of 500 eV. All lattice properties, atomic locations and internal coordinates in the system are entirely loosened during structural optimization. Using the PHONOPY code [30], we calculated the phonon frequencies for ideal C11_b-VSi₂ and VSi₂ including varying vacancies to assess the dynamic stability.

3. Results and Discussion

3.1. Structural Properties for C11_b-VSi₂ with Different Atomic Vacancies

Here, VSi₂ is a tetragonal (C11_b) structure. To investigate the relationship between different vacancies and properties, we constructed a crystal structure containing nine V atoms and ten Si atoms as shown in Figure 1 and discussed C11_b-VSi₂ with different vacancies, such as V-va1, V-va2, Si-va1 and Si-va2 [31]. By discussing the forming energies of perfect C11_b-VSi₂ and different atomic vacancies C11_b-VSi₂, the thermodynamic stability of C11_b-VSi₂ containing vacancies is studied. The vacancy formation energy (E_f), which is the energy necessary to generate a vacancy and is dependent on the relationship between the atom and its nearby atoms whenever the element is removed, is used to determine the thermodynamic stability of vacancies [32]. Only when the formation energy is negative, it is known that a solid material is thermodynamically stable. In particular, the lower the value of the vacancy formation energy, the better the thermodynamic stability [33]. When atomic vacancies are present, C11_b-VSi₂ with different atomic vacancies have different levels of thermodynamic stability. Therefore, to investigate the structural stability of C11_b-VSi₂, we estimated the atomic vacancy formation energy as follows:

$$E_f=E_{Vac}^M-E_{Perf}^{Total}+{\mu }_M$$

(1)

The total energy of C11_b-VSi₂ with an M (M = V and Si) vacancy and a C11_b-VSi₂ perfect crystal without a vacancy are $E_{Vac}^M$ and $E_{Perf}^{Total}$, respectively. μ_M is the chemical potential of the M atom that has been eliminated.

The vacancy formation energies and lattice parameters of perfect C11_b-VSi₂ and C11_b-VSi₂ containing different atomic vacancies are listed in

Table 1. The lattice parameters and vacancy formation energy E_f (eV/atom) of C11_b-VSi₂ with various vacancies.

Phase	a	b	c	Ef	Refs.
VSi2	3.178	3.178	7.7	−0.297	present
-	3.155	3.155	7.729	-	[10]
V-va1	3.134	3.133	7.78	−0.186	present
V-va2	3.134	3.134	7.779	−0.185	present
Si-va1	3.117	3.116	7.851	−0.151	present
Si-va2	3.118	3.118	7.851	−0.152	present

, as well as the other theoretical results [10]. Because the vacancy formation energies in C11_b-VSi₂ are less than zero, these vacancies are thermodynamically stable in the ground state [34]. Moreover, with the introduction of vacancies, the vacancy formation energy of C11_b-VSi₂ containing vacancies is slightly larger compared with that of the perfect C11_b-VSi₂. It shows that the introduction of vacancies increases the thermodynamic instability of C11_b-VSi₂. The calculated vacancy formation energies of V vacancies are smaller than those of Si vacancies, as shown in Table 1, indicating that the V vacancies are more thermodynamically stable than the Si-vacancies. As a result, C11_b-VSi₂ is more inclined to form V vacancies.

Moreover, the lattice parameters of perfect C11_b-VSi₂ and C11_b-VSi₂ containing different vacancies are discussed in depth, which have also been discussed in NbSi₂ [35]. In general, the lattice parameters of C11_b-VSi₂ containing vacancies are slightly smaller than that of perfect C11_b-VSi₂ along the a-axis and the b-axis, but slightly larger than that of perfect C11_b-VSi₂ along the c-axis. The differences in the lattice parameters determine the elastic properties. For V vacancies, the lattice constants of V-va1 are approximately equal to those of V-va2. Additionally, C11_b-VSi₂ with V vacancies have larger values of lattice parameters a and b than Si vacancies; however, the c values are smaller than those of Si vacancies. For Si vacancies, the lattice constant values for Si-va2 vacancies are larger than those for Si-va1 in the a and b directions, and the c values are approximately equal. The position and atomic type of the removed atoms lead to differences in the variation of the lattice constants.

The computed phonon dispersion of perfect C11_b-VSi₂ and C11_b-VSi₂ with vacancies in the Brillouin zone along the high symmetric direction is shown in Figure 2. The phonon spectrum is an important indicator to judge whether the system is stable [36]. If the phonon frequencies are all above the zero point, it means that the material does not have imaginary frequencies, indicating that the material is stable [37]. From Figure 2, the phonon frequencies of both perfect C11_b-VSi₂ and C11_b-VSi₂ in the presence of vacancies have no imaginary frequencies, indicating that they are dynamically stable. This yields results that are consistent with those obtained for the vacancy formation energy.

3.2. Elastic Properties

Before analyzing the elastic properties of C11_b-VSi₂ with various vacancies, it is necessary to consider the mechanical stability of the material [38]. The elastic constants (C_ij) and the elastic flexibility matrix (S_ij) are calculated for perfect C11_b-VSi₂ and C11_b-VSi₂ containing different vacancies in

Table 2. The elastic constants C_ij (GPa) of C11_b-VSi₂ with various vacancies.

Phase	C11	C12	C13	C33	C44	C66
VSi2	198.7	125.7	97.2	390.9	140.9	95.7
Si-va1	224	148.2	92	346.4	155.1	145.9
Si-va2	222.2	149.5	91.2	341.2	152.7	147.5
V-va1	201	159.8	87.2	345.6	161.3	149
V-va2	203.1	155.6	84.8	348	162	149.1

and

Table 3. Calculated elastic compliance constants S_ij (GPa⁻¹) of C11_b-VSi₂.

Phase	S11	S12	S13	S33	S44	S66
VSi2	0.00867	−0.00504	−0.00091	0.00301	0.00710	0.01045
Si-va1	0.00818	−0.00511	−0.00082	0.00332	0.00645	0.00685
Si-va2	0.00845	−0.00536	−0.00082	0.00338	0.00655	0.00678
V-va1	0.01369	−0.01052	−0.00081	0.00329	0.00620	0.00671
V-va2	0.01220	−0.00909	−0.00073	0.00325	0.00617	0.00671

, respectively. The generalized mechanical stability criterion for the tetragonal system is C₁₁ > 0, C₃₃ > 0, C₄₄ > 0, C₆₆ > 0, C₁₁ > C₁₂, C₁₁ + C₃₃ − 2C₁₃ > 0, 2C₁₁ + C₃₃ + 2C₁₂ + 4C₁₃ > 0 [39]. Clearly, the perfect C11_b-VSi₂ and C11_b-VSi₂ containing different vacancies satisfy the stability criterion for all elastic constants. Therefore, they have mechanical stability.

The linear compression resistance along the a-axis, b-axis, and c-axis is represented by C₁₁, C₂₂, and C₃₃, respectively [40]. Table 2 shows that the C₁₁ values for perfect C11_b-VSi₂ and C11_b-VSi₂ containing atomic vacancies are generally lower than the C₃₃ values, indicating that the deformation resistance along the a-axis is lower than the deformation resistance along the c-axis. The C₁₁ values of C11_b-VSi₂ with different vacancies are larger than those of perfect C11_b-VSi₂, yet the C₃₃ values are lower than those of perfect C11_b-VSi₂. It shows that the vacancies created by the removal of atoms increase the resistance to deformation along the a-axis but weaken resistance to deformation along the c-axis. Additionally, since the C₁₁ value of C11_b-VSi₂ containing Si vacancies is larger than that of C11_b-VSi₂ containing V vacancies, the deformation resistance produced by the removal of Si atoms is larger than that produced by the removal of V atoms in the direction along the a-axis. In addition to this, it is well known that C₄₄ and C₆₆ are correlated to shear modulus, with larger values of C₄₄ and C₆₆ corresponding to larger shear modulus. Table 2 shows that the C₄₄ and C₆₆ values of C11_b-VSi₂ with different vacancies are larger than those of perfect C11_b-VSi₂, indicating that the atom vacancies significantly enhance the shear deformation resistance of C11_b-VSi₂. Moreover, the V-va2 vacancies has the largest C₃₃, C₄₄ and C₆₆ values compared to other vacancy types, which implies that C11_b-VSi₂ with V-va2 vacancies have a greater resistance to deformation.

The mechanical properties of C11_b-VSi₂ containing these vacancies are investigated in order to establish the link between vacancies and mechanical performance. The elastic moduli (including bulk modulus B, Poisson’s ratio v, Young’s modulus E, and shear modulus G) of perfect C11_b-VSi₂ and C11_b-VSi₂ with different atomic vacancies for the C11_b structure obtained by the Voigt–Reuss–Hill method are presented in

Table 4. Calculated bulk modulus, B (GPa), shear modulus, G (GPa), Young’s modulus, E (GPa) Poisson’s ratio, v and Vickers hardness, H_V (GPa) of C11_b-VSi₂ with various vacancies.

Phase	BV	BR	BH	GV	GR	GH	E	ν	BH/GH	HV
VSi2	159.0	150.6	154.8	106.6	82.1	94.3	235.2	0.247	1.642	13.1
Si-va1	161.8	160.4	161.1	122.3	90.9	106.7	262.1	0.229	1.510	15.7
Si-va2	161.1	159.8	160.4	121.2	88.9	105	258.7	0.231	1.528	15.3
V-va1	157.3	155.7	156.5	121.9	65.7	93.8	234.5	0.250	1.668	12.8
V-va2	156.0	154.3	155.2	122.7	71.1	96.9	240.7	0.241	1.600	13.7

, where B_V (G_V) and B_R (G_R) are B (G) in the Voigt and Reuss approximations, respectively, and the expressions are as follows [41,42,43,44,45]:

$$B=\frac{B_V+B_R}{2}$$

(2)

$$G=\frac{G_V+G_R}{2}$$

(3)

$$E=\frac{9\mathit{GB}}{G+3B}$$

(4)

$$\nu =\frac{3B-E}{6B}$$

(5)

$$H_V=2{\left(\frac{G^3}{B^2}\right)}^{0.585}-3$$

(6)

The ratio of bulk modulus to shear modulus (B/G) validates the solid material’s brittleness/ductility [46]. If the B/G ratio is higher than 1.75, the solid material will exhibit ductility. Otherwise, the material will exhibit brittleness. The higher the B/G value, the more ductile the solid material is [47]. Based on the B/G ratios of C11_b-VSi₂ with varying vacancies as well as perfect C11_b-VSi₂ shown in Table 4, the B/G of perfect C11_b-VSi₂ is 1.642, which is less than 1.75, showing that C11_b-VSi₂ is brittle. Meanwhile, with the introduction of vacancies, the calculated B/G values of C11_b-VSi₂ containing these vacancies generally decrease, except the B/G value of V_-va1 (1.668). Notably, the B/G values of C11_b-VSi₂ with V vacancies are all larger than those of C11_b-VSi₂ with Si vacancies, and thus, C11_b-VSi₂ containing V vacancies have less brittleness.

Bulk modulus is a macroscopic property of a material that reflects its resistance to external compression, i.e., incompressibility, as is widely known. As seen in Table 4, the bulk modulus of C11_b-VSi₂ containing different vacancies is larger than that of perfect C11_b-VSi₂, indicating that the removal of atoms instead makes the material more incompressible. The force that resists form change under shear stress is known as the shear modulus (G). Shear stress is more closely related to hardness than bulk modulus, and shear modulus is a more appropriate predictor of hardness [48]. In Table 4, C11_b-VSi₂ with the introduction of V-va1 vacancy has the smallest shear modulus and is smaller than perfect C11_b-VSi₂. Therefore, its ductility is slightly better than that of perfect C11_b-VSi₂. Additionally, all the remaining types of vacancies increase the shear modulus, compared to perfect C11_b-VSi₂. Thus, the hardness values of C11_b-VSi₂ with V-va2, Si-va1, and Si-va2 increase, and this variation of hardness with vacancy has been found in chromium silicide [49]. In addition, the stiffness of a solid can be described using Young’s modulus. The stiffness of a solid is proportional to its Young’s modulus. As can be seen from Table 4, the introduction of atomic vacancies increases the Young’s modulus of C11_b-VSi₂, with the largest being for Si-va1 vacancies and only those containing V-va1 vacancies having lower Young’s modulus. Accordingly, removing atoms increases the elastic stiffness of VSi₂. It is congruent with the results of the elastic constants. From the point of view of interatomic interactions, elastic modulus is closely related to the interaction between atoms. The increased elastic modulus of C11_b-VSi₂ containing different vacancies, except for the V-va1 vacancy, can be tentatively judged that the removal of atoms enhances the atomic interaction of C11_b-VSi₂.

Poisson’s ratio is a well-known method for determining a solid’s stability under shear deformation [50]. Solids having a Poisson’s ratio of −1 to 0.5 are relatively stable under shear deformation in general. Furthermore, a larger Poisson’s ratio indicates that the solid is more malleable. As shown in Table 4, the Poisson’s ratios of perfect C11_b-VSi₂ and C11_b-VSi₂ containing different vacancies is in the range of 0.229~0.25, which is between −1 and 0.5, indicating that they are stable solids. Among them, VSi₂ containing V-va1 vacancies has the largest Poisson’s ratio, indicating that it has better plasticity than the other compounds. It also corroborates that the Vickers hardness of VSi₂ containing V-va1 vacancies is the smallest among them. Poisson’s ratio also can be used to estimate the solid’s brittleness and ductility. If the solid has a ν > 0.33, the material is ductile; otherwise, it is fragile. From the calculations in Table 4, ν < 0.33 for perfect C11_b-VSi₂ and C11_b-VSi₂ with different vacancies; thus, they are brittle. This result agrees well with that of B/G.

3.3. Elastic Anisotropy

The elastic anisotropy is an important index reflecting the mechanical anisotropy of materials and plays a very large role in the generation of microcracks. To better describe the elastic anisotropy, in this paper, we use elastic anisotropy indices such as the universal anisotropy index (A^U), compression and shear anisotropy percentages (A_comp and A_shear), and shear anisotropy factors (A₁, A₂, and A₃) to investigate the elastic anisotropy [51]. The calculation equations are as follows:

$$A^{\mathrm{U}}=\frac{B_V}{B_R}+5\frac{G_V}{G_R}-6$$

(7)

$$A_{\mathrm{comp}}=\frac{B_V-B_R}{B_V{+B}_R}\times 100\%$$

(8)

$$A_{\mathrm{shear}}=\frac{G_V-G_R}{G_V{+G}_R}\times 100\%$$

(9)

Furthermore, due to the fact that the VSi₂ we study in this paper is a tetragonal crystal structure, considering its shear anisotropy factor:

$$A_1=\frac{{4C}_{44}}{C_{11}+C_{33}-2C_{13}}$$

(10)

$$A_2=\frac{{4C}_{55}}{C_{22}+C_{33}-{2C}_{23}}$$

(11)

$$A_3=\frac{4C_{66}}{C_{11}+C_{22}-2C_{12}}$$

(12)

A₁, A₂, and A₃ represent the degree of shear anisotropy corresponding to the (100), (010) and (001) planes, respectively. For these elastic anisotropy indices, if A₁ = A₂ = A₃ = 1 and A^U = A_comp = A_shear = 0 [52], the solid shows elastic isotropy; otherwise, it is anisotropic. At the same time, when the solid has larger values of A^U, A_comp and A_shear, it has higher elastic anisotropy.

The calculated elastic anisotropic indices are listed in

Table 5. Calculated elastic anisotropy indeces of C11_b-VSi₂.

Phase	AU	Acomp	Ashear	A1	A2	A3
VSi2	1.5479	0.0271	0.1298	1.4261	1.4305	2.6219
Si-va1	1.7359	0.0043	0.1473	1.6056	1.6255	3.8829
Si-va2	1.8248	0.0040	0.1537	1.6032	1.6359	4.0859
V-va1	4.2873	0.0051	0.2996	1.7335	1.7281	7.2068
V-va2	3.6397	0.0055	0.2663	1.6986	1.6894	6.3582

. From Table 5, both perfect C11_b-VSi₂ and C11_b-VSi₂ with different vacancies are anisotropic, since their elastic anisotropy indices deviate from 0. The greater the deviation, the greater the anisotropy [53]. It is obvious that the changes corresponding to A_comp and A_shear are different. A_shear increases with the introduction of vacancies. However, A_comp is decreasing. This may be due to the difference in the variation of bulk and shear moduli. Therefore, using A_comp and A_shear alone to evaluate the material elastic anisotropy has limitations. On the contrary, A^U is more accurate to evaluate the elastic anisotropy by considering both bulk modulus and shear modulus. As can be seen from Table 5, the A^U values of C11_b-VSi₂ containing different vacancies are larger than those of VSi₂, and in contrast, the A^U values of C11_b-VSi₂₂ containing V vacancies increase sharply. It indicates that the introduction of vacancies improves the anisotropy of the material, especially the V vacancies. Among them, the elastic anisotropy of vacancy V-va1 is the largest. In addition, the shear anisotropy indices A₁, A₂ and A₃ are not equal to 1, indicating that they are all anisotropic. The introduction of the vacancy increases A₁, A₂ and A₃, meaning that the vacancy increases the shear anisotropy of C11_b-VSi₂ in the (100), (010) and (001) plane all. For A₃, the increase is the largest. It shows that the vacancy enhances the shear anisotropy of C11_b-VSi₂ in the (001) plane. the increase in V-va1 is particularly significant.

The elastic anisotropy of the crystal can also be visualized by the three-dimension (3D) surface construction diagram of the elastic modulus. When the crystal is isotropic, the 3D surface construction diagram is perfectly spherical. On the contrary, it is anisotropic. In this work, we focus on the elastic anisotropy of the bulk modulus B and Young’s modulus E through the 3D surface construction diagram. For the tetragonal structure, the inverse of the bulk modulus B and Young’s modulus E are calculated as follows [54]:

$$\frac{1}{B}=\left(S_{11}{+S}_{12}{+S}_{13}\right){-(S}_{11}{+S}_{12}{-S}_{13}{-S}_{33}{)l}_3^2$$

(13)

$$\frac{1}{E}=S_{11}\left(l_1^4+l_2^4\right)+\left(2S_{13}+S_{44}\right)\left(l_1^2{l}_3^2+l_2^2{l}_3^2\right)+S_{33}{l}_3^4+\left(2S_{12}+S_{66}\right)l_1^2{l}_2^2$$

(14)

Here, S_ij is the elastic compliance constant listed in Table 3, and l₁, l₂, and l₃ are the direction cosines.

Figure 3 and Figure 4 show the 3D surface construction diagrams of the bulk modulus and Young’s modulus, respectively. When the deviation of the 3D surface construction diagram from the spherical shape is greater, the anisotropy of the solid is greater [55]. As can be seen from Figure 3, there is no significant difference between the 3D diagrams of bulk modulus of C11_b-VSi₂ containing vacancies, then it is necessary to consider the 3D diagram of Young’s modulus. It is obvious from Figure 4 that the 3D diagrams of Young’s modulus of C11_b-VSi₂ containing different vacancies and perfect C11_b-VSi₂ are non-spherical. Compared with the 3D diagrams of bulk modulus, the 3D diagrams of Young’s modulus are more irregular in shape, indicating that they are all anisotropic, and the anisotropic feature is greater than that of the bulk modulus. Moreover, with the introduction of vacancies, the deviation of 3D diagrams of Young’s modulus from the spherical shape is greater, indicating that the different atomic vacancies enhance the anisotropy of Young’s modulus. This result is in good agreement with the results corresponding to A^U in Table 5.

However, the 3D construction diagrams in Figure 4 do not clearly show the subtle differences in the elastic modulus anisotropy of C11_b-VSi₂ containing different vacancies. Therefore, to see more details of the elastic anisotropy, the two-dimension (2D) projections of bulk modulus and Young’s modulus in the (001) and (100) planes, which have been employed successfully [56], for both perfect C11_b-VSi₂ and C11_b-VSi₂ containing different atomic vacancies are shown in Figure 5. From Figure 5, in general, the introduction of vacancies increases the anisotropy of the elastic modulus. For the Young’s modulus, Figure 5a shows that the shape of C11_b-VSi₂ containing different atomic vacancies in the (001) plane is more deviated from circular than that of perfect C11_b-VSi₂. It means that the vacancies make Young’s modulus more anisotropic. As shown in Figure 5c, the 2D projection of C11_b-VSi₂ containing Si atomic vacancies in the (100) plane is approximately or slightly anisotropic to that of perfect C11_b-VSi₂. In contrast, C11_b-VSi₂ containing V-atom vacancies is more obviously irregular in shape and more anisotropic. This result is in agreement with the one corresponding to A^U. For the bulk modulus, in the (001) plane, the shapes of the graphs are all circular and the crystals show no significant anisotropy. However, in the (100) plane, the shape of C11_b-VSi₂ containing different vacancies deviates significantly from a circle and shows a higher anisotropy.

Table 6. Calculated directional elastic moduli in three principal directions (in GPa) of C11_b-VSi₂.

Phase	E			B
	[100]	[010]	[001]	[100]	[010]	[001]
VSi2	115.3	115.3	332.2	121.3	121.3	277.3
Si-va1	122.2	122.2	301.2	146.7	146.7	196.4
Si-va2	118.3	118.3	295.9	145.4	145.4	189.7
V-va1	73.0	73.0	303.9	139.8	139.8	197.6
V-va2	82.0	82.0	307.7	138.7	138.7	184.4

shows the elastic moduli in the [100], [010] and [001] directions obtained according to Figure 5. It is obvious from Table 6 that the values of Young’s modulus and bulk modulus in the [100] and [010] directions are smaller than those in the [001] direction because the value of C₁₁ is smaller than that of C₃₃, making the crystal easier to compress along the a-axis [57].

3.4. Electronic Properties

As mentioned above, the reason for the change in the elastic modulus of C11_b-VSi₂ containing atomic vacancies is related to the inter-charge interactions, where the type and position of the atoms determine the inter-charge interactions [58]. The removal of V and Si atoms in C11_b-VSi₂ due to inter-charge interactions changes the electron equilibrium concentration between adjacent atoms and alters the chemical bonding in VSi₂. There are two different types of chemical bonds in perfect C11_b-VSi₂, the V-Si bond and the Si-Si covalent bond. The bond lengths of its Si-Si covalent and V-Si bonds were calculated to be 2.529 Å and 2.579 Å, respectively, which is consistent with earlier theoretical predictions. The bond lengths of Si-Si covalent and V-Si bonds for perfect C11_b-VSi₂ and C11_b-VSi₂ with different vacancies are shown in

Table 7. The bond lengths of Si-Si covalent bonds and V-Si bonds for perfect C11_b-VSi₂ and C11_b-VSi₂ with different atomic vacancies.

Lengths (Å)	VSi2	Si-va1	Si-va2	V-va1	V-va2
Si-Si	2.529	2.555	2.552	2.472	2.458
V-Si	2.578	2.530	2.531	2.528	2.526

. As shown in Table 7, the insertion of Si vacancies boosts the charge interaction between the V atom and the Si atom while weakening the charge interaction between the Si atom and the other Si atom. As a result, the single cell cohesion energy is weakened due to electron collapse, leading to lattice contraction. This is the major reason why Si vacancies have a greater elastic modulus than perfect C11_b-VSi₂.

For the V-atom vacancy, the presence of the V atom vacancy shortens the bond length of both the V-Si bond and the Si-Si covalent bond. In other words, the removal of the V atom enhances the interatomic charge interactions. We suggest that the loss of the V-Si bond on the shear surface causes the enhancement of the elastic modulus.

The electron density difference can be sued to further investigate electron transfer in chemical bonding [59]. Figure 6 shows the electron density difference for perfect C11_b-VSi₂ and C11_b-VSi₂ with V and Si atomic vacancies. In each graph, the electron density difference ranges from −0.1e/ų to 0.1e/ų. The blue color implies the maximum localization of electrons in this region, and the red color indicates the maximum delocalization of electrons in this region. In Figure 6a, intact V-Si bonds and Si-Si covalent bonds in perfect C11_b-VSi₂ are observed. The removal of atoms clearly disrupts the localized hybridization between adjacent atoms compared to the figure with atomic vacancies. By comparison, we find that the electron delocalization of V vacancies and Si vacancies is weaker than that of perfect C11_b-VSi₂. In other words, the localized hybridization between V and Si atoms is improved by the removed atoms, which has been confirmed in TaSi₂ that the removals of Ta and Si can enhance the electron hybridization between Ta and Si [60]. As can be seen in Figure 6, for C11_b-VSi₂ containing vacancies, the inter-atomic charge interactions are stronger than the corresponding charge interactions for perfect C11_b-VSi₂. This is the reason that the elastic modulus of C11_b-VSi₂ containing vacancies increases.

In Figure 6b,c, for Si vacancies, the removal of Si atoms enhances the local hybridization between V and Si atoms, but the Si-Si atom interactions are weakened. This is consistent with the results obtained from the discussion of bond lengths above. Thus, the enhanced elastic modulus comes from the electronic bonding properties. From Figure 6d,e, the removal of V atoms when V vacancies are present also results in enhanced Si-Si inter-atomic interactions.

We calculated the total and partial densities of states (DOS and PDOS) to discuss the bonding properties of perfect C11_b-VSi₂ and C11_b-VSi₂ containing different atomic vacancies in Figure 7. DOS denotes the number of electron states per unit energy interval when the electron energy levels are quasi-continuously distributed [61,62]. The dashed line with zero energy represents the Fermi energy level [63]. The valley at the Fermi energy level is called the pseudogap, which indicates the presence of a stable phase. From the comparison in Figure 7, the pseudogap of C11_b-VSi₂ with V vacancies is the smallest, proving that V vacancies are more stable than Si vacancies. This result is in line with the vacancy formation energy conclusion. The TDOS profile of C11_b-VSi₂ around the Fermi energy level is predominantly from the V-3d state and the Si-3p state, showing substantial hybridization between V and Si atoms, as seen in Figure 7. Figure 7a shows that the Si-3p state is divided into two parts by the Si-3s state, indicating the formation of hybridization between Si and Si atoms in VSi₂. However, the PDOS profile resulting from the removal of V and Si atoms is slightly different from that of the perfect C11_b-VSi₂. The energy of the V vacancies is determined by a tiny change around the Fermi energy level. As can be seen in Figure 7b,c, for the V vacancy, the PDOS produces some small peaks near the Fermi energy level. The charge energy of V and Si atoms increases compared to that of perfect C11_b-VSi₂. The introduction of Si vacancies has a similar feature. The DOS profile of C11_b-VSi₂ with Si vacancies is similar to that of C11_b-VSi₂ with V vacancies, while more electrons are transferred from the lower energy region to the Fermi energy level compared with the density of states of V vacancies. It leads to stronger local hybridization between the Si-3s and Si-3p states, forming Si-Si covalent bonds. The loss of Si atoms, in particular, causes a charge transfer from the Si-3s state to the Si-3p state.

3.5. Debye Temperature

The Debye temperature (θ_D) is a fundamental thermodynamic parameter that is affected by melting point, coefficient of thermal expansion, specific heat, and other properties, and describes the lattice vibrations and changes in specific heat [64]. In order to assess the overall performance of high-temperature alloys, the Debye temperatures of perfect C11_b-VSi₂ and C11_b-VSi₂ with different atomic vacancies need to be considered [65]. Under the Debye model, the Debye temperature of a solid is given by the modulus of elasticity, which is calculated from the sound velocity, and the expression for the Debye temperature is as follows:

$${\theta }_{\mathrm{D}}=\frac{h}{K_{\mathrm{B}}}{\left[\frac{3n}{4\pi }\left(\frac{N_{\mathrm{A}}\rho }{M}\right)\right]}^{\frac{1}{3}}{v}_{\mathrm{m}}$$

(15)

In the formula, k_B is Boltzmann’s constant, n is the total number of atoms in the molecule, N_A is Avogadro’s constant, h is Planck’s constant, ρ is the density, and M is the molecular weight.

The longitudinal and transverse sound velocities are denoted by ν_l and ν_t, respectively, and the formulas are as follows:

$$v_{\mathrm{l}}={\left[\left(B+\frac{4G}{3}\right)/\rho \right]}^{\frac{1}{2}}$$

(16)

$$v_{\mathrm{t}}={\left(\frac{G}{\rho }\right)}^{\frac{1}{2}}$$

(17)

The mean sound velocity can be calculated from ν_l and ν_t with the following equations:

$$v_{\mathrm{m}}={\left[\frac{1}{3}\left(\frac{2}{v_{\mathrm{t}}^3}+\frac{1}{v_{\mathrm{l}}^3}\right)\right]}^{-\frac{1}{3}}$$

(18)

Table 8. The density ρ, sound velocities (longitudinal ν_l, transverse ν_t and mean ν_m) and Debye temperature θ_D.

Phase	ρ (g/cm3)	νl (m/s)	νt (m/s)	νm (m/s)	θD (K)
VSi2	4.561	7842.4	4546.8	3570.6	452
V-va1	4.379	8018.8	4628.3	3635.1	457
V-va2	4.378	8061.2	4706.8	3695.3	465
Si-va1	4.513	8198.9	4862.5	3815.6	480
Si-va2	4.509	8161.8	4825.3	3786.9	476

lists the density, longitudinal sound velocity, transverse sound velocity, mean sound velocity and Debye temperature for perfect C11_b-VSi₂ and C11_b-VSi₂ with different vacancies. It is known that the properties of the chemical bonds determine the Debye temperature. The higher the Debye temperature, the higher the interatomic forces and the higher the chemical bond strength. As can be seen from Table 8, the Debye temperature containing V vacancies and Si vacancies is larger than that of perfect C11_b-VSi₂, indicating that the introduction of vacancies enhances the chemical bond strength. Among them, the highest Debye temperature for the Si-va1 vacancy results in the strongest chemical bond strength. This is consistent with the conclusions obtained in the previous section.

4. Conclusions

In summary, we used first-principle calculations to study the structural properties, thermal stability, elastic properties, electronic properties and Debye temperature of perfect C11_b-VSi₂ and C11_b-VSi₂ with different vacancies and have concluded the following:

(1) The introduction of V vacancies and Si vacancies in C11_b-VSi₂ is thermodynamically stable at the ground state. Moreover, the V vacancies are more thermodynamically stable than the Si vacancies in C11_b-VSi₂.

(2) The introduction of vacancies enhanced the bulk modulus, shear modulus and Young’s modulus of C11_b-VSi₂, which significantly improved the mechanical behavior of C11_b-VSi₂.

(3) The elastic anisotropy results indicate that the introduction of vacancies enhances the elastic anisotropy of C11_b-VSi_2, and the V-atom vacancies are the most pronounced.

(4) The difference in electron density difference and density of states between perfect C11_b-VSi₂ and C11_b-VSi₂ with different vacancies suggests that the introduction of vacancies enhances the interactions between charges.

(5) Calculations on the Debye temperature show that the introduction of vacancies enhances the chemical bond strength of C11_b-VSi₂.